Laetitia Le Pourhiet

Benjamin Huet

Dave A May

Loic Labrousse

Laurent Jolivet

Kinematic interpretation of the 3D shapes of metamorphic core complexes

  L'archive ouverte pluridisciplinaire

Introduction

Metamorphic Core Complexes (MCCs) form dome shaped structures of metamorphic rocks exhumed from the middle to lower crustal depth under a detachment fault. These domes are present in many post-orogenic contexts. The dome axis, which is recognized by mapping the schistosity pattern in the exhumed metamorphic rocks, can either be elongated parallel to the stretching lineation (a-type domes Figure 1, orange structure), or normal to this direction (b-type domes Figure 1, blue structure). In the first case, the schistosity strikes mainly parallel to the direction of stretching while in the second case it strikes normal to it. Hence, b-type domes are cylindrical structures and may be modeled in 2D plane strain cross-section while the a-type domes are intrinsically 3D structures since local horizontal constriction occurs normal to stretching. This a-type versus b-type classification, first defined in [START_REF] Jolivet | Strain localization during crustal-scale boudinage to form extensional metamorphic domes in the Aegean Sea[END_REF]

based on

Mediterranean examples, relies on a quantitative description of the finite-strain field. The interpretation of these structures in terms of larger scale kinematics is an ill posed inverse problem, which reaches by definition, a non-unique answer.

Structural geologists generally unambiguously interpret b-type domes as the mark of extension. This interpretation could be questioned, however, since the structure is cylindrical, it constitutes the simplest conceptual model. The Bitterroot Eocene core complex in the Northern Rockies, Idaho [START_REF] Hyndman | Bitterroot dome -Sapphire tectonic block, an example of a plutonic core-gneiss-dome complex with its detached suprastructure. Cordilleran Metamorphic Core Complexes[END_REF] and the North Cyclades Detachment System (NCDS), Greece [START_REF] Jolivet | The North Cycladic Detachment System[END_REF] constitutes a set of excellent natural examples of b-type structures, although in both cases, strike-slip faults are inferred in the vicinity of these domes [START_REF] Foster | Kinematics and timing of exhumation of Eocene metamorphic core complexes along the Lewis and Clark fault zone, northern Rocky Mountains, USA[END_REF][START_REF] Philippon | Deciphering subduction from exhumation in the segmented Cycladic Blueschist Unit (Central Aegean, Greece)[END_REF]].

However, the kinematic interpretation of a-type domes is not straightforward and depends on the geodynamic context. In some cases like Kesebir-Kardamos, Bulgaria [START_REF] Bonev | Mesozoic--Tertiary structural evolution of an extensional gneiss dome-the Kesebir--Kardamos dome, eastern Rhodope (Bulgaria--Greece)[END_REF], Buckskin Rawhide, Harcuvar; Arizona, US [Reynolds and Spencer,1985;[START_REF] Howard | Crustal extension along a rooted system of imbricate low-angle faults: Colorado River extensional corridor, California and Arizona[END_REF], Naxos, Greece [START_REF] Kruckenberg | Flow of partially molten crust and the internal dynamics of a migmatite dome, Naxos, Greece[END_REF] or the Western Gneiss Region, Norway [START_REF] Labrousse | Pressure-temperature-time deformation history of the exhumation of ultra-high pressure rocks in the Western Gneiss Region, Norway[END_REF]], the domes are mainly associated to extension with a slight component of strike-slip in the kinematics. In other cases, these domes are principally associated with strikeslip kinematics with a component of transtension like the Dai Nui Con Voi dome, Vietnam, [START_REF] Jolivet | Oligo-Miocene midcrustal subhorizontal shear zone in Indochina[END_REF] or a component of transpression like along the Red River Fault system [START_REF] Leloup | The Ailao Shan-Red River shear zone (Yunnan, China), Tertiary transform boundary of Indochina[END_REF]. However, in the Betics Cordillera [Martinez-Martinez et al. 2002;[START_REF] Augier | P-T-D-t retrograde evolution of the Nevado-Filabride complex (SE Spain) : new insights from in situ 40Ar/39 ages and metamorphic petrology[END_REF] or in the Variscan orogeny (e.g. Montagne noire, France, [Echtler et Malavielle, 1990]), the large-scale kinematics that accompanies their formation is deeply debated. Whether these structures are interpreted as the result of shortening followed by perpendicular stretching, of coeval stretching and lateral shortening, or of pure strike-slip wrenching changes the interpretation of these structures in term of the large scale kinematics at the time of their formation.

Since [START_REF] Buck | Modes of continental lithospheric extension[END_REF], MCCs are known to form in a local extensional setting when the viscous diffusivity of the flow in the lower crust is small, i.e. when lower crust is thick and/or weak. Several processes leading to either, or both, of these conditions have been proposed and modeled: (i) thermal equilibration of a thickened crust at high Moho temperatures [START_REF] Gaudemer | Thermal control on post-orogenic extension in collision belts[END_REF][START_REF] Block | Core complex geometries and regional scale flow in the lower crust: Tectonics[END_REF][START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF][START_REF] Tirel | Sequential development of interfering metamorphic core complexes: Numerical experiments and comparison with the Cyclades, Greece, in: Collision and Collapse at the Africa-Arabia-Eurasia Subduction Zone[END_REF], (ii) partial melting of the lower crust [START_REF] Lister | Plutonism and the origin of metamorphic core complexes[END_REF][START_REF] Vanderhaeghe | Partial melting and flow of orogens[END_REF][START_REF] Rey | Extension rates, crustal melting, and core complex dynamics[END_REF],

(iii) adjunction of water [START_REF] Jackson | Conditions for flow in the continental crust[END_REF] or (iv) underthrusting of weaker units below stronger ones prior to extension [Huet et al. 2011a[Huet et al. , 2011b]]. All the models stated above only accounted for cylindrical structures and therefore are theoretically only applicable to b-type domes.

Here, we first attempt to understand the impact of 3D kinematic extensional boundary conditions on the 3D geometry of MCCs, i.e. the relative geometry of the schistosity pattern versus the stretching lineation pattern within the exhumed crust. In the second part of the paper, we compare the results of the model to the MCCs of the Cyclades, which provides examples of the two types of domes formed during the same, well constrained, geodynamic event. Lastly we draw more general conclusions about the kinematic implications of the observation of a-type domes.

Modeling approach

Choice of initial conditions

The models have a fixed sized geometry of 200 x 200 x 100 km and evolve in time for 12 Myr, with a constant extension rate of 1cm/yr. Given the assumed size and the grid resolution (3 km) of the models, we expect to resolve the first order impact of large scale, 3D kinematic constraints. The models consist of three horizontal layers (See Figure 2 for the vertical profile).

Since the thermo-mechanical initial conditions for forming MCCs are well understood from 2D modeling, and our study focus on the impact of 3D boundary conditions on their resulting shape, the model design accounts for most of the factors that are known to favor the occurrence of MCCs. This includes an initially thickened crust of 50km [START_REF] Buck | Modes of continental lithospheric extension[END_REF][START_REF] Gaudemer | Thermal control on post-orogenic extension in collision belts[END_REF][START_REF] Block | Core complex geometries and regional scale flow in the lower crust: Tectonics[END_REF]] and an initial thermal gradient which is set to 17.5°C/km. We note that this gradient yields a Moho temperature of 875°C, which is higher than the 800°C limit proposed by [START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF]. At asthenospheric depth, the initial temperature is clamped not to exceed 1300°C (Figure 2). The crust itself consists of two layers of 25 km each. At a given temperature, the top layer (brown shades, Figure 2) is mechanically stronger than the lower one (blue shades, Figure 2), imposing a step in the strength profile at the interface between the two layers. These models thus resemble the set up of the pioneering analogue models of MCCs [START_REF] Brun | Analogue modeling of detachment fault systems and core complexes[END_REF], except they account for the re-strengthening of material as it gets exhumed and cools.

The initial weakness of the lower layer can be equally interpreted as the result of adjunction of water in the lower crust [START_REF] Jackson | Conditions for flow in the continental crust[END_REF], as the occurrence of partial melt in the lower crust [START_REF] Lister | Plutonism and the origin of metamorphic core complexes[END_REF][START_REF] Rey | Extension rates, crustal melting, and core complex dynamics[END_REF] or as the presence of underplated soft sediments or continental upper crust [Huet et al. 2011a]. We do not dismiss that the formation of MCCs can be of diapiric nature [START_REF] Lister | Plutonism and the origin of metamorphic core complexes[END_REF], however, in order to keep the model general, we assign a constant density to the whole crust of 2800 kg.m -3 so that buoyancy forces are limited to those arising from the thermal expansion coefficient which is set constant to 3x10 -5 K -1 .

Treatment of the rheology

At the scale of these models, the rheology in any given volume has no reason to be one of pure quartz, pure olivine or plagioclase and must be influenced by the layering of the rocks, their anisotropy, structural softening and hardening due to small scale boudinage and folding. In practice it is not yet possible to account for all this complexity in crustal or lithospheric scale modeling. As the viscous strength of rock depends to first order on temperature, we approximate it by a Newtonian flow rule based on Frank-Kamenetskii, for which we need to provide two parameters, a reference viscosity 0 and a characteristic temperature , in order to compute the effective viscosity,

eff = 0 e T (1)
At high stress, this viscosity is lowered by the action of a Drucker Prager visco-plastic flow rule such that

eff = 1 2 y ˙ II , (2)
where

˙ II = 1 2 ˙ ij ˙
ij is the second invariant of the strain rate tensor ( ˙ ij ) and the yield stress,

, depends on pressure P and plastic strain accumulated by the particles following

y = P sin( p ) + C 0 cos( p ) . (3)
Accounting for plastic strain softening is necessary [START_REF] Lemiale | Shear banding analysis of plastic models formulated for incompressible viscous flows[END_REF] to localize plastic shear bands in visco-plastic codes such as GALE [START_REF] Moresi | A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials[END_REF]. We impose softening on friction in the upper and lower crust with a friction angle decreasing from 30° to 10° with plastic strain varying from 0 to 20%, given the element size (3km), softening is therefore achieved for a displacement of 500 m in nature.

The rheological parameters used for the study are listed in Table 1 and were chosen because they resulted in a yield strength envelope similar to pre-existing 2D models of Huet et al. [2011a], thereby allowing us to verify (or benchmark) our results with pre-existing models run with a other numerical code, namely FLAMAR [START_REF] Yamato | Burial and exhumation in a subduction wedge: Mutual constraints from thermomechanical modeling and natural P-T-t data (Schistes Lustrés, western Alps)[END_REF]]. In the parameterization, we have chosen to keep 0 constant and to vary the characteristic temperature 1/. This choice was governed by the idea that when the temperature is zero at the surface, both the lower and upper crust must have the same effective viscosity in order to enable yielding.

However, with this choice, the viscosity of the lower crust may drop to a non-realistic value at the base of the lower crust. Therefore we adopted a lower cut-off for the viscosity in the models, which we set to 10 19 Pa.s, as in previously published 2D models [Huet et al. 2011a[START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF]].

The previous study by Huet et al. [2011a] investigated the role of the rheological layering and temperature on the distribution, or localization of the deformation without imposing a priori the location of faults. To avoid strain localization on the side of the box, they used random damage all through the upper crust except in the vicinity 30 km from the boundary. The model domain was longer and show that with a 25 km thick brittle crust, strain started to localize by distributed necking with a wavelength of 70 km. Depending on the rheological layering of the crust, either the strain remained distributed and no exhumation of lower crust was occurring or only one of the initial grabens was still active after 2 Myr before turning into a MCC. In the present study, we investigate the role of 3D kinematics on the shape of MCC. We have therefore chosen a rheological layering which favor their formation and imposed a 50 km wide zone with random pre-existing damage in the center of the model to avoid unwanted border effects.

Numerics

All the numerical experiments were run with GALE 1.6.1, an open source code, which solves incompressible Stokes flow coupled with heat conservation via a finite element, particlein-cell method in three dimensions [START_REF] Moresi | A Lagrangian integration point finite element method for large deformation modeling of viscoelastic geomaterials[END_REF]. The computational mesh consisted of 64x64x32 Q1 (trilinear) elements. Evolution in time is obtained through advection of particles (tracking lithology) and the thermal energy equation, with time steps limited to 10% of Courant criterion.

As GALE uses Q1Q1 elements, it is necessary to scale the viscosity of the model to be close to one in order to obtain accurate results and convergence. The scaling used to run the code is provided in Table 2. Additional passive markers have been included to compute the finite strain field in the lower crust and to track P-T paths. An Octave script used for post-processing the finite strain ellipsoid and exporting it into a VTK format for visualization purpose are included as supplementary material. The method used to compute the finite strain tensor in 3D from passive markers is detailed in appendix A.

Each run consists of 200-300 time steps, requiring approximately 72 hr of computation using a parallel geometric multi-grid (GMG) solver on 16 CPU (dual core AMD Opteron 2.6

GHz, each with 1 GB of RAM). A full description of the multi-grid method used in GALE is provided in Appendix B. We refer the reader, who would like to reproduce the results or reuse the input files, to the supplementary material. It includes the input files of the models, solver options and patches for GALE to reproduce the specific boundary conditions and viscous cutoff implemented for this study.

Impact of boundary conditions on the shape of MCCs

Three types of extensional boundary conditions are considered. The boundary conditions are represented by black arrows that indicate the velocity applied on the front (pink) and back (cyan) boundaries of the models (Figure 3a). First, 3D cylindrical extension is applied by imposing a constant normal velocity on one side of the model (the cyan boundary, Figure 3a and4). This model is used both for verifying the modeling approach as compared to existing 2D 10 models and to compare the results with a non-cylindrical, 3D boundary conditions. The second model considers extension occurring at an extensional step-over between two left lateral strikeslip faults applied on the front (pink) and back (cyan) boundaries of the model (Figure 3b and5).

The third model considers the case of a transtensional fault propagator which accommodates a left-lateral step in the rate of extension at the front of the model (pink) (Figure 3c and6). The two lateral boundaries are free-slip boundaries.

Surface deformation

The surface deformation pattern after 12 Myr of extension is outlined in map view, by the light brown stripes, which were originally forming square (Figure 3). The most striking feature at that point is that the cylindrical and the fault propagator kinematics lead to the exhumation of deep crustal material to the surface (blue material on Figure 3a, 3c), whilst the step-over kinematics model does not (Figure 3b). In the case of cylindrical boundary conditions (Figure 3a), the model produces one structure, elongated normal to the direction of stretching. In the case of fault propagator kinematics, the initially deep crustal material is exhumed along a curvilinear trend which follows the fault propagator and turns to become orthogonal to the free slip side of the model. This trend is therefore parallel to the stretching direction close to the back boundary, whereas it is perpendicular to the stretching direction close to the left hand side boundary.

Looking in more detail, one sees that in between the two branches, there are less exhumed rocks and that the strike-slip part appears to be elongated further towards the back of the model than the location of the cylindrical part of the exhumed structure (Figure 3c).

Topography

Colored isolines indicate the surface topography (Figure 3). All the models produce relatively low topographic expression with an altitude ranging from -800 m to 800 m for the model 1 (Figure 3a), from -1200m to 1000m for model 2 (Figure 3b) and from -1200 m to 1300 m for model 3 (Figure 3c). The exhumed deeper crustal rocks tend to be located below topographic lows rather than beneath the topographic highs. Noticeably, the two models with a strike-slip component produce significantly more topographic expression than the cylindrical model. The steepest gradients are observed along the strike-slip faults. Along the front side, the topography may be over estimated due to border effects, but in the center of the model, the topographic gradients remains significant. In the regions exhibiting strike-slip deformation, the topographic lows tend to be aligned in the direction of stretching. Comparing the cylindrical model (Figure 3a) with the cylindrical part of the fault propagator model (Figure 3c), the presence of the strike-slip fault changes the geometry of topography, since a secondary topographic low that is absent from the purely cylindrical model appears.

3D-geometry and finite-strain inside the domes

Along cross-sections through all the models, one sees that the Moho remains strikingly flat (Figure 4, 5 and 6) while the deeper, ductile, part of the crust is exhumed within dome shaped structures, even in model 2 where it does not reach the surface. The finite strain field calculated locally in the lower crust is represented by the stretching lineation superposed to the geometry of the models (cylinders colored in function of the strike of the lineation). The lineation and the pole of the schistosity are also plotted on stereo-diagrams.

In the case of cylindrical extension (Figure 4), the dome is approximately symmetric and elongated normal to the stretching direction. The detachment fault at the top of the lower crust is outlined by nearly horizontal lineations parallel to the stretching direction (Figure 4a and red dots on Figure 4b), whilst vertical lineations delineate the deep core of the dome. The poles of the schistosity are aligned on a vertical great circle parallel to the stretching direction (black dots on Figure 4b). This pattern reflects the cylindrical nature of the dome. However, a few vertical schistosity planes that strike parallel to the stretching direction mark local constriction. The similarity between this cylindrical model and the 2D models on which we have based the model set up [Huet et al., 2011a] provides confirmation that although the models were run with a lower resolution in 3D, the first order kinematics is maintained. The dome formed in this simulation displays all the characteristics of a b-type dome.

The extensional step-over kinematics (Figure 5) leads to the localization of a dome in the center of the model. In this simulation, the deeper ductile part of the crust is exhumed with the shape of a dome beneath the topographic depression formed by normal and strike-slip faults in the upper brittle part of the crust. This dome is elongated in the direction of stretching imposed from the boundary condition. It is a non-cylindrical structure because the lineations (Figure 5a) are pointing out of plane when considering a cross section taken normal to the elongation of the dome. The lineations (red dots Figure 5b) within the dome are systematically horizontal and the schistosity planes (black dots Figure 5b) are approximately vertical, except at the very top of the dome forming an antiform structure. The particles show that this structure is slightly asymmetric with its steepest dipping limb located on the side of the strike-slip fault.

In this simulation, the strain is highly partitioned with depth. In the upper crust, the deformation takes place along steep grabens or pull apart basins (Figure 3b) which reflects the extensional nature of the boundary conditions. However, in the deeper part of the crust, the lineation forms a horizontal sigmoid compatible with pure left-lateral deformation, while the ductile exhumed lower crust within the dome deforms by constriction and forms an antiform (Figure 5a-b). The deformation in the exhumed part of the deeper crust is characteristic of an atype dome (Figure 1) with a strong local component of constriction, which was not imposed in the boundary conditions (Figure 3b). Dynamically, the formation of this structure is due to the local depression that forms above the extensional step-over. The large topographic gradient drives the flow in the ductile crust, attracting it towards the pull apart, where it rises to form the dome. As the dome grows, the strength of the upper crust drops locally, leading to strong strain localization.

The fault propagator kinematics (Figure 6) display the formation of two independent domes linked by a transfer zone. The first dome forms along the gray boundary with an almost cylindrical shape, similar to the dome of Figure 4a, although the lineations at the front of the structure are slightly oblique to the stretching direction (indicated by yellow colored cylinders instead of gray cylinders). The second dome forms at the location of the imposed strike-slip boundary condition. As for the a-type dome of Figure 5a, the structure is an asymmetric fold with a steeper dipping limb close to the fault. However, the lineations in the deeper part of the crust (Figure 6a) tend to strike oblique to the direction of stretching by an angle of 25 to 50° (red dots on Figure 6b) and indicate that this dome has been fed from both the back (cyan) and the right (black) side of the model.

This model with mixed boundary condition shares many common characteristics with model 1 and 2. However, in this mixed mode, the orientation of the stretching lineation in the first 10 km of the crust is much more spread out than for purely a-or b-type domes. In this mixed mode, the azimuth of the lineation may locally rotate by 40-50° as compared to the direction of extension or shear imposed from the boundary conditions.

P-T paths

The computation of synthetic P-T paths is common practice to validate models with observations in orogenic settings [START_REF] Gerya | Two-dimensional numerical modeling of pressure-temperature-time paths for the exhumation of some granulite facies terrains in the Precambrian[END_REF][START_REF] Yamato | Burial and exhumation in a subduction wedge: Mutual constraints from thermomechanical modeling and natural P-T-t data (Schistes Lustrés, western Alps)[END_REF]]. In both the case of the cylindrical and the fault propagator models (Figure 4c, 6c), despite the rather different shape of the domes in map view and the different kinematic imposed at the boundary, the P-T paths recovered from the two models are very similar. In both cases, the rocks are exhumed from all depths along the initial thermal gradient of 17°/km and starting from there, the exhumation takes place along isothermal or slightly heating P-T paths. In both models the final cooling of the exhumed rocks follows a linear gradient of c.a. 50°C/km which reflects the thinning of the crust.

P-T paths of the step-over model (Figure 4c) contrast with the first two. They do not display heating P-T paths, and the decompression arises with significant cooling. The thermal gradient after 12 Myr reaches only 33°/km as a result of the relatively small amount of thinning of the lower crust involved in these models as compared to model 1 and 3.

Summary and discussion of modeling results

Impact of initial conditions

An important feature of all these models is the rheological step that was introduced at mid crustal level. As discussed earlier, this step may represent a change in lithology due to nappe stacking as in Huet et al. [2011a;2011b], or could be caused by the presence of fluid [McKenzie and [START_REF] Jackson | Conditions for flow in the continental crust[END_REF] or be the limit of partially molten crust [START_REF] Lister | Plutonism and the origin of metamorphic core complexes[END_REF][START_REF] Vanderhaeghe | Partial melting and flow of orogens[END_REF][START_REF] Rey | Extension rates, crustal melting, and core complex dynamics[END_REF]. This step in viscosity is important for the outcome of the models as it allows localizing the deformation on a single detachment zone as shown and discussed in detail in Huet et al. [2011a] and it facilitates the rise of ductile material beneath the pull-apart in model 2. Similarly, we did not explore the effect of the initial depth of the step in viscosity, but we contend it would affect the initial spacing of the faults and therefore potentially change the outcome of the models by increasing the number of domes if the step was located at a shallower depth [START_REF] Bullard | Gravity measurements in East Africa[END_REF][START_REF] Meinesz | Les graben africains, résultat de compression ou de tension dans la croûte terrestre[END_REF][START_REF] Spadini | Spacing of consecutive normal faulting in the lithosphere: A dynamic model for rift axis jumping,Tyrrhenian Sea[END_REF]. If this step was to be at a lower depth, several domes would probably form initially, but at larger strain one would localize strain more efficiently [START_REF] Lavier | Factors controlling normal fault offset in an ideal brittle layer[END_REF][START_REF] Wijns | Mode of crustal extension determined by rheological layering[END_REF]. If several steps were to be included, we can predict that the localization of the deformation would probably jump from the deepest one to the shallower one over time as was observed in Huet et al. [2011a]

Rate of extension versus rate of exhumation

We observe that the model with transtensional splays (model 3, Figure 6c) exhumes lower crustal rocks to the surface along warmer geothermal gradients compared to pull-apart models (model 2, Figure 5c) for the same rate of extension. Part of this discrepancy is due to a geometrical effect of the design of the model because the pull-apart modelled here is infinite and therefore the source of material is limited to the left and right side, whilst in the case of a splay fault, significant amount of material is pumped from the back side (cyan side, Figure 6) of the model. Moreover, as the crust thins much more in the case of a splay fault (Figure 6) than in the case of a pull-apart, and considering that the strength of the lower crust or the mantle is not sufficient to maintain steps in the Moho, the lower crustal material located on the side with no extension is forced to flow towards the a-type structures. The second reason for the discrepancy is that as the dome gets wider, the diffusivity of the flow entering the a-type part of the dome increases as the cube of the thickness of the channel [START_REF] Buck | Modes of continental lithospheric extension[END_REF], enhancing exhumation rates.

One may also question the effect of changing the rate of extension on the outcome of our models. Based on analogue models of sand and silicone putty, [START_REF] Brun | Analogue modeling of detachment fault systems and core complexes[END_REF] and Brun

[1999] suggested that increasing the rate of extension decreases the rate of exhumation. This rate dependence is mainly due to viscous coupling which competes with the diapiric ascent of light silicone putty into the sand. As it was discussed in [START_REF] Wijns | Mode of crustal extension determined by rheological layering[END_REF], the change of the exhumation rate observed in analogue models is largely reduced when the buoyancy forces are small, which is the case in the model presented in this paper where buoyancy forces are kept to minimal and most of the exhumation is kinematically driven. Therefore, if one would decrease the extensional velocity, the rate of exhumation would probably decrease (but not drastically), in response to both kinematic forcing and reduced buoyancy forces in response to cooling by diffusion. The P-T paths would reach a cooler final thermal gradient, but as the rheological step remains, the overall kinematics would not change drastically. Similarly, increasing the rate of extension would probably increase the rate of exhumation and render the P-T path warmer in response to increased buoyancy forces, as was observed in models including melts [START_REF] Rey | Extension rates, crustal melting, and core complex dynamics[END_REF]. In any case, the structures and their orientations, which are central to our analysis, would remain unchanged. The presence of the rheological step thus render the overall model kinematics, or finite strain pattern, almost independent of the rate of extension.

Significance of a-type domes

We did not manage to create a-type domes with purely cylindrical boundary conditions, even though local components of constriction were observed in model 1. We note however that this could be an issue with the numerical grid resolution used in these models. It is possible that a-type structures could form at segments of the boundaries between domes if we considered similar instabilities as the one proposed by [START_REF] Gerya | Dynamical Instability Produces Transform Faults at Mid-Ocean Ridges[END_REF] for the formation of transform faults at mid oceanic ridges. We can also posit that similar dome shapes would occur if we had imposed the formation of two offset domes in the initial conditions such as those used in [START_REF] Choi | Thermomechanics of mid-ocean ridge segmentation[END_REF][START_REF] Allken | Three-dimensional numerical modeling of upper crustal extensional systems[END_REF]].

In any case, this does not contradict the main conclusion of our modelling study: when a transtensional strike-slip component is involved in the deformation of hot post-orogenic (i.e.

previously thickened) crust, a-type domes form. The most important result is that the constrictive ductile deformation within these domes is the mark of strike-slip transtensive faulting in the upper crust, rather than shortening normal to the stretching direction. We agree that at a first glance these structures resemble folds, however, the contact between the upper and the lower units is used as a detachment fault to exhume deeper units within the dome. Finally, the modelled topographic pattern indicates that sedimentary basin should form on the lateral side of a-type domes.

Comparison of the models to the Aegean MCCs

The Aegean Sea is a post-orogenic domain thinned in the back arc of the Hellenic subduction zone [START_REF] Le Pichon | The Aegean Sea[END_REF]]. In the Cyclades, the two types of Metamorphic Core Complexes have been exhumed with Naxos and Tinos being considered as representative of a-type and b-type domes, respectively [START_REF] Jolivet | Strain localization during crustal-scale boudinage to form extensional metamorphic domes in the Aegean Sea[END_REF] (Figure 7a). In this section, we discuss the applicability of our models to this domain and the implications for our understanding of the Aegean dynamics.

The P-T paths obtained are not sufficiently different from the natural P-T paths of Naxos [START_REF] Martin | Signification des ages U-Pb sur zircon dans l'histoire métamorphique de Naxos et Ikaria (Cyclades, Grèce)[END_REF][START_REF] Duchène | Pressure-temperature time evolution of metamorphic rocks from Naxos (Cyclades, Greece): Constraints from thermobarometry and Rb/Sr dating[END_REF] or Tinos [START_REF] Parra | Relation between the intensity of deformation and retrogression in blueschist metapelites of Tinos Island (Greece) evidenced by chlorite-mica local equilibria[END_REF] to reject either model based solely on thermo-barometric validation. However, in the fault propagator model, the behavior of the two parts of the dome differs. Whilst rocks are exhumed from all depths close to the free-slip boundary (dark gray, Figure 3f), within the strike-slip part of the dome, only high temperature paths are sampled (light gray, Figure 3f).

The lineations measured in the MCCs of the Cyclades are reported on Figure 7a. The outline of the modeled area is projected on the map so that the strikes of the lineation along the North Cycladic Detachment System (green arrows, [START_REF] Jolivet | The North Cycladic Detachment System[END_REF] are parallel to the direction of stretching of the models. In that frame, the roughly north-south lineations of Naxos, Paros and Ios (red arrows) [START_REF] Huet | Thrust or detachment? Exhumation processes in the Aegean: insight from a field study on Ios (Cyclades, Greece)[END_REF][START_REF] Gautier | Crustal-scale geometry and kinematics of late-orogenic extension in the central Aegean (Cyclades and Ewia Island)[END_REF]] form a 30-60° angle with the direction of stretching of the model and the late lineations along the South Cycladic Detachment System appear with intermediate strike (orange arrows, [Isgleder et al., 2009]). Not only is the geographic repartition of the lineation in good agreement with the fault propagator model (Figure 3c) but also the domes of Naxos, Paros and Ios fall actually in the class of domes that are elongated in the direction of the lineation. This observation, together with the warmer P-T paths in these domes and the migration to the south with time of from top-to-the-north-northeast detachments to a top-to-the-south detachment [START_REF] Grasemann | Miocene bivergent crustal extension in the Cyclades (Greece)[END_REF], leads us to conclude that the transtensional fault propagator model is a valid model for the formation of the Cyclades.

The obliquity of the lineation in the central Cyclades was previously attributed to a late phase of solid rotation around a vertical axis [START_REF] Morris | First palaeomagnetic results from the Cycladic Massif, Greece, and their implications for Miocene extension directions and tectonic models in the Aegean[END_REF]. The sinistral fault propagator model includes this obliquity into the continuum of the exhumation of the MCCs.

A strike-slip boundary of the Aegean-west Anatolian extensional system may have started ~25 Ma ago [START_REF] Jolivet | Aegean tectonics: strain localisation, slab tearing and trench retreat[END_REF]. Between 25 and 15 Ma, the total surface displacement of the strike-slip boundary can be measured by the amount of extension in the northern Menderes Massif resulting in some 50-60 km [START_REF] Van Hinsbergen | Reconciling the geological history of western Turkey with plate circuits and mantle tomography[END_REF]. Until now, a lot of putative structures have been proposed to accommodate strike-slip displacement in the area with different orientation and kinematics [START_REF] Gautier | Crustal-scale geometry and kinematics of late-orogenic extension in the central Aegean (Cyclades and Ewia Island)[END_REF], Walcott and White,1998, Philipon et al. 2012], however no large strike-slip structure has ever been clearly evidenced from field data or seismic reflection data. Our numerical experiments show that the Naxos structure is better explained as a transtensional structure rather than a purely extensional one. We argue that the shape of the dome of Naxos [START_REF] Vanderhaeghe | Structural development of the Naxos migmatite dome, in Gneiss Domes in Orogeny[END_REF], Kruckenberg et al., 2011] was acquired as it plays the role of a lower crustal root of a strike-slip fault which was active in the Miocene times.

Whether the a-type shape of the dome of Naxos marks a local segment jump between the et al., 2012] show that in the Miocene time, Anatolia and the Cyclades were probably further east than the tears current location. The strike-slip wrenching recorded in the a-type structure of Naxos could well be the lower crustal signature of a step fault accommodating the steep gradient in slab retreat rate caused by the slab tear (Fig. 8c). Within that hypothesis, the tear in the slab must initiate prior to 15 Myr, which corresponds to the late evolution of that dome [START_REF] Duchène | Pressure-temperature time evolution of metamorphic rocks from Naxos (Cyclades, Greece): Constraints from thermobarometry and Rb/Sr dating[END_REF].

General discussion and conclusions

Previous models for the formation of non-cylindrical domes in the crust involved components of normal shortening [START_REF] Thompson | Thermal evolution and exhumation in obliquely convergent (transpressive) orogens[END_REF], normal stretching [START_REF] Rey | Viscous collision in channel explains double domes in metamorphic core complexes[END_REF] or any complex poly-phased history which would include sequences of wrenching and folding or a late refolding of a pre-existing b-type dome with a shortening direction normal to the direction of extension [Philipon et al 2012]. All these equally valid models emphasize the non-uniqueness of the solution when trying to reconstruct the kinematics from finite-strain markers. Here, we show that the using a strike-slip fault propagator, or an extensional step-over in the boundary conditions of numerical experiments, together with initial conditions which favor the formation of MCCs, provides yet another valid alternative model for the emplacement of non-cylindrical domes. However, the main point we would like to emphasize is that the presence of constrictive ductile dome structures is not a sufficient argument to infer transpressive boundary conditions, or far field shortening in the absence of a geochronological argument for a more complex polyphased deformation history, as was proposed in earlier models.

21

Comparing the numerical experiments with the Aegean, we propose an alternative model to the formation of the Cyclades and relate the N-S axis of the Naxos dome to a phase of transtension. A possible source for the transtension could be the presence of a tear in the Aegean slab, but other hypothesis such as a transfer fault between the Cyclades and the Menderes massif (Western Turkey) remain as a possibility. In any case, we argue that the elongation and the structure of the Naxos dome is an indicator that the dome acted as a strike-slip structure during its activity. An argument to support this model is that on the opposite side of the Mediterranean, a-type domes are clearly observed in the Betics [START_REF] Augier | P-T-D-t retrograde evolution of the Nevado-Filabride complex (SE Spain) : new insights from in situ 40Ar/39 ages and metamorphic petrology[END_REF] and are again aligned to a tear propagating in a slab [Spackman and Wortel, 2004].

The numerical experiments and conclusions drawn in this paper are directly applicable to MCCs other than those in the Aegean. However, given the availability of data from the Aegean, this region represents an excellent location to validate our ideas and verify our 3D thermomechanical models. Our numerical experiments constitute a first order study of the mode of 3D deformation of thickened crust. Each of these experiments would require a systematic parametric study varying rates of extension or shear, relative timing between the onset of shear and extension, depth of the viscosity step, length scale of the step over in order to be applied and validated with data in different geological settings and geodynamic contexts.

The aim of this paper was not to perform such parametric studies, but to open the door to alternative models and interpretations for the 3D ductile deformation of post-orogenic, relatively thickened crust. With our simulations, we have shown that it is possible to form folds in response to local constriction within extensional step-over between strike-slip faults (model 2, Figure 5).

This model setup may well apply to the Red River Fault system along which many discontinuous exhumed domes are out-cropping [START_REF] Leloup | The Ailao Shan-Red River shear zone (Yunnan, China), Tertiary transform boundary of Indochina[END_REF]. Within this model, these domes would be the mark of extensional step-over within a purely strike-slip system. Similarly, many of the Variscan gneiss domes which form en echellon are usually interpreted as transpressional structures (e.g. Pyrenees [START_REF] Denèle | The Hospitalet gneiss dome (Pyrenees) revisited: lateral flow during Variscan transpression in the middle crust[END_REF]). Our model shows they could alternatively have formed below an extensional step over within a purely strike-slip, large-scale kinematic boundary conditions. The formation of these domes during the early stage of wrenching could also participate in weakening the strike-slip faults by advecting hot/weak material within the fault zone and locally reducing the effective elastic thickness of the fault zone as was proposed by [Chery, 2008] to explain the GPS signal around the San Andreas Fault.

Further work is needed to completely understand the dynamics of strike-slip systems in the ductile part of the crust, however as 3D thermo-mechanical modeling is now becoming readily available, we can expect much progress to be made in the near future.

APPENDIX

A Computation of synthetic tectoglyphs

The method used for the computation of the synthetic finite-strain field is an extension to 3D of a classical 2D method of Ramsay and Huber [1987]. The synthetic finite-strain field is computed over a set of passive markers, which is meshed with tetrahedra. The schistosity and lineation are computed for all tetrahedra in the following way.

Each tetrahedron is described by three non-colinear vectors that connect the four vertices.

These vectors are noted

a i = x i , y i ,z i T , i = 1,2,3 (A.1)
in the initial (undeformed) configuration and

K G D C u p = f h (B.1)
The velocity pressure (u,p) solution is obtained via Schur Complement Reduction. This entails applying a Krylov method to the Schur complement system

DK 1 G C p = ˆ f , (B.2)
Where

ˆ f = DK 1 f h , (B.3)
to obtain the solution for p. Following this, we apply another Krylov method to

K u = f Gp , (B.4)
in order to obtain the solution for u. In GALE, the Schur complement system in (B.2) is solved using the Conjugate Gradient method. Iterations are terminated when the initial residual has been reduced by a factor of 1e5. Systems involving K are solved using FGMRES, using a stopping condition requiring that the initial residual be reduced by a factor of 1e6. The same Krylov method used to define the matrix-vector product in Eqn. The performance of a Krylov method is strongly linked to the preconditioner used. The ideal preconditioner is both scalable and optimal; meaning that number of iterations required to reach convergence is constant as the finite element mesh is refined, and that the CPU time required to reach convergence scales linearly with the number of unknowns. For the system in Eqn. (B.2), GALE uses a preconditioner defined by the mass matrix, scaled by the inverse of the element-wise effective viscosity. For Q1-Q1 velocity-pressure spaces, this preconditioner has been proven [START_REF] Grinevich | An Iterative Method for the Stokes Type Problem with Variable Viscosity[END_REF] and demonstrated for geodynamic applications [START_REF] Geenen | Scalable robust solvers for unstructured FE geodynamic modeling applications: Solving the Stokes equation for models with large localized viscosity contrasts[END_REF][START_REF] Burstedde | Parallel scalable adjoint-based adaptive solution of variable-viscosity Stokes flow problems[END_REF] to produce iteration counts for the Schur complement system that are independent of the mesh resolution.

To achieve an overall optimal solution strategy for the complete Stokes solver, GALE employs a geometric multi-grid (GMG) preconditioner for systems of the form K y = x. A geometric multi-grid preconditioner utilizes a hierarchical representation of the discrete operator K. The simplest hierarchy contains two levels, consisting of K (fine level) and K c , a "coarse" grid operator which is defined on a coarser finite element mesh. In the classical multi-grid strategy, one removes high frequency components of the error in the solution by applying several iterations of an iterative method, such as Richardson+Jacobi or Richardson+Gauss-Seidel. Due to their properties, these methods are called "smoothers". Since only long wavelengths remain, the smoothed error can be represented on a coarser grid. The inter-grid transfer from the fine-tocoarse level is called "restriction". On the coarse grid, the exact error can be computed using a direct solver, as the coarse problem contains fewer degrees of freedom than the fine grid. The low-frequency errors on the coarse grid are transferred to the fine grid in order to correct the fine grid approximation for the low-frequency errors. The inter-grid transfer from the coarse-to-fine grid is called "prolongation". The smoother is also applied to the corrected fine grid approximation. We refer to [START_REF] Wesseling | An Introduction to Multigrid Methods[END_REF], [START_REF] Briggs | A multigrid tutorial[END_REF] and [START_REF] Trottenberg | Multigrid[END_REF] for an in-depth discussion of multi-grid theory.

To obtain "text-book" efficiency multi-grid performance, the coarse grid operator has to be a meaningful approximation of the fine grid operator K. For partial differential equations (PDEs) containing coefficients which are highly spatially variable (continuous or discontinuous)

and heterogeneous (e.g. viscosity), this is difficult within the framework of GMG. Primarily this is due to the inability of the coarse grid to adequately resolve fine scale structures. Additionally, in practice, the classical smoothers are found to be largely ineffective for problems which possess large variations in coefficients. We note that in the instance of discontinuous coefficients, if the coarse grid exactly resolves the coefficient jump, all the aforementioned issues are eliminated. However, for most practical geodynamic models employing markers -this scenario will rarely ever occur. In summary, the construction of the coarse grid operator and the choice of the smoother are extremely important for the development of efficient GMG preconditioners when they are applied to problems which exhibit large variations in coefficients.

In GALE, the restriction operation R is defined via bilinear interpolation between the fine and coarse grid levels. The prolongation operator P is defined as the transpose of R. Galerkin coarse grids operators are utilized on all levels (except the finest) in which the coarse grid operator is constructed via K c = RK R T . This type of coarse grid can be interpreted as a projection of the entire fine grid matrix onto the coarse grid. In general, practitioners often don't utilize Galerkin coarse grid operators due to the programming complexity (particularly in parallel) of defining the triple matrix product required. One alternative to Galerkin coarse grid operators is simply to project the coefficients onto the coarse grid and re-discretize the original PDE. However, in practical applications with large variations in coefficients, the Galerkin coarse grid operator produces a much more robust preconditioner. In addition, the "smoother" used in GALE is a fully-fledged Krylov method (FGMRES) equipped with a block Jacobi+ILU preconditioner. On both the down and upward sweep in the multi-grid V-cycle, we apply 8 iterations of this Krylov method. In contrast to the classical Richardson+(Jacoboi/Gauss Seidel/SOR) type smoothers, our practical experience of using such "heavy" smoothers have proved to be robust when applied to problems involving large variations in viscosity. Furthermore, on the coarsest grid level, we utilize a fully parallel LU factorization SuperLU_Dist [Xioye, 2005]. All linear algebra, Krylov methods and preconditioners are provided via PETSc v3.0 [START_REF] Balay | PETSc home page[END_REF].

FIGURE CAPTIONS

Figure 1: On the right, a 3D sketch (modified from [START_REF] Jolivet | Strain localization during crustal-scale boudinage to form extensional metamorphic domes in the Aegean Sea[END_REF]) illustrates two classes of domes. On the left, a sketch illustrates the stereo-plot projections of the lineation (L) and foliation (S) for the two kinds of dome. In the a-type dome, constriction is important and the foliation is folded with axis aligned with the direction of stretching. In the b-type dome, the foliation is folded with an axis normal to the direction of extension. The a-type domes are not cylindrical and cannot be modeled in 2D. stands for upper crust, lower crust and mantle respectively; c. Yield strength envelopes drawn in extension for a background strain rate equal to the mean strain rate (black line, 1.5x10 -15 s -1 ) and ten times the mean strain rate (red line, 1.5x10 -14 s -1 ). [START_REF] Jolivet | Strain localization during crustal-scale boudinage to form extensional metamorphic domes in the Aegean Sea[END_REF] ) illustrates two classes of domes. On the left, a sketch illustrates the stereo-plot projections of the lineation (L) and foliation (S) for the two kinds of dome. In the a-type dome, constriction is important and the foliation is folded with axis aligned with the direction of stretching. In the b-type dome, the foliation is folded with an axis normal to the direction of extension. The a-type dome are not cylindrical and can not be modeled in 2D.

Figure 2 : a. Initial geotherm is linear with a thermal gradient of 17°C/km down to the 1300°C isotherm and constant deeper; b. lithostratigraphic column with density u.c., l.c. and m. stands for upper crust, lower crust and mantle respectively c. Yield Strength Envelopes drawn in extension for a background strain rate equal to the mean strain rate (black line, 1.5x10-15 s-1) and ten times the mean strain rate (red line, 1.5x10-14 s-1). [START_REF] Parra | Relation between the intensity of deformation and retrogression in blueschist metapelites of Tinos Island (Greece) evidenced by chlorite-mica local equilibria[END_REF] and Naxos [START_REF] Duchène | Pressure-temperature time evolution of metamorphic rocks from Naxos (Cyclades, Greece): Constraints from thermobarometry and Rb/Sr dating[END_REF]

  Menderes massif and the Northern Cycladic Detachment System, or whether it is driven by global geodynamics is a matter of debate which we cannot quantitatively assess based on the results of our numerical experiments. However, to conclude on the Aegean dynamics, we would like to propose a completely qualitative and conceptual, yet deliberately provocative model for the formation of the splay fault that we imposed in the boundary condition of our preferred Aegean model. A tear within the Aegean slab is documented by Berk-Biryol et al. [2011] under western Turkey and propagating along the Pliny and Strabo trenches south east of Crete and corresponds today to a steep lateral gradient in trench retreat rate. The timing of the initiation of this tear is debated. Pe-Piper et Piper [2007] and Jolivet et al. [2009] evaluate it at c.a. 20 Myr based on adakites ages (Fig. 8 a. and b.) while van Hinsbergen et al. [2010] posit this tear is not older than 15 Myr based on palinspatic plate reconstruction. Reconstruction of the Aegean area [Philippon

  (B.2) is used in both Eqn. (B.3) and Eqn. (B.4).

Figure 2

 2 Figure 2 : a. Initial geotherm is linear with a thermal gradient of 17°C/km down to the 1300°C

Figure 3 :

 3 Figure 3: Kinematic boundary conditions, surface deformation and topography with equidistance

Figure 4 :

 4 Figure 4: Results for cylindrical extension (Model 1). a. Internal deformation of the model

Figure 5 :Figure 6 :

 56 Figure5: Results for extensional step over (Model 2), legend is the same as for Figure4

Figure 8 :Figure 1 :

 81 Figure 8: a and b Paleogeographic reconstructions of the Cyclades and slab tear simplified from

Figure 3 :Figure 4 :

 34 Figure 3: Kinematic boundary conditions, surface deformation and topography with equidistance of 100m for the 3 models after 12 Myr evolution. Brown lines initially formed squares. a. Model 1 assumes cylindrical extension b. Model 2 corresponds to a 60 km wide extensional step over c. Model 3 represents a trans-tensional fault propagator.

Figure 5 :

 5 Figure 5: Results for extensional step over, Model 2, legend is the same as for Figure 4

Figure 6 :Figure 7

 67 Figure6: Results for trans-tensional fault propagator, Model 3, legend is the same as for Figure4

Figure 8

 8 Figure7: a. Topographic map of the Cyclades with outline of the Cycladic MCCs (purple area), as well as stretching lineations and sense of shear (arrows drawn after[START_REF] Martin | Signification des ages U-Pb sur zircon dans l'histoire métamorphique de Naxos et Ikaria (Cyclades, Grèce)[END_REF],[START_REF] Huet | Thrust or detachment? Exhumation processes in the Aegean: insight from a field study on Ios (Cyclades, Greece)[END_REF],[START_REF] Iglseder | I and S-type plutonism on Serifos (W-Cyclades, Greece)[END_REF],[START_REF] Jolivet | The North Cycladic Detachment System[END_REF] andRing et al. [2011]). The di erent colors outline di erent provinces associated with typical strikes with a similar color code as in Fig 4.The modeled area is reported with white dashed lines, black boxes denotes where the cylindrical part and non cylindrical part where sampled in the model 3 to produce the part c and d of the gure; b: geological sketch representing the main structural features of model 3; c: the blue and magenta lines are the natural P-T paths of Tinos[START_REF] Parra | Relation between the intensity of deformation and retrogression in blueschist metapelites of Tinos Island (Greece) evidenced by chlorite-mica local equilibria[END_REF] and Naxos[START_REF] Duchène | Pressure-temperature time evolution of metamorphic rocks from Naxos (Cyclades, Greece): Constraints from thermobarometry and Rb/Sr dating[END_REF] respectively. Light and dark gray lines denote paths located within non-cylindrical and cylindrical part of model 3 (see a); d: stereo-plot of the stretching lineation (red) and pole of the foliation (black) in the cylindrical (left) and and non-cylindrical (right) part of model 3 oriented in the same frame as the map in a. The cylindrical part is a typical b-type dome while the non cylindrical part resemble a-type with non negligible number of lineation oriented N-S like in Naxos.

TABLES 761 Table 1 :

 7611 Model parameters

	762
	Notation/Unit	Lower crust Upper crust Mantle
	0 (Pa.s)	10 26	10 26	10 26
	(°C -1 )	0.035	0.020	0.012
	eff at 300°C (Pa.s) 3x10 21	2x10 23	3x10
	eff at 700°C (Pa.s) 10 19	8x10 19	2x10
	eff at 1100°C	10 19	10 19	2x10
	(Pa.s)
	Co (Pa)	2x10 7	2x10 7	2x10 3
	Co inf (Pa)	2x10 7	2x10 7	2x10 3
	(°)	30	30	15
	inf (°)	10	10	15
	(K -1 )	3x10 -5	3x10 -5	3x10 -5
	(m 2 s -1 )	10 -6	10 -6	10 -6
	(kg m -3 )	2800	2800	3300
	763

Table 2 : Scaling used for computations

 2

	764
	Notation/Unit Model Nature (SI) Scaling (SI)
	Velocity [LT -1 ] 10 5	10 -10	10 -15
	Distance [L]

ACKNOWLEDGMENTS

GALE is software hosted by the Computational Infrastructure for Geodynamics (CIG) and developed by CIG, Monash University and the Victorian Partnership for Advanced Computing (VPAC). This research was supported by ANR EGEO. Author D.A.M. was supported by the ETH Zurich Postdoctoral Fellowship Program. The topographic map was generated using http://www.geomapapp.org

in the final (deformed) configuration, where the sign T designates the transpose operator. The two sets of vectors are gathered in two square matrices,

for the initial configuration and

for the final configuration.

The deformation matrix D, which relates the initial coordinates of the tetrahedron vectors to the final ones, links the matrices u and U:

The finite strain ellipsoid is characterized by the eigenvalues I 1 I 2 I 3 and the corresponding eigenvectors v 1 , v 2 and v 3 of matrix

The major (resp. intermediate and minor) principal axis of the finite strain ellipsoid is parallel to v 1 (resp. v 2 and v 3 ). The length of the major (resp. intermediate and minor) principal axis is I 1 1 2

(resp. I 2 1 2 and I 3 1 2 ). Finally, the stretching lineation is a line parallel to v 1 and the schistosity is a plane normal to v 3 . All these computations can be done using the script tectoglyphs.m provided with the manuscript.

B Multi-grid in GALE

Within GALE, the discretized Stokes equations are represented as follows;


Emmanuel Duguet
email: emmanuel.duguet@u-pec.fr

Loïc

Pascale Petit

Keywords: JEL classification: C93, J16, J61, J71 gender, origin, hiring discrimination, first order stochastic dominance

Introduction

Correspondence testing is an appropriate technique for evaluating hiring discrimination in the labor market. Economists have been using this methodology since the end of the seventies [START_REF] Firth | Racial discrimination in the British labor market[END_REF][START_REF] Firth | Sex Discrimination in Job Opportunities for Women[END_REF]. The technique is a controlled experiment that consists of fabricating written applications (CVs and cover letters) of fictitious candidates. The applications are identical, except for one characteristic a priori not tied to productivity, like the applicant's gender or origin. The applications are sent to the same job offers, and one examines whether the applicants obtained a similar access to hiring interviews. [START_REF] Riach | Field Experiments of Discrimination in the Market Place[END_REF] highlighted two trends emerging from the results of the correspondence testings conducted in order to evaluate the degree of gender discrimination in hiring. First, women suffer from discrimination in hiring for well-paid positions involving responsibilities. In Great Britain, for instance, [START_REF] Firth | Sex Discrimination in Job Opportunities for Women[END_REF] found that women were less likely than men to obtain a position as a qualified accountant; in Philadelphia, NEUMARK et al. (1996) found significant discrimination against women for server jobs in high-price restaurants and significant discrimination against men for server jobs in low-price restaurants. Second, discrimination can be observed in jobs in which one gender is over-represented. For example, there is discrimination against men in secretary jobs traditionally occupied by women; conversely, there is discrimination in hiring against women in activities traditionally occupied by men, like mechanics (WEICHSELBAUMER, 2004). Lastly, [START_REF] Bertrand | Are Emily and Greg More Employable Than Lakisha and Jamal? A Field Experiment on Labor Market Discrimination[END_REF] stressed that the scale of discrimination in hiring varies according to the ethnic origin of the candidates. Their testing evaluated the scale of racial discrimination according to the gender of the job applicants. The authors found that candidates with Caucasian-sounding names receive 50% more callbacks than their counterpart with African-American sounding names.

From a methodological viewpoint, the correspondence testing literature consists in sending resumes to a recruiter that invites the candidates to an interview or not. Consider the case of two candidates, A and B. Four answering cases are possible: no candidate is invited, only candidate A is invited, only candidate B is invited or both candidates are invited. The standard discrimination measure is a calculation of the following discrimination coefficient. One computes the difference between the percentage of offers for which candidate A had been invited without candidate B and the percentage of offers for which candidate B had been invited without candidate A. With this method, the cases were both candidates had been invited are not used in the measurement of discrimination because both candidates are considered to have been equally treated.

Our aim is to provide a more general measure of hiring discrimination. We claim that the cases in which both candidates have been invited should also be included in the analysis when information is available about the order in which the candidates have been called. The justification is that in the standard testing procedures, the fictitious candidates have instructions to respond that they have already found a job when they are called, so that they always decline the job offer. Therefore, the recruiters decide to go further down the short list of candidates, and thus reveal a ranking of the candidates. In this paper we propose using the condition of first order stochastic dominance (henceforth, FOSD) on the candidates' ranking. Consider the case with ≥ 3 candidates. Candidate A FOSD candidate B when he/she has higher probabilities to rank first than B, to rank among the first two candidates than B, to rank among the first three candidates than B etc. We interpret the condition of FOSD as the existence of strong discrimination against candidate B in this paper, because we show that FOSD involves a stronger form of discrimination that the one measured by the standard discrimination coefficient.

We apply this new method to a correspondence testing conducted in 2009 which aimed to measure both gender and ethnic origin discrimination among French young workers in the Paris region. We fabricated eight similar resumés: four origins for each gender (French, Moroccan, Senegalese, and Vietnamese). The origin is indicated by the sound of the first and last names. The eight candidates all hold French nationality, live in Paris, and apply for the same jobs as software developers at the Master's degree level. Applications have been sent simultaneously to the same 310 job offers. This protocol enables us to examine first, whether hiring discrimination is based on origin, and second, whether discrimination against women varies with their origin.

We show that, while the standard discrimination coefficients point in the right direction in most of the cases, they tend to underestimate the importance of discrimination in the labor market. The next section presents the methodology, and the following presents the application and summarizes our results.

Methodology

Standard discrimination coefficient

Consider a recruiter with preferences for the candidates A and B represented by the utilities and . These utilities are specific to each recruiter and result from pre-conceptions about the candidates, because the candidates are equally productive by construction of the testing experiment.

Each recruiter has a reservation utility level above which the candidates are invited to an interview. We define the relative utility levels = -and = -. The four potential answering cases can be represented in the following way. If < 0 and < 0, no candidate is invited to an interview. When < 0 < , only candidate B is invited; when < 0 < , only candidate A is invited. Finally, when > 0 and > 0, both candidates are. These cases are illustrated in Figure 1. The standard measure of discrimination against candidate B, used in the literature, considers only cases in which only one of the two candidates is invited. These cases are illustrated by the North-West and the South-East quadrants of Figure 1. We denote this discrimination coefficient as ∆ A, B :

∆ A, B = P < 0 < -P < 0 < = P A invited, B uninvited -P B invited, A uninvited

According to this measure, there is no discrimination when both candidates have equal chances to be invited, and a positive number indicates that candidate A is -on average -preferred to candidate B. 1

Also notice the property that ∆ A, B = -∆ B, A .

First order stochastic dominance and strong discrimination

We propose to extend the standard measure of discrimination ∆ to the ranking of the candidates when both are invited, which is equivalent to consider all the quadrants in Figure 1. 1 Since ∆ measures discrimination when only one of the two candidates is invited to the interview, we use the subscript 1.

45°

Both invited: Man preferred

> ! > 0 Both invited: Woman preferred ! > > 0 Woman invited < 0 < !

None invited

! and < 0 ! < 0 <

Man invited

!

An important point to mention is that we need to consider the ranking among ≥ 3 candidates, since our application includes 8 candidates. Therefore we develop our analysis in the general case of an undetermined number of candidates k. Among these k candidates, we perform all the pairwise comparisons available (here :

-1 /2 = 28).

In order to compare the rankings of two candidates, we use the concept of first order stochastic dominance. Suppose that the k candidates are ranked according to the recruiter's utilities. The candidates that have not been invited satisfy the condition that ' = ' -< 0. The ranking of the candidates (from 1 st to k th ) results from the recruiter's tastes for the candidates. The highest utility corresponds to the candidate ranked first, and negative utilities correspond to the candidates that have not been invited. In order to perform our analysis, we need to separate the candidates that have not been invited from the others by creating a rank k+1. This additional rank is required because the candidates that have not been called cannot be ranked between themselves. This also allows us to defined ∆ directly. We know only that the uninvited candidates' utilities are below the recruiter's reservation utility levels and therefore that they are ranked behind the candidates that have been invited.

Consider first the case for which all the candidates have been invited. Using the order statistic2 , we obtain the ranking of the utilities of the candidates, 0 < ≤ ) ≤ ⋯ ≤ + , that corresponds to the ranking , -1, … ,1. When only j candidates are invited, we have the ranking ≤ … ≤ +-' < 0 < +-'. … ≤ + that corresponds to the ranking + 1, … , + 1, 0, … ,1.

The first stochastic dominance of candidate A over candidate B is defined as: This relationship is especially easy to interpret when is set at the reservation utility level of the recruiter, since it means that candidate A has a higher probability to be invited to the interview than candidate B. We also see that FOSD covers more cases than the standard discrimination measure because it makes use of all possible utility thresholds.

Pr ≥ ≥ Pr ≥ ∀ ,
For practical reasons we will work with the ranks, since they are observable while the utilities are not. We just need to reverse the inequalities inside the probabilities, since the higher the utility the lower the rank (rank 1 for the most preferred candidate with utility + :

Pr : ≤ : ≥ Pr : ≤ : ∀ : ∈ <1, … , + 1= and ∃ :̅ such that Pr : ≤ :̅ > Pr : ≤ :̅ Consider the case : = 1. Then Pr : ≤ 1 = Pr : = 1 , which gives the probability to be ranked first. If the corresponding inequality holds, the candidate A has a higher probability to be ranked first than the candidate B. Now set : = 2. We conclude that candidate A has a higher probability to be ranked among the two first candidates than candidate B. Performing the comparisons up to : = , Pr : ≤ is the probability that candidate A is invited to the interview. Therefore candidate A has a higher probability to be invited to the interview than candidate B. In summary, when A FOSD B, the candidate A always has a higher probability to be in the leading group than the candidate B, whatever the definition of the leading group is. This definition is especially relevant for the (b)The converse of this implication is false.

In order to simplify the exposition, we introduce the following notation, considering candidates A and B among k candidates:

Pr : = @ ∩ : = B = Pr @, B , and we use the two following conventions:

1) Each of the k candidates has k+1 possible ranks: 1 st ,2 nd ,…,k th or uninvited which we denote (k+1) th .

2) When there are no ties, we set Pr @, @ = 0, which means that two candidates cannot have the same rank. The existence of ties does not affect the property. 3 The CDF defined over the ranks is equivalent to the survival function defined on the utilities.

On the other hand, it is always true that Pr : C ≤ : = ∑ Pr :

C = E F G , H ∈ <I, J=,
and

Pr : C = E = ∑ Pr E, 0 +. 'G , so that Pr : C ≤ : = ∑ ∑ Pr E, 0 +. 'G F G
.

Proof:

(a) Notice that Pr : C ≤ is the probability that candidate i is invited, regardless of the situation of the other candidates. FOSD implies that:

Pr : ≤ ≥ Pr : ≤ ⇔ ∑ ∑ Pr E, 0 +. 'G + G ≥ ∑ ∑ Pr 0, E +. 'G + G ⇔ ∑ ∑ Pr E, 0 + 'G + G + ∑ Pr E, + 1 + G ≥ ∑ ∑ Pr 0, E + 'G + G + ∑ Pr + 1, E + G
(1), and the sum ∑ ∑ Pr E, 0

+ 'G + G
adds the probabilities of all the possible couples E, 0 between 1 and k. Therefore :

∑ ∑ Pr E, 0 + 'G + G = ∑ ∑ Pr 0, E + 'G + G
and inequality (1) simplifies to:

∑ Pr E, + 1 ≥ ∑ Pr + 1, E + G + G ⇔ Pr A invited, B uninvited ≥ Pr B invited, A uninvited ⇔ ∆ A, B ≥ 0 (b)
The second part of the proof is easier. We have shown in part (a) that Δ A, B ≥ 0 is a consequence of the last inequality that defines FOSD, and therefore it does not imply FOSD.

Q.E.D.

This allows us to define strong discrimination against candidate B as A FOSD B, and weak discrimination as ∆ A, B ≥ 0 or equivalently ∆ B, A ≤ 0. Weak discrimination against candidate B simply states that the probability to be invited to an interview alone is higher for candidate A than for candidate B, while strong discrimination implies that the probability that candidate A has a good ranking is higher than for candidate B. Therefore weak discrimination uses only the extreme case where one candidate only is invited, while strong discrimination also includes all the cases where both candidates are invited. This has three interesting consequences.

First, two candidates with similar values for ∆ might or might not display FOSD. This opens the possibility to classify candidates with similar discrimination coefficients according to the FOSD criterion. Second, the importance of ∆ is not always a measure of the importance of discrimination, because a candidate experiencing strong discrimination can have a lower discrimination coefficient than one experiencing weak discrimination. We find evidence for both properties in our application.

Third, it is possible to provide an overall ranking of the discrimination faced by the candidates by using the conditions of weak and strong discrimination together. We will first rank the candidates according to the strong discrimination criterion and, when they do not FOSD each other, we will rank them according to the weak discrimination criterion. Combining these two measures, we are able to indicate each candidate's position on an ordinal discrimination scale.

Figures 2 and3 provide two examples of the presence or absence of FOSD taken from our data. In

Figure 2, we compare the Vietnamese origin candidates, and find that women dominate men. But it is also possible that two candidates do not dominate each other; this happens when one CDF crosses the other. In Figure 3, we compare the French origin candidates and find that none dominates the other. The CDFs of all the candidates are presented in Table 4, and the FOSD analysis in Tables 5A to 5C (28 pairwise comparisons).

Application and results

We performed a correspondence testing in the Paris region. It has been designed to assess the effects of gender and ethnic origin on access to hiring interviews. Eight CVs of fictitious candidates with similar characteristics regarding productivity were constructed for young software developers with Masters' degrees. These candidates, all of French nationality, differed solely by gender and ethnic origin indicated by their names (French, Moroccan, Senegalese or Vietnamese).

Data collection

This stage was designed with a few challenges in mind. First, we wanted to limit the probability of detection when sending eight CVs simultaneously. Second, we also needed to consider a job with a high labor demand in order to have high response rates. This methodological precaution proved to be particularly useful in a context of economic recession. For this purpose, we consulted the Historic Data File of the French unemployment agency (Pôle Emploi) in order to find a profession where both the number of candidates and the labor demand was high. The profession that met these criteria is software development.

The 8 fictitious candidates explicitly state their French nationality on their CVs; their first and last names indicate their gender and ethnic origin. The applications were sent in response to the same job offers, and were equivalent in terms of productivity-related characteristics. The candidates possessed identical diplomas, professional careers and experience, and the same computing and language skills. None of the candidates had been unemployed, and they all had jobs at the time of their applications. In addition, we had these fictitious CVs checked by recognized professional experts in the field. More details on the experimental design are given in Appendix 1.

Because these applications were sent at the same time in response to the same job offers, they had to include some elements of differentiation. These differences involved the presentation of the CV:

font type and size, page layout, while still remaining standard. The candidates declared professional experience acquired in real companies, which were different but comparable (in terms of activity, size and market power). The candidates hobbies were also different, while remaining very standard and impersonal (sport, cinema, reading, music, etc.). The short cover letters accompanying the CVs were also worded differently, while remaining standard. Each candidate was attributed a postal address, a mobile telephone number, and an e-mail address.

To ensure that the style or content of a particular application did not systematically influence companies' choices with respect to a particular candidate -in spite of the precautions taken in constructing them -we used a random permutation of the CVs between the identities of the fictitious candidates. The CVs and cover letters were thus alternated between the different candidates.

Applications for each job offer were sent on the day the offer was posted on the Internet, at an interval of several minutes, by e-mail from each candidate's mailbox. The response is considered positive when the recruiter invites a candidate to an interview or asks for further details about the candidate's current situation or qualifications. On the contrary, the response is considered negative if the recruiter ignores or formally rejects the application. When a candidate is invited to an interview, he/she declares that he/she already found a job and declines the offer.

Results

Since the candidates share a high level of qualification and a similar experience in a growing line of business, we did not expect to find a significant degree of hiring discrimination. Even if there were some discriminatory practices, they should be restrained by the difficulty in finding workers. That, however, is not at all what we find; discrimination clearly drives the data.

First, discrimination by origin is the most important type both by the range of candidates it applies to and by its degree. Second, gender discrimination is present for African origin women but is clearly weak or absent for women of French and Vietnamese origin.

We present the standard discrimination coefficients first and show how accounting for strong discrimination (FOSD) puts some light on these first-pass results. This will clearly show the contribution of this new measure to the analysis of hiring discrimination.

Table 2 presents the gender standard discrimination coefficients separately for each origin. We compare the success rate of each woman to the success rate of the man with the same origin. A positive coefficient indicates discrimination against women, while a negative coefficient indicates discrimination against men. We find that women with a French or African (Senegal, Morocco) origin are discriminated against, while men with a Vietnamese origin are discriminated against women of the same origin. The magnitudes of the discrimination coefficients are not comparable, however, since they use a different reference (the man of same origin varies from one coefficient to another).

For instance French-origin women have a discrimination coefficient of 5,2%, and Senegalese-origin women have a coefficient of 5,5%, but they do not face the same situation at all. If we use the strong discrimination criterion, we find that French-origin women are not FOSD by French-origin men (Table 5A), while Senegalese-origin women are FOSD by Senegalese-origin men (Table 5C). This is one case in which the new measure is useful. The two remaining cases (Moroccan and Vietnamese origins) fit with the FOSD criterion, so that Moroccan-origin women (Table 5B) and Vietnamese-origin men face strong discrimination (Table 5C). We can therefore propose the following summary: French-origin women face weak gender discrimination, while African-origin women and Vietnamese-origin men face strong gender discrimination. We will show later that it is possible to provide a full ranking of the candidates according to the weak and strong discrimination criteria.

For now, let A ≽ N B denote a strong preference for candidate A (a strong discrimination against candidate B), and A ≽ O B denote a weak preference for candidate A (a weak discrimination against candidate B). We have found that:

French-origin men ≽ O French-origin women Moroccan-origin men ≽ N Moroccan-origin women Senegalese-origin men ≽ N Senegalese-origin women Vietnamese-origin women ≽ N Vietnamese-origin men One interesting property of the strong and weak discrimination criteria is that they allow ranking all the candidates of the testing along a single discrimination ordinal scale.

Table 3 presents the origin discrimination coefficients separately for each gender. Whether we consider women or men, one result never changes: foreign-origin candidates are always discriminated against to the benefit of the French-origin candidates. But to obtain this result, we need the FOSD criterion. Indeed, a close look at the discrimination coefficients reveals that the coefficient of Vietnamese-origin women is positive (2,6%) but not significant, so that, according to this coefficient, there would be not discrimination against Vietnamese-origin women. But a look at Table 5A shows that French-origin women FOSD Vietnamese-origin women, so that there is strong discrimination against Vietnamese-origin women. If we compare all of the candidates according to the strong discrimination criterion, we obtain a partial ordering:

(French-origin men, French-origin women) ≽ N (Moroccan-origin men, Vietnamese-origin women)

≽ N (Senegalese-origin men, Vietnamese-origin men)

≽ N Moroccan-origin women ≽ N Senegalese-origin women This partial ordering provided by FOSD can be completed by applying the weak discrimination criterion. We have already seen that French-origin men are preferred to French-origin women. To complete the ordering, we just need to compare Moroccan-origin men to Vietnamese-origin women, and Senegalese-origin men to Vietnamese-origin men. This can be done by using the additive property of weak discrimination: ∆ A, C = ∆ A, B + ∆ B, C . Using this property and the information contained in Tables 2 and3, we obtain the complete ranking of the candidates by the recruiters:

French-origin men ≽ O French-origin women A potential critique that could apply to our methodology is the possible use of alphabetical ordering by recruiters as opposed to a discriminatory ordering. Therefore, we have checked whether alphabetical ordering had been used for each job offer and found that the data patterns are compatible with alphabetical ordering in only 11.5% of the cases (Table A1 in appendix 2). Therefore, the results that we present may reflect the recruiters' ranking rather than the neutral alphabetical ordering.

Conclusion

We have performed a correspondence testing in a line of business with a high job vacancy to applicant ratio. Considering young workers with a Master's degree and no period of unemployment, we did not expect to find evidence of significant discriminatory practices. The experiment indicates the contrary: overall, we do find evidence of significant discrimination by gender and origin.

Moreover, in most cases, the degree of discrimination found is strong, since it corresponds to a pattern of stochastic dominance. While the standard discrimination coefficients go in the expected direction in most of the cases, they tend to underestimate the importance of discrimination in the French labor market. French-origin candidates of both genders clearly have more chances to be invited to a hiring interview than foreign-origin candidates, and African-origin women have the lowest chances to be invited.

Introducing strong discrimination (FOSD) brings value added to the analysis for two reasons. First, we find that candidates with similar standard discrimination coefficients do not face the same situation in the labor market, since some are weakly discriminated against and others are strongly To be continued 1 -0,045 0,003 -0,013 -0,048 -0,016 0,048 0,032 -0,003 0,029 2 -0,048 0,010 -0,048 -0,061 -0,016 0,058 0,000 -0,013 0,032 3 -0,077 0,019 -0,045 -0,077 -0,010 0,097 0,032 0,000 0,068 4 -0,074 0,026 -0,032 -0,100 -0,013 0,100 0,042 -0,026 0,061 5 -0,087 0,035 -0,035 -0,097 -0,003 0,123 0,052 -0,010 0,084 6 -0,094 0,032 -0,032 -0,094 -0,006 0,126 0,061 0,000 0,087 7 -0,106 0,026 -0,035 -0,094 -0,013 0,132 0,071 0,013 0,094 8 -0,106 0,019 -0,035 -0,097 -0,019 0,126 0,071 0,010 0,087 To be continued 1 -0,016 -0,052 -0,019 -0,035 -0,003 0,032 2 -0,058 -0,071 -0,026 -0,013 0,032 0,045 3 -0,065 -0,097 -0,029 -0,032 0,035 0,068 4 -0,058 -0,126 -0,039 -0,068 0,019 0,087 5 -0,071 -0,132 -0,039 -0,061 0,032 0,094 6 -0,065 -0,126 -0,039 -0,061 0,026 0,087 7 -0,061 -0,119 -0,039 -0,058 0,023 0,081 8 -0,055 -0,116 -0,039 -0,061 0,016 0,077 Stochastic dominance (6) ≽ N (5) (7) ≽ N (5) (8) ≽ N (5) (7) ≽ N (6) No (7) ≽ N (8)

First

  Order Stochastic Dominance and the Measurement of Hiring Discrimination: A ranking extension of correspondence testings with an application to gender and origin Emmanuel Duguet, Loic Du Parquet, Yannick L'Horty, Pascale Petit To cite this version: Emmanuel Duguet, Loic Du Parquet, Yannick L'Horty, Pascale Petit. First Order Stochastic Dominance and the Measurement of Hiring Discrimination: A ranking extension of correspondence testings with an application to gender and origin. 2012. halshs-00731005

Figure 1 :

 1 Figure 1: Interview decisions (gender case)

  measurement of discrimination, and this is what motivates our use of FOSD. Graphically, FOSD means that the CDF of candidate A -defined on ranks -stands above the CDF of candidate B. 3 We demonstrate the following property for the comparison of two candidates, A and B, in a testing with k candidates. Property: (a) A FOSD B ⇒ Pr A invited, B uninvited ≥ Pr B invited, A uninvited , or equivalently : A FOSD B ⇒ ∆ A, B ≥ 0.

FigureFigure 3 :

 3 Figure 2: Vietnamese origin women FOSD Vietnamese origin men

≽N

  Moroccan-origin men ≽ O Vietnamese-origin women ≽ N Senegalese-origin men ≽ O Vietnamese-origin men ≽ N Moroccan-origin women ≽ N Senegalese-origin womenThis ranking provides a position of the candidates along an ordinal discrimination scale, which is illustrated by Figure4. In order to simplify the reading, we have separated genders. We obtain three main results. First, the candidates can clearly be separated into several groups. Among men, those of French and Moroccan origins face the weakest discrimination, while those of Vietnamese and Senegalese origins face a strong degree of discrimination from recruiters. Second, among women, those of French and Vietnamese origins are clearly preferred, while African-origin candidates face the strongest discrimination. Second, women are more often discriminated against than men, which is illustrated by an upward shift along the discrimination ordinal scale compared to men.

Figure 4 :

 4 Figure 4: Discrimination ordinal scale

  (3) (1) ≽ N (4) (1) ≽ N (5) (1) ≽ N (6) (1) ≽ N (7) (1) ≽ N (8) (2) ≽ N (3) (2) ≽ N (4) (2) ≽ N (5) (2) ≽ N (6) (2) ≽ N (7) (2) ≽ N (8)

  given the origins of the candidates, and the last names associated with each origin are among the most widespread. The candidates live in inner Paris, and in similar districts in terms of geographic position and demographic composition. The individual characteristics of the 8 fictitious candidates are presented in Table1. In the case of the Vietnamese origin candidates, we had to indicate the gender explicitly on the CV, because the first names (Tien Hiep and Minh Trang) were not widespread enough to identify gender.

	CDF
	0.3
	0.2
			Table 1: Characteristics of the 8 candidates
			Candidates for job applications			Gender		Ethnic origin
	0.1		Anthony BERTRAND Place of residence: Paris 14 th district Sophie MOREAU Place of residence: Paris 14 th district			Male Female		French
	0.0	0	1 Abdallah ZALEGH 2 Place of residence: Paris 13 3 th district 4 Jamila KAIDI Place of residence: Paris 14 th district	5	Male Female	6	7 Moroccan	8	9
			Gender Place of residence: Paris 13 Amadou DIALLO th district Fatou DIOUF Place of residence: Paris 13 th district	Rank Man		Woman Male Female	Senegalese
			Tien Hiep PHAM Place of residence: Paris 13 th district Minh Trang TRAN Place of residence: Paris 13 th district	Male (indicated on CV) Female (indicated on CV)	Vietnamese

The first names given to the candidates are common,

Table 2 : Gender discrimination for each origin

 2 The Student statistics have been computed by the block bootstrap with 10000 replications. *: significant at the 5% level.

	Discrimination
	Senegal
	Morocco
		Vietnam
		Senegal
	Vietnamese
		Morocco
	French
		French
		Gender
	Women	Men

Table 3 : Origin discrimination for each gender

 3

	The Student statistics have been computed by the block bootstrap with 10000 replications.
	*: significant at the 5% level.
	Gender :	Men	Women
		Origin : Morocco
	No to both	65,8%	72,9%
	Yes to French origin, no to foreign origin	13,2%	16,8%
	Yes to foreign origin, no to French origin	6,5%	4,5%
	Yes to both	14,5%	5,8%
	∆	6,8%*	12,3%*
	Student	2,70	4,90
		Origin : Senegal
	No to both	67,7%	72,6%
	Yes to French origin, no to foreign origin	18,4%	19,0%
	Yes to foreign origin, no to French origin	4,5%	4,8%
	Yes to both	9,4%	3,5%
	∆	13,9%*	14,2%*
	Student	5,37	5,31
		Origin : Vietnam
	No to both	68,7%	65,8%
	Yes to French origin, no to foreign origin	19,0%	14,2%
	Yes to foreign origin, no to French origin	3,5%	11,6%
	Yes to both	8,7%	8,4%
	∆	15,5%*	2,6%
	Student	6,06	0,90

Table 4 : CDF of candidates' ranking

 4 The table gives Pr : C ≤ : , with : ∈ <1, … ,9=.

	Gender	Woman	Man	Woman	Man	Woman	Man	Woman	Man
	Origin	France	France	Morocco Morocco Senegal	Senegal	Vietnam Vietnam
	Rank	(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)
	1	0,158	0,094	0,013	0,058	0,010	0,026	0,061	0,029
	2	0,200	0,165	0,042	0,090	0,032	0,090	0,103	0,058
	3	0,213	0,203	0,065	0,142	0,045	0,110	0,142	0,074
	4	0,219	0,219	0,084	0,158	0,058	0,116	0,184	0,097
	5	0,219	0,239	0,100	0,187	0,065	0,135	0,197	0,103
	6	0,223	0,261	0,103	0,197	0,071	0,135	0,197	0,110
	7	0,226	0,274	0,103	0,210	0,077	0,139	0,197	0,116
	8	0,226	0,277	0,103	0,210	0,084	0,139	0,200	0,123
	9	1	1	1	1	1	1	1	1

The order statistic sorts the utility values in ascending order so that the utility vector , … , + becomes 8 , … , + 9 where is the smallest utility and + is the highest utility.

discriminated against. Second, candidates that do not seem to be discriminated against, because their discrimination coefficient is not significant although positive, can in fact face strong discrimination. Therefore combining strong and weak discrimination measures allows for a better assessment of the degree of discrimination that exists against the candidates.

Table 5A: Stochastic dominance analysis

The table presents the differences between the CDFs of the candidates given in Table 4. Negative numbers have been underlined. Reading: French origin candidates dominate all the other candidates, and do not dominate each other.

Reference :

Woman, France 8) 8)

Table 5C: Stochastic dominance analysis (continued from table 5B)

The table presents the differences between the CDFs of the candidates given in Table 4.

Appendix 1 -Experimental Design

The field of the study covers all offers of full-time jobs as a developer, with either fixed-term or permanent contracts, located in the Paris region. We tested all the job offers that came to our attention between the beginning of February 2009 and the beginning of April 2009. Overall, 310 job offers were tested, corresponding to the mailing of 2480 applications (8 x 310).

The 8 candidates differ only by their origin and gender. The CVs of all of the candidates were sent for each job offer. They indicated their age (25 years old), their nationality (French) and their family situation (single, no children) on their CVs. They followed a similar school and university curriculum: a scientific baccalauréat (A level), followed by a degree in computing, and finally a Master's degree in computing from one of the following Eventually, the CVs have been rotated.

Appendix 2 -Checking alphabetical ordering

We have investigated whether the ranking of the candidates on each job offer is compatible with alphabetical ordering. The results are presented in Table A1, and we find that this is not the case.