Laurent Michon
email: laurent.michon@univ-reunion.fr

Valérie Cayol

Ludovic Letourneur

Aline Peltier

Nicolas Villeneuve

Thomas Staudacher

Edifice growth, deformation and rift zone development in basaltic setting: Insights from Piton de la Fournaise shield volcano (Réunion Island)

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Basaltic shield volcanoes grow by the combined effect of endogenous and exogenous processes (e.g., [START_REF] Annen | The long-term growth of volcanic edifices: numerical modelling of the role of dyke intrusion and lava-flow emplacement[END_REF]. Endogenous growth is mainly related to dyke intrusions, which preferentially propagate along the existing rift zones. The volcano shape consequently depends on the intensity of their activity and may vary from elongated, e.g., Karthala in Grande Comore, and Mona Loa and Kilauea in Hawaii, to sub-circular, e.g., volcanoes of Western Galapagos. Exogenous growth results from the accumulation of lava flows and pyroclastic cones [START_REF] Naumann | Physical volcanology and structural development of Cerro Azul Volcano, Isabela Island, Galapagos: implications for the development of Galapagos-type shield volcanoes[END_REF][START_REF] Rowland | Slopes of oceanic basalt volcanoes[END_REF]. Long lava flows tend to flatten the volcano morphology, whereas short and thick ones may encourage the development of steep edifices [START_REF] Naumann | Physical volcanology and structural development of Cerro Azul Volcano, Isabela Island, Galapagos: implications for the development of Galapagos-type shield volcanoes[END_REF][START_REF] Annen | The long-term growth of volcanic edifices: numerical modelling of the role of dyke intrusion and lava-flow emplacement[END_REF].

At Piton de la Fournaise, Réunion Island, however, the magmatic activity led to the development of a steep central cone in the Enclos caldera despite the predominance of fissure eruptions that feed fluid basaltic lava flows. Moreover the cone is almost circular whereas most of the eruptions develop along the N25-30 and N120 rift zones (Fig. 1; Michon et al., 2007a). To answer these paradoxes, [START_REF] Annen | The long-term growth of volcanic edifices: numerical modelling of the role of dyke intrusion and lava-flow emplacement[END_REF] proposed, using numerical models, that the geometry of the cone was resulting from a predominant endogenous growth with only 30% of the edifice volume directly related to short lava flows. Although the models roughly fit the current cone's geometry, the input and natural data partly disagree. First, the simulated dykes (mean thickness of 3 m) are not as thick in nature as in the models. Indeed, recent works showed that the mean dyke thickness related to the post-1997 eruptions was around 0.5-0.7 m [START_REF] Fukushima | Transferts de magma au volcan du Piton de la Fournaise determinés par la modélisation 3D des données d'interférométrie radar entre 1998 et[END_REF][START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF]. Second, the model of [START_REF] Annen | The long-term growth of volcanic edifices: numerical modelling of the role of dyke intrusion and lava-flow emplacement[END_REF] simulates 10,000 dyke intrusions, 7% of them leading to an eruption. Such an intrusion/eruption ratio is the opposite of that observed since 1972, 81 eruptions (92%) and 7 intrusions (8%). Models with more realistic parameters show that the central cone should be much flatter [START_REF] Annen | The long-term growth of volcanic edifices: numerical modelling of the role of dyke intrusion and lava-flow emplacement[END_REF].

This paper aims at understanding the origin of the atypical geometry of the central cone of Piton de la Fournaise. It examines the morphology and inner structure of the cone. It also considers the role of the magmatic and volcano-tectonic structures in both the growth and the deformation of the cone, by combining field data and GPS data, as well as numerical models of rift dyke intrusions. Such a multi-disciplinary approach highlights the existing links between the deformation observed in the field, the displacement measured by GPS, and the recurrent intrusions of the dykes along preferential paths. The origin of the central cone and its deformation are interpreted in the light of both endogenous and exogenous growths.

Morphology and inner structure of the central cone

Piton de la Fournaise is the active volcano of Réunion Island (Fig. 1). The eruptive centre is located in the upper part of a large Ushaped caldera, the Enclos-Grand Brûlé structure, where most of the eruptions have occurred in the past 4.5 ka [START_REF] Bachèlery | Le Piton de la Fournaise (Ile de la Réunion). Etude volcanologique, structurale et pétrologique[END_REF]. The concentration of activity in the upper part of the caldera led to the construction of a 400-m-high cone. Magma intrusions, which originate from a magma chamber located at about sea level below the cone [START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF], are concentrated along two N25-30 and N120 rift zones (Fig. 1; Michon et al., 2007a). Outside the Enclos caldera, magma intrusions may continue along the NE and SE rift zones. The present-day cone summit shows two collapse structures: Bory crater in the West, which is currently inactive, and Dolomieu in the East, which experienced a caldera collapse during the large April 2007 eruption (Fig. 1; Michon et al., 2007b). Before this collapse, the elongated shape of the pre-existing Dolomieu was the result of the coalescence of several pit craters [START_REF] Lénat | Structure and dynamics of the central zone of Piton de la Fournaise volcano[END_REF][START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF]. As a whole, Bory and Dolomieu confer a global E-W elongation to the summit of the central cone.

The use of a 25-m step digital elevation model (DEM) provided by the Institut Géographique National allows a precise description of the limits of the cone, the distribution of the slopes and the structure of the summit zone. The cone's lower limit corresponds to an almost continuous sharp break-in-slope that separates the caldera floor and the steep flanks of the cone (Fig. 2). The local outlined discontinuities are due to the Puy Mi-Côte parasitic cone in the North, and the N120 topographic ridge in the SE, which is formed by the alignment of large pyroclastic cones along the N120 rift zone. These two structures disregarded, the base of the cone presents a sub-circular shape with a radius of about 1.6 km (Fig. 2a). The flanks of the cone present steep slopes, which contrast with the classical view of a basaltic effusive edifice. The western part of the cone, i.e., west of a N15-trending line centred on Dolomieu, is characterised by relatively homogeneous slopes ranging between 15°and 25°. East of the N15 boundary, the cone shows steeper slopes (25-30°), which locally reach 35°. Although the N120 topographic ridge gives to the southeast flank a complex slope distribution, the low-pass filtered DEM strikingly shows that the steep slopes of the cone's eastern half are not circumferentially distributed, but present two linear trends in the N150 and N55 directions (Fig. 2b). It is noteworthy that the N55 trend is aligned with a 500-m-wide and 700-m-long zone characterised by a network of parallel lineaments (Fig. 2c). Field observations and aerial photographs reveal that the lineaments correspond neither to the limit of lava flows, nor to eruptive fissures. Therefore, it is likely that the lineaments coincide to the surface trace of a fault zone, the activity of which may be related to the dynamics of the magmatic system. Finally, the summit morphology of the cone is characterised by present collapse structures, Bory and Dolomieu, surrounded by relatively flat areas which were interpreted as fully filled pit craters (Fig. 2b;Bachèlery,1981;[START_REF] Lénat | Structure and dynamics of the central zone of Piton de la Fournaise volcano[END_REF]. Past and present collapse structures are restricted to a slightly elongated zone, the centre of which is nearly superimposed to that of the 3.2 km-across circle marking the base of the cone (Fig. 2b).

The April 2007 caldera, which deeply cut the central cone, formed new scarps that reveal its internal structure (Fig. 3). The northern and southern caldera walls exhibit three main structural features:

(1) The northwestern part of the cone built up above a pre-Bory pit crater, the existence of which was already hypothesised by [START_REF] Lénat | Structure and dynamics of the central zone of Piton de la Fournaise volcano[END_REF] from the occurrence of gentle slopes north of Bory (Figs. 2b and3a). The limits of this collapse structure are underlined in the caldera wall by a strong lithological contrast characterised by thick lava units, and the succession of thin lava flows and scoria layers, inside and outside the pit respectively. (2) The volcanic units cut by the pre-Bory pit crater are subhorizontal and parallel to the present day surface, suggesting that the current geometry of the cone is inherited from a stage predating the development of the pre-Bory pit crater. This evolution disagrees with the historical reports of the 18th and 19th centuries in which the cone is considered to be symmetric and centred on the Bory pit crater until 1760, and becoming asymmetric when the activity shifted to the east flank of the cone in 1766 (Bory de Saint-Vincent, 1804).

(3) The caldera walls reveal a relatively constant lithological pile with a succession of scoria and scoria-rich units underneath a pile of thin lava flows (Fig. 3). The scoria materials suggest a period of predominantly explosive activity, which probably led to the construction of a large pyroclastic cone.

Fractures of the central cone

During the last decades, several works have addressed the origin of the tectonic structures of the central cone [START_REF] Bachèlery | Le Piton de la Fournaise (Ile de la Réunion). Etude volcanologique, structurale et pétrologique[END_REF][START_REF] Bachèlery | Caractères structuraux des eruptions historiques du Piton de la Fournaise (Ile de la Réunion)[END_REF][START_REF] Lénat | Structure and dynamics of the central zone of Piton de la Fournaise volcano[END_REF][START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF]. The different analyses, which were essentially or strictly restricted to the summit of the central cone, revealed the coexistence of concentric and radial volcano-tectonic structures, consisting mainly of extensional fractures and eruptive fissures, respectively. The concentric fractures are restricted within the first tens of meters from the rim of Dolomieu and Bory. Their formation is interpreted as related to the current and paleo pit crater development, and the dynamics of the superficial magma chamber. Only one fracture network which trends obliquely to the northeast rim of Dolomieu, is linked to another process, i.e., the slumping of the eastern flank [START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF]. The volcanotectonic structures described on the flanks of the central cone correspond to eruptive fissures, which usually present a radial strike and open en échelon during dyke injection. This distribution locally changes in the east flank where the eruptive fissures mimic conjugate systems [START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF]. The en échelon distribution on the northern and southern flanks of the central cone was initially interpreted as indicative of co-intrusive left-lateral and right-lateral displacements in the respective flanks, causing the eastward motion of the eastern flank of Piton de la Fournaise [START_REF] Bachèlery | Caractères structuraux des eruptions historiques du Piton de la Fournaise (Ile de la Réunion)[END_REF]. However, Michon et al. (2007a) recently showed that the en échelon pattern most likely resulted from a reorientation of the magma intrusion close to the surface due to the rotation of the main principle stress σ 1 from vertical to downslope, rather than a shear deformation during dyke intrusion.

Since the last structural work at the cone scale in 1990, Piton de la Fournaise underwent 27 eruptions and 3 intrusions, each of them leading to deformation of the summit zone. A new structural analysis was realised on the central cone and the adjacent Enclos floor in order to determine if tectonic structures developed and accommodated the deformation measured from GPS and interferometry [START_REF] Briole | Deformation at Piton de la Fournaise: evolution of the monitoring techniques and knowledge acquired in the last five years[END_REF][START_REF] Sigmundsson | Opening of an eruptive fissure and seaward displacement at Piton de la Fournaise volcano measured by RADARSAT satellite radar interferometry[END_REF][START_REF] Froger | The deformation field of the August 2003 eruption at Piton de la Fournaise, Reunion Island, mapped by ASAR interferometry[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF]. Our field observations reveal a wide range of fractures affecting the eruptive system. In a way similar to [START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF], the present work only considers the fractures cutting at least two adjacent geological units. Fractures restricted to single lava flows or cones, and related to their own dynamics, such as the syn-cooling subsidence or collapse, were disregarded. According to the fracture geometry, size and relationship with the geological formations, three different types of volcano-tectonic structures were distinguished, among which the two first correspond to the radial and concentric fractures already described by [START_REF] Bachèlery | Caractères structuraux des eruptions historiques du Piton de la Fournaise (Ile de la Réunion)[END_REF] and [START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF]:

(1) Eruptive fissures mainly open perpendicular to the slope (Fig. 4a).

Their mean length inferred from 68 post-1997 eruptive fissures is 180 m. Only few eruptive fissures are oblique to the slope. Their length is significantly higher (i.e., 2.5 to 11 times longer) than that of the downslope eruptive fissures. They are located in the SSE and east flanks of the cone, where they strike N20-25 and N65 respectively. In the east flank, the recent N65trending eruptive fissure of the August 2004 eruption is connected to the southeastern corner of Dolomieu. Aerial photographs of 1936 reveal that older eruptive fissures developed at the same location with a similar geometry. The geometric association of the N65 eruptive fissures with the radial ones forms the conjugate-like systems described by [START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF]. Occurrence of long and oblique eruptive fissures in the flanks suggests that besides the sub-surface downslope σ 1 , which controls the orientation of most of the eruptive fissures (Michon et al., 2007a), local stress fields or structures may influence the magma migration through the east flank. The local structural control is corroborated by the development of N80 and N50 eruptive fissures close to the southern rim of Dolomieu and northwest of Bory, respectively. Indeed, their orientation is parallel to the southern fault of Dolomieu and at the limit of the pre-Bory pit-crater. (2) Concentric fractures consist in tens to hundreds of metres long extensional fractures either parallel to the scarps of Bory and Dolomieu (Fig. 5a), or at the limit of paleo pit craters. Fracture width evolves from few centimetres to 1-2 m close to the limits of the pit-craters. The maximum fracture density is observed in Bory, southwest and north of Dolomieu and at the limit of Petit Plateau. [START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF] recently proposed that these concentric fractures result from annular extension related to pit-crater collapses. Analysis of the fractures before and after the April 2007 caldera reveals that most of the fractures form during the collapse events and are reactivated by radial extension stresses during the filling of the collapse structure (Michon et al., 2009-this issue).

(3) The third type of fracture corresponds to linear structures that usually affect several geological units, i.e., lava flows or pyroclastic cones. They differ from concentric fractures either around 450 and 640 m long in which fractures indicate an almost exclusive extension. Only a small part of these deformation zones corresponds to the N150 R' shear network described by [START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF]. Contrary to these authors, we believe that the shear indicators and the lava tube offsets, which suggest minor opposite lateral displacements, result from small differential movements between the lava blocks during a general extension rather than a succession of opposite senses of shear.

Co-intrusive deformation

Types of eruptions and the related magmatic paths

Activity of Piton de la Fournaise is characterised by fissure eruptions fed by a magma reservoir located at about sea level [START_REF] Nercessian | Internal structure of Piton de la Fournaise volcano from seismic wave propagation and earthquake distribution[END_REF][START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF]. Considering the location of the eruption site, three types of eruptions can be distinguished [START_REF] Peltier | Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion Island) inferred from deformation and geochemical data[END_REF]:

(1) Summit eruptions that start and remain in Dolomieu. (2) Proximal eruptions, which may start in the summit but progress to the flanks of the central cone and usually propagate downslope to the Enclos caldera floor. (3) Distal eruptions that develop away from the central cone, starting in the Enclos caldera floor. Numerical models suggest that the elevation of the first eruptive fissure, which opens during an eruption corresponds to the maximum elevation reached by the feeding dyke [START_REF] Fukushima | Transferts de magma au volcan du Piton de la Fournaise determinés par la modélisation 3D des données d'interférométrie radar entre 1998 et[END_REF][START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF]. Then, only the dykes related to summit and proximal eruptions (68% of the total eruptions) intrude into the central cone. We determined the dyke orientations related to the 25 proximal eruptions from the distribution of the related eruptive fissures for the period between 1981 and 2007 (Fig. 4b). Note that summit eruptions were disregarded since their dyke orientation cannot be accurately determined. Our compilation shows that 92% of the proximal dykes follow the N25-30 and N120 rift zones described by Michon et al. (2007a). 74% of them occur along the N25-30 axis with the same number of dyking events to the North and the South, and 26% of the eruptions developed on the SE fla n ka l o n gt h eN 1 2 0r i f tz o n e( Fig. 4b).

Co-intrusive displacements inferred from GPS

At Piton de la Fournaise, dyke intrusions usually lead to an asymmetric deformation of the volcano characterised by a concentration of displacements east of the dykes [START_REF] Lénat | The beginning of the 1985-1987 eruptive cycle at Piton de la Fournaise (La Réunion); New insights in the magmatic and volcano-tectonic systems[END_REF][START_REF] Sigmundsson | Opening of an eruptive fissure and seaward displacement at Piton de la Fournaise volcano measured by RADARSAT satellite radar interferometry[END_REF][START_REF] Froger | The deformation field of the August 2003 eruption at Piton de la Fournaise, Reunion Island, mapped by ASAR interferometry[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF]. GPS measurements carried out in 1981 and 1995 indicate that the eastern part of the cone summit experienced a 1.5-2-m eastward displacement during this period, whereas the western half remained stable [START_REF] Briole | Deformation at Piton de la Fournaise: evolution of the monitoring techniques and knowledge acquired in the last five years[END_REF]. A denser GPS network composed of about 80 stainless steel benchmarks cemented around the crater, on the flanks and at the base of the summit cone was installed by the Piton de la Fournaise Volcano Observatory in order to better constrain co-intrusive deformation of the edifice. Since 2002, this network was repeatedly measured after each eruption. The location of each benchmark was determined with respect to a reference receiver situated north of the Enclos caldera. Sampling-rate of both receivers was of 1 s and benchmark locations were recorded during a 7-min stop on each.

As demonstrated above, between 1981 and 2007, dykes intruding the central cone followed the N120 rift zone, and the northern and southern segments of the N25-30 rift zone, in roughly the same proportion. Between April 2003 and November 2003, 3 eruptions and 1 intrusion occurred at the summit and along each rift zone, making the co-intrusive deformation of this period representative of the general behaviour during summit and proximal eruptions. GPS data clearly show the deformation decoupling between the eastern and western parts of the cone (Fig. 7a). West of a N-S decoupling axis located between Bory and Dolomieu, displacements are small at the summit (b22 cm) and negligible at the base of the cone and in the Enclos (Table 1). East of the axis, displacements rotate from a northeastward to a south-eastward motion from north to south, showing an average eastward motion of the eastern half of the central cone. Total displacement values of the summit are three times larger than in the West. GPS data reveal that vertical displacements are restricted to the cone, reaching a maximum at the summit. In contrast, both the cone and the proximal part of the Enclos caldera floor experience significant outward horizontal displacements. Finally, it is interesting to note that the northern and western flanks present a steady decrease of displacement values, whereas the south-east and east flanks are characterised by a step wise decrease, correlated with variations of the topography. For the E flank, the GPS data located in the steepest part (data points 1a to 4a in Table 1a and Fig. 7a) indicate a slight decrease of vertical displacements, i.e., between 39.5 and 34.2 cm, whereas they sharply decrease from 24.2 to 3.3 cm over the same distance in the Enclos floor (data points 5a to 7a).

As a consequence, the recurrent dyke intrusions up to the summit and along the intra-cone segments of the N25-30 and N120 rift zones lead to heterogeneous growth of the central cone. 1c). Moreover, GPS data related to the January 2004 distal eruption indicate the existence of different displacement domains separated by the N55 deformation zone described above. The N55 fracture zone limits a differential uplift, which is twice as high west of the fracture zone as east.

In summary, our data clearly show two different co-intrusive modes of deformation. The cumulated deformation related to summit and proximal intrusions is centred on the summit craters and is mostly restricted to the central cone. In contrast, the distal intrusions led to a widespread deformation, the limits of which cannot be determined with the current GPS network.

Modelling of co-intrusive deformations

We use a three-dimensional boundary element method [START_REF] Cayol | 3D mixed boundary elements for elastostatic deformation field analysis[END_REF] to investigate: (1) the injection-related endogenous growth of the central cone; (2) the interaction between dyke injections and fractures; and (3) stress changes due to dyke injection along the northern and southern segments of the N25-30 rift zone. Models consider an elastic medium. Young's modulus is 5 GPa and Poisson ration is 0.25 [START_REF] Cayol | Three-dimensional modeling of the 1983-1984 eruption at Piton de la Fournaise volcano, Réunion Island[END_REF]. For the modelling, a mesh of the topography is constructed from a Digital Elevation Model of the volcano (Fig. 8). This mesh is dense close to the deformation sources and coarse further away. In order to limit errors caused by the finite extension of the ground surface mesh, its extension is chosen to be 5 times the source dimensions.

Three rift zone dykes were modelled, two along the N25-30 rift zone, to the North and to the South of the summit cone, and one along the N120 rift zone, east of the summit cone. Length of the dykes, heights and dips (Fig. 8) are typical of intrusions in rift zones [START_REF] Sigmundsson | Opening of an eruptive fissure and seaward displacement at Piton de la Fournaise volcano measured by RADARSAT satellite radar interferometry[END_REF][START_REF] Froger | The deformation field of the August 2003 eruption at Piton de la Fournaise, Reunion Island, mapped by ASAR interferometry[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF][START_REF] Peltier | Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion Island) inferred from deformation and geochemical data[END_REF]. The overpressure of the dykes is assumed to be 3 MPa, leading an average dyke opening of 2.3 m at the surface. As Piton de la Fournaise eruptions typically correspond to average dyke openings of 30 cm at the surface, these models are representative of more than 20 successive dyke intrusions, a number similar to the amount of intrusions that occurred between 1998 and 2007. Models show that dyke intrusions trigger an asymmetric deformation of the central cone (Fig. 9a). Flanks steepen east of the N25-30 rift zone of around 0.05°and do not significantly change in the cone's western half. Only the summit part located north of Dolomieu and Dolomieu itself flatten. It is worth noting that the N150-160 fracture zone, which shows predominant extensional deformation, is located at the transition between domains showing contrasting deformation, i.e. increasing and decreasing slopes (Fig. 9b).

We computed changes of Coulomb stresses on N55 vertical faults in order to investigate the influence of intrusions in the N25-30 and N120 rift zones on the N55 fractures network. Changes of Coulomb stress (ΔS) depend on normal and shear stress changes at given fault plane orientations for a given slip directions. They are defined by ΔS=Δτ-µ'Δσ [START_REF] King | Static stress changes and the triggering of earthquake[END_REF], where Δτ is the shear stress change on a given failure plane (positive in the direction of fault slip), Δσ is the change in normal stress on the plane (positive in compression) and µ', the apparent coefficient of friction, includes the effect of pore pressure changes. Here, we set µ' = 0. Coulomb stress changes on optimally oriented vertical planes were also calculated to determine whether a N55 fault zone could result from repeated rift zone intrusions. Following [START_REF] Cayol | Three-dimensional modeling of the 1983-1984 eruption at Piton de la Fournaise volcano, Réunion Island[END_REF],w e assumed a roughly isotropic regional stress in the horizontal plane. Models reveal that one kilometre east of Dolomieu the direction of the optimally oriented planes is consistent with the observed N55-65 lineaments (Fig. 11). It suggests that failure can be reached along N55-65 trending faults. However, the moderate Coulomb stress changes in the area of the N55-65 fracture zone suggest that stresses related to repeated dyke injections are probably not important enough to create faults.

Finally, computations were carried out to unravel principal stress variations due to repeated injections in the northern and southern segments of the N25-30 rift zone. Dyke propagation is controlled by the distribution of the total stress σ t , which corresponds to σ t =Δσ + σ r , where Δσ is the stress change due to magma injections and σ r is the pre-existing regional stress. The pre-existing regional stress being quantitatively unassessable, we solely determinedthestressvariationscausedby dyke intrusions. Fig. 12a reveals that minimum principal stresses changes, Δσ 3 , are subvertical and parallel to the rift zone, while intermediate principal stresses changes, Δσ 2 , are subhorizontal and parallel to t h er i f td y k e( Fig. 12a). [START_REF] Letourneur | The effects of rock heterogeneities on dyke paths and asymmetric ground deformation: the example of Piton de la Fournaise (Réunion Island)[END_REF] showed that the maximum principal stress σ r 1 of the regional stress is vertical. Consequently, the total minimum principal stress σ t 3 is probably horizontal and parallel to the N25-30 rift zone. As dykes intrude perpendicular to the total minimum principal stress, σ t 3 , the stress field induced by N25-30 rift injections promotes N120 dyke injections (Fig. 12b).

Discussion

Origin of the steep central cone

The central cone of Piton de la Fournaise presents abnormally steep flanks for a basaltic volcano characterised by fissure eruptions and fluid lava flows. [START_REF] Annen | The long-term growth of volcanic edifices: numerical modelling of the role of dyke intrusion and lava-flow emplacement[END_REF] proposed that this morphology mainly results from predominant endogenous processes with recurrent thick and short intrusions originating from a small and shallow magma chamber, i.e., 200 m in diameter and located 500 m below the summit. Besides the disagreement between the input and natural data concerning the intrusion vs. eruption ratio (see the description in Section 1), the existence of a small shallow magma chamber is questionable. Indeed, the present deformation and seismic data suggest that the feeding dykes originate from a relatively deep magma chamber located at about sea level [START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF][START_REF] Peltier | Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion Island) inferred from deformation and geochemical data[END_REF]. Furthermore, a small magma reservoir cannot explain the development of collapse structures such as the pre-Bory pit crater or the recent Dolomieu caldera, which are much larger than the hypothesized shallow magma chamber. Hence, we consider that the geometry of the central cone of Piton de la Fournaise does not result from the predominant endogenous processes described by [START_REF] Annen | The long-term growth of volcanic edifices: numerical modelling of the role of dyke intrusion and lava-flow emplacement[END_REF].

The April 2007 caldera collapse, which was coeval to the largest historical eruption deeply cut the central cone (Michon et al., 2007b;[START_REF] Urai | Size and volume evaluation of the caldera collapse on Piton de la Fournaise volcano during the April 2007 eruption using ASTER stereo imagery[END_REF]. Observations of the caldera walls reveal that the cone predating the pre-Bory pit crater was essentially made of red scoria. It was subsequently covered by several tens of metres of thin lava flows. We propose that the scoria formations are evidences of the occurrence of a hidden strombolian pyroclastic cone below the upper thin lava flows. Strombolian cinder cones are usually characterised by slopes between 25°to 30°, corresponding to the repose angle of the material [START_REF] Porter | Distribution, morphology and size frequency of cinder cones on Mauna Kea volcano, Hawaii[END_REF][START_REF] Wood | Morphometric evolution of cinder cones[END_REF]. At Piton de la Fournaise, the east and southeast flanks of the central cone present such slope values. However, our GPS data, which clearly show progressive slope tilting during the recurrent intracone intrusions, suggest that part of the present slope values result from deformation processes. In the cone's western half, where the intrusionrelated deformation is limited, slope values are 10°less than that of a classical cinder cone. Consequently, we interpret the present cone's geometry as resulting from a twofold evolution: a first period during which a cinder cone built up in the Enclos caldera floor, and a second period of predominant effusive activity overlapping and smoothing the initial cone geometry. A similar scenario was proposed to explain the steep morphology of Nyiragongo (Democratic Republic of Congo), where an initial explosive phase built up a steep cone and a lava lake activity with periodic overflowing veneered over the pyroclastic core [START_REF] Mcdonald | Volcanoes[END_REF][START_REF] Demant | Volcanological and petrological evolution of Nyiragongo volcano, Virunga volcanic field, Zaire[END_REF].

Tectonic structures, co-intrusive deformation and rift zones

Since 1981, all the structural analyses and geodetic measurements have shown that the eastern part of the central cone of Piton de la Fournaise was affected by an eastward motion [START_REF] Bachèlery | Le Piton de la Fournaise (Ile de la Réunion). Etude volcanologique, structurale et pétrologique[END_REF][START_REF] Bachèlery | Caractères structuraux des eruptions historiques du Piton de la Fournaise (Ile de la Réunion)[END_REF][START_REF] Delorme | March 1986 eruptive episodes at Piton de la Fournaise volcano (Reunion Island)[END_REF][START_REF] Lénat | The beginning of the 1985-1987 eruptive cycle at Piton de la Fournaise (La Réunion); New insights in the magmatic and volcano-tectonic systems[END_REF][START_REF] Zlotnicki | Eruptive mechanisms on Piton de la Fournaise volcano associed with the December 4, 1983, and January 18, 1984, eruptions from ground deformation monitoring and photogrammetric surveys[END_REF][START_REF] Briole | Deformation at Piton de la Fournaise: evolution of the monitoring techniques and knowledge acquired in the last five years[END_REF][START_REF] Sigmundsson | Opening of an eruptive fissure and seaward displacement at Piton de la Fournaise volcano measured by RADARSAT satellite radar interferometry[END_REF][START_REF] Froger | The deformation field of the August 2003 eruption at Piton de la Fournaise, Reunion Island, mapped by ASAR interferometry[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF][START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF][START_REF] Tinard | Caractérisation et modélisation des déplacements du sol associés à l'activité volcanique du Piton de la Fournaise, île de La Réunion, à partir de données interférométriques[END_REF]. It has been first proposed that the entire eastern flank was involved in the deformation process (e.g., [START_REF] Bachèlery | Le Piton de la Fournaise (Ile de la Réunion). Etude volcanologique, structurale et pétrologique[END_REF][START_REF] Lénat | The beginning of the 1985-1987 eruptive cycle at Piton de la Fournaise (La Réunion); New insights in the magmatic and volcano-tectonic systems[END_REF]. The NE and SE rift zones were then considered as sinistral and dextral shear zones accommodating the lateral displacements of the volcano flank. However, recent GPS and radar interferometry data do not confirm such a general process and reveal that the deformation is directly related to dyke intrusions and concentrated around the central cone [START_REF] Sigmundsson | Opening of an eruptive fissure and seaward displacement at Piton de la Fournaise volcano measured by RADARSAT satellite radar interferometry[END_REF][START_REF] Froger | The deformation field of the August 2003 eruption at Piton de la Fournaise, Reunion Island, mapped by ASAR interferometry[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF]. If the deformation pattern is well constrained with these types of data, little is known about the effect of recurrent intrusions in the development, or interaction with tectonic structures. Using field data collected around the summit collapse structures, i.e., the break-in-slope in the upper eastern flank, the fractures in the northern part of the N150-160 fracture zone and the conjugate-like geometry of the eruptive fissures, [START_REF] Carter | Pits, rifts and slumps: the summit structure of Piton de la Fournaise[END_REF] proposed an additional deformation process. According to these authors, the east flank of the central cone undergoes a progressive slump above a low strength layer, which is favoured by the recurrent intrusions along the N10 and N170 rift zones. We disagree with this model for the following reasons: (1) the fractures used to interpret the northern limit of the slump as a dextral shear zone do indicate a predominant extension, and correspond to a small part of a much longer extension structure, i.e., the N150-160 fracture zone; (2) the conjugate-like geometry of eruptive fissures results from two sets of eruptive fissures that are the common radial fissures and the long N55 trending fissures which develop in the east flank only; (3) the eastern base of the cone lacks any thrust fault or concentric compressional structures which would result from the slump; (4) The summit collapses disregarded, the deformation of the cone can be exclusively due to the dyke intrusions [START_REF] Froger | The deformation field of the August 2003 eruption at Piton de la Fournaise, Reunion Island, mapped by ASAR interferometry[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF] and to a minor pre-eruptive inflation centred on thecone [START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF][START_REF] Peltier | Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion Island) inferred from deformation and geochemical data[END_REF].

This brief review reveals a lack of any global model in which both field observations and GPS or radar interferometry data are integrated. In the following, we propose to combine our multi-disciplinary data in order to determine: (i) the origin of the different tectonic structures (except for the summit collapses and the concentric fractures that are clearly linked to the magma withdrawal; [START_REF] Hirn | Episodes of pit-crater collapse documented by seismology at Piton de la Fournaise[END_REF][START_REF] Longpré | The November 2002 eruption at Piton de la Fournaise volcano, La Réunion Island: ground deformation, seismicity, and pit crater collapse[END_REF]Michon et al., 2007b;Michon et al., 2009-this issue); (ii) the role of the recurrent dyke intrusions in the morphology of the central cone; and (iii) the origin of the N25-30 and N120 rift zones and their consequences in the cone's evolution.

Our structural analysis of the central cone allows identification of two sets of tectonic structures. First, the different spatial distributions of flank fractures, eruptive fissures and dykes, especially in the western part of the cone (Figs. 4 and6a), confirm that the development of the flank fractures does not correspond to the dyke-induced proximal deformation. In contrast, the superposition of the areas affected by these fractures and the co-intrusive deformation (Figs. 6a, 7a and 9a) suggests that the flank fractures are linked to the outward inflating tilting of the eastern part of the cone. Hence, we propose that the flank fractures, which are essentially characterised by extension, correspond to relatively shallow (maximum depth of few hundreds of metres) tension cracks that accommodate the endogenous growth during summit and proximal injections. The N20-30-trending fractures may result from the recurrent intrusions along the N25-30 rift zone whereas the N150-160 fractures may originate from the combined effect of intrusions along the N120 rift zone and the northern segment of the N25-30 rift zone (Fig. 9b andc).

The N55 fracture zone represents the second set of fractures. It is underlined by: (i) long, oblique eruptive fissures in the east flank of the cone (Fig. 4); (ii) a dense fracture network (Fig. 6c); (iii) striking lineaments on the DEM (Fig. 2c); and (iv) linear very steep slopes in the SE flank (Fig. 2b). As a whole, it corresponds to a 3-km-long structure, which is parallel to a larger tectonic structure that accommodated the collapse of the southern part of the Enclos caldera [START_REF] Michon | The morphology of Piton de la Fournaise basaltic shield volcano (La Réunion island): characterization and implication in the volcano evolution[END_REF]. Given the characteristics of the N55 fracture zone and the occurrence of volcano-tectonic events along or close to the lineaments [START_REF] Sapin | Stress, failure and fluid flow deduced from earthquakes accompanying eruptions at Piton de la Fournaise volcano[END_REF], we interpret it as a fault zone. The direction of the eruptive fissures in the eastern flank suggests that the fault zone controls the magma intrusion in the eastern part of the cone. Numerical models show that the N55-65 fault zone may be punctually reactivated by stresses related to dyke intrusions along the rift zones. They also indicate that the development of the fault zone does probably not result from intrusion-related stresses. Hence, we propose that the N55-65 fault zone is a tectonic structure inherited from the collapse of the Enclos Fouqué caldera, which is reactivated by the present-day dynamics of the magmatic system.

Our numerical models and GPS data indicate that the recurrent intrusions along the rift zones induce a heterogeneous deformation of the cone with (i) a significant inflation of the eastern part of the summit, (ii) the steepening of the southeast and east flanks, and (iii) very slight tilting of the western part (Fig. 9a). The superposition of the tilted zones and the N55 and N150 very steep slope zones strongly suggests that the slope value difference of about 5-10°between the east and southeast flanks, and the west and north flanks results from the endogenous growth associated with intra-cone intrusions along the N25-30 and N120 rift zones (Fig. 9). Hence, if most of the present morphology of the central cone results from an exogenous growth, the specific morphology of the cone's eastern part is clearly due to endogenous processes.

The absence of any elongation along the N25-30 rift zone despite the substantial concentration of the intrusions along this structure contrasts with the morphologies of Hawaiian volcanoes and Karthala. [START_REF] Annen | The long-term growth of volcanic edifices: numerical modelling of the role of dyke intrusion and lava-flow emplacement[END_REF] showed that the geometry of the volcanoes was sensitive to the dyke length/source depth ratio. Volcanoes with ratios lower than 1 present an elongation perpendicular to the rift zones, whereas the elongation is parallel to the rift zones for ratios greater than 1. For ratios close to one, the edifice remains circular despite the presence of radial rift zones. At Piton de la Fournaise, the mean dyke length along the northern and southern segments of the N25-30 rift zone and the source depth are of the same range (about 2000 m). Then, the development of short dykes along the N25-30 and N120 rift zones may explain the subcircular geometry of the central cone.

The N25-30 and N120 rift zones were determined from the orientation and location of the dykes intruded since 1981. Michon et al. (2007a) proposed that the N25-30 rift zone, on which the cone built up, results from recurrent magma intrusions along a large, 10 kmlong structure in the Enclos caldera floor. The origin of the N120 rift zone is less constrained. It could be related to the regional N120-trending volcanic zone between Piton des Neiges and Piton de la Fournaise, the development of which has been interpreted as to be controlled by N120 crustal faults (Fig. 1; [START_REF] Chevalier | Volcano-structural evolution of Piton des Neiges, Reunion Island, Indian Ocean[END_REF][START_REF] Michel | Self-potential and magnetic surveying of La Fournaise volcano (Réunion Island: correlations with faulting, fluid circulation, and eruption[END_REF]Michon et al., 2007a). However, the lack of any N120trending intrusions in the Enclos caldera, west of the N25-30 rift zone, and the contrasted geometries between the narrow rift zone east of the cone and the wide volcanic zone between Piton de la Fournaise and Piton des Neiges suggest a different origin of both structures. Our numerical modelling reveals that recurrent intrusions along the northern and southern segments of the N25-30 rift zone encourage intrusions in the N120 direction (Fig. 12a). Hence, we explain both the existence of the N120 rift zone East of the N25-30 rift zone and the lack of N120 trending dykes West of it, as a consequence of the general eastward dip of the N25-30 dykes (Fig. 12b). Such a geometry, which has been deduced from several sets of deformation data [START_REF] Sigmundsson | Opening of an eruptive fissure and seaward displacement at Piton de la Fournaise volcano measured by RADARSAT satellite radar interferometry[END_REF][START_REF] Froger | The deformation field of the August 2003 eruption at Piton de la Fournaise, Reunion Island, mapped by ASAR interferometry[END_REF][START_REF] Fukushima | Finding realistic dike models from interferometric synthetic aperture radar data: The February 2000 eruption at Piton de la Fournaise[END_REF][START_REF] Peltier | Constraints on magma transfers and structures involved in the 2003 activity at Piton de La Fournaise from displacement data[END_REF][START_REF] Peltier | Cyclic magma storages and transfers at Piton de La Fournaise volcano (La Réunion Island) inferred from deformation and geochemical data[END_REF], induces asymmetric deformations associated with much larger stress changes East of the N25-30 rift zone than West. The small length of the N120 rift zone may originate from a rapid decrease of the stress increase due to the N25-30 intrusions away from the plumbing system. This rift system, which has been determined from the summit and proximal intrusions is not fully connected with the NE and SE rift zones outside the Enclos caldera. Indeed, recent geophysical data [START_REF] Brenguier | 3-D surface wave tomography of the Piton de la Fournaise volcano using seismic noise correlations[END_REF] suggest that the NE and SE rift zones along which the distal eruptions occur are related to tectonic structures that are located below the magma chamber that feeds the summit and proximal eruptions.

Conclusion

Our work corresponds to the first analysis of the central cone of Piton de la Fournaise, which combines structural and GPS data, and numerical modelling. The main results can be summarised as follows:

(1) Similarly to Nyiragongo [START_REF] Demant | Volcanological and petrological evolution of Nyiragongo volcano, Virunga volcanic field, Zaire[END_REF], the central cone is formed by a pyroclastic core overlapped by few tens of metres of lava flows. The general steep geometry of the cone does consequently not result from endogenous processes but from an initial phase during which a pyroclastic cone formed.

(2) The cone undergoes a contrasted co-intrusive deformation during summit and proximal eruptions along two N25-30 and N120 rift zones. Most of the co-intrusive displacements are concentrated in its eastern half where they induce a progressive outward inflation and the tilting of the eastern and southeastern flanks. A dense network of flank fractures, which differs from the eruptive fissures and the concentric fractures linked to the summit collapses, accommodate the progressive endogenous growth.

(3) The dynamics of the plumbing system also reactivates a N55 fault zone, which mainly controls magma intrusions in the eastern flank and decouples the deformation of the Enclos caldera floor during the north-eastward distal eruptions. (4) The N25-30 rift zone is related to a large, 10-km-long fracture zone, whereas the short, narrow N120 rift zone results from a N120-trending stress accumulation west of the N25-30 rift zone, which guides magma intrusions.

Fig. 1 .

 1 Fig. 1. 25-m step Digital Elevation Model of Piton de la Fournaise. The edifice is characterised by two NE and SE rift zones outside the Enclos caldera. In the Enclos caldera, the recent eruptions are concentrated along two N25-30 and N120 rift zones (Michon et al., 2007a). The central cone is cut by two summit collapse structures, Bory and Dolomieu. Coordinates in UTM WGS84 (zone 40S).

Fig. 2 .

 2 Fig. 2. (a) Slope map of the central cone calculated from the 25-m step DEM. The cone, which is sub-circular with a radius of 1.6 km, is coalescent to a parasitic cone in the north and a N120 topographic ridge in the NE, formed by several pyroclastic cones. The western and eastern parts of the cone present different slope values. (b) Slope map of the calculated from a low-pass filtered DEM in order to determine the main characteristics of the morphology. The very steep slopes of the cone's eastern half are distributed in two linear steep slope zones trending in the N55 and N150 directions. PP: Petit Plateau; PBPC: Pre Bory pit crater. c-Shaded relief representation of a 7-m step DEM calculated by stereophotogrammetry from a set of aerial photographs performed in 1989. This high-resolution DEM highlights clear, parallel lineaments in the northeastern continuity of the N55 steep slope zone. Coordinates in UTM WGS84 (zone 40S).

Fig. 3 .

 3 Fig. 3. Panoramas and interpretations of the northern (a) and southern (b) scarps of the April 2007 caldera (photo b: Frederic Massin).

Fig. 4 .

 4 Fig. 4. (a) Distribution of the eruptive fissures. (b) List and location of the dykes associated to proximal eruptions since 1981. Coordinates in UTM WGS84 (zone 40S).

Fig. 5 .

 5 Fig. 5. (a) Picture of the concentric fractures associated to the collapse structures. (b) Pictures of the flank fractures that affect the eastern half of the central cone. Here, the N150 deformation zone.

  6. (a) Spatial distribution of the flank fractures. (b) Rose diagram illustrating the three main fracture trends: N20-30, N55-65 and N150-160. (c) Spatial distribution of the three main fracture trends. Coordinates in UTM WGS84 (zone 40S).

  Between April and November 2003, the cone underwent a very slight steepening of the western flank (≈4×10 -3 °/yr), and a significant inflation of the eastern part, leading to a steepening of the south and east flanks (≈1.5 × 10 -2 °/ yr). Moreover, GPS data suggest that each N55 and N150 trending steep slope zone undergoes a homogeneous deformation. This overall deformation pattern inferred for the period between April 2003 and November 2003 fully agrees with the slope change determined by Fukushima (2005) between March 1998 and June 2000, during which 5 eruptions occurred along the S and N segments of the N25-30 rift zone and along the N120 rift zone. The striking similarities between the two different periods suggest a relatively constant mode of deformation of the cone when summit and proximal eruptions occur. Most of the deformation determined for the period between March 2002 and January 2006 is similar to that of April to November 2003 in terms of (i) decoupling between the western and eastern parts and (ii) and progressive decrease of displacement amplitudes from the top to thebaseoftheconeinthenorth,southwestandwestflanks (Fig. 7b). It only differs in one point, which is the regular increase of total displacement values from 146.5 to 180.4 cm between the eastern part o ft h es u m m i ta n dt h eE n c l o sfloor (Table 1b). Contrary to the periods between March 1998 and June 2000 and between April 2003 and November 2003, during which only summit and proximal eruptions occurred, the period between March 2002 and January 2006 was characterised by summit, proximal and also 3 distal eruptions in the Plaine des Osmondes (see Fig. 1 for location). The displacement pattern related to one of these eruptions, in January 2004, suggests that the progressive increase of displacements when going east of the summit is induced by distal intrusions toward the Plaine des Osmondes (Fig. 7cand Table

  GPS for periods between April and November 2003 (a), March 2002 and January 2006 (b) and caused by the magma migration linked to the January 2004 distal eruption (c).See Fig.8for the location of the data points.

Table 1

 1 Displacements measured by GPS for periods between April and November 2003 (a), March 2002 and January 2006 (b) and caused by the magma migration linked to the January 2004 distal eruption (c).See Fig.8for the location of the data points.

Fig. 7 .

 7 Fig. 7. Horizontal (left) and vertical (right) co-intrusive displacements recorded by GPS. (a) Between April and November 2003. The dashed lines indicate the location of the dykes. (b) Between March 2002 and January 2006. (c) During the January 2004 eruption. Coordinates in UTM WGS84 (zone 40S).

  4, corresponding to laboratory values and small fluid pressure. Positive values of Coulomb stress changes bring faults closer to failure whereas negative values inhibit failure. We considered Coulomb stress changes induced by single dyke intrusion in each segment of the rift zones and by combined intrusions (Fig. 10). Dyke injections in the northern segment of the N25-30 rift zone promote right-lateral fault slip in the area of the N55-65 fractures, whereas injections in the southern segment of the N25-30 rift zone and in the N120 rift zone encourage left-lateral slip in the N55-65 fracture zone. Combinations of injections in the northern and southern segments of the N25-30 rift zone as well as in the N120 rift zone favour right lateral slip of N55 faults located 1 km East of Dolomieu crater.

Fig. 8 .

 8 Fig. 8. Map view, E-W and N-S cross sections, and perspective view of the Boundary Element Model computed to test the influence of repeated dyke injections into the N120 and N25-30 rift zones. All dykes are 1000 m high. The N120 dyke is represented by a blue mesh. It is vertical and 2.5 km long. The northern segment is represented by a red mesh. It dips 70°to the East and is 2.3 km long. The southern segment is represented by a green mesh. It dips 70°to the East and is 1.2 km long. The perspective view only shows a portion of the surface of the topography mesh, which extends five times further than shown.

Fig. 9 .

 9 Fig. 9. (a) Modelled slope changes due to repeated dyke intrusions into the N25-30 and N120 rift zones (see Fig. 8 for the dyke intrusion geometries). Dyke surface traces are represented with a white thick solid line. Thin lines represent elevation contours at 100 m intervals. (b) Slope map suggesting the effect of dyke intrusions in the development of the eastern and south-eastern steep flanks. (c) E-W cross section presenting the outward inflation of the cone's eastern half.

Fig. 10 .

 10 Fig. 10. Coulomb stress changes (Mpa) at a 1500 m elevation induced on N55 vertical faults by intrusions in the N25-30 and N120 rift zones (see Fig. 8 for the dyke intrusion geometries). Stress changes for right-lateral (left column) and left lateral (right column) slip are compared. Surface traces of the intruded dykes are shown as bold white lines. Null Coulomb stress changes are indicated by a light black line.

Fig. 11 .

 11 Fig. 11. Coulomb stresses changes (MPa) at a 1500 m elevation on optimally oriented vertical faults induced by repeated dyke intrusions in the N25-30 and N120 rift zones (see Fig. 8 for dyke intrusion geometries). Surface traces of the intruded dykes are shown as bold white lines. The directions of optimal slip planes are indicated in black and white, where the black and white colours indicate opposite slip sense.

Fig. 12 .

 12 Fig. 12. (a) Principal stress changes Δσ corresponding to repeated magma injections in the northern and southern segments of the N25-30 rift zone (see Fig. 8 for dyke intrusion geometries). Δσ 1 , Δσ 2 and Δσ 3 correspond to the maximum, intermediate and minimum principal stress changes, respectively. (b) Interpretative sketch of the short, narrow N120 rift zone as resulting from the stress accumulation induced by recurrent intrusion along the N25-30 rift zone.

Acknowledgments

The authors thank Ruth Andrew for improving an initial version of the manuscript. Two anonymous reviewers are kindly thanked for their constructive reviews, which helped to significantly improve the manuscript. This is IPGP contribution 2437.


W Román Acevedo

S Di Napoli

F Romano

G Rodríguez Ruiz

P Nukala

C Quinteros

J Lecourt

U Lüders

V Vildosola

D Rubi

Polaron formation in Bi-deficient BaBiO 3

Keywords:

BaBiO3 is a charged ordered Peierls-like perovskite well known for its superconducting properties upon K or Pb doping. We present a study on the transport and electronic properties of BaBiO3 perovskite with strong Bi-deficiency. We show that it is possible to synthesize BaBiO3 thin layers with Bi-vacancies above 8-10% by depositing an yttrium-stabilized zirconia capping layer. By combining transport measurements with ab initio calculations we propose an scenario where the Bi-vacancies give rise to the formation of polarons and suggest that the electrical transport is dominated by the migration of these polarons trapped at Bi 3+ sites. Our work shows that cation vacancies engineering -hardly explored to date-appears as a promising pathway to tune the electronic and functional properties of perovskites.

I. INTRODUCTION

Perovskite oxides display a plethora of physical properties, such as magnetism 1 , ferroelectricity 2 and/or superconductivity 3 , making them attractive materials for the development of sensing 4 , information storage 5 or computing [START_REF] Yu | Neuro-inspiring computing using resistive synaptic devices[END_REF] devices, among others. The possibility of implementing these devices has boosted in the last decades as oxide thin films and heterostructures with excellent crystalline properties, chemistry and atomically controlled surfaces and interfaces have been developed. Defects are ubiquitously found in oxide thin films, oxygen vacancies being the most common ones [START_REF] Gunkel | [END_REF] . Cationic vacancies have been less studied due to their lower concentration in comparison with oxygen vacancies 8 , but they might have a significant impact on the physical properties of the films as they could strongly modify the oxide electronic structure. BaBiO 3 (BBO) is a perovskite which becomes superconducting upon Pb or K doping 9,10 with high critical temperatures up to ∼ 30 K. Pristine BBO is an insulator due to the existence of a charge ordered state comprising Bi ions with formal +3 and +5 valences 11 . This is accompanied with structural distortions involving both the so-called breathing distortions and tiltings of BiO 6 octhaedra 12 , which reduces the high temperature cubic crystal symmetry to a monoclinic structure at room temperature. Bi-O bond distances have been reported as 2.11 and 2.29 Å, respectively 12 . Spectroscopic experiments have failed to resolve valence differences between Bi ions in both crystallographic sites 13,14 , while theoretical work suggests a common ≈ +3 state, which has been rationalized in terms of spatially extended Bi-O bonds related to hybridized Bi-6s and O-2p atomic orbitals [15][16][17] . BBO has been proposed to behave as a topological insulator upon electron doping 18 , and it was also theoretically sug-gested to present a 2D electron gas in its Bi-terminated (001) surface 19 . Moreover, 2D superconductivity has been shown at the BBO/BaPbO 3 interface 20,21 . Reports on BBO thin films are scarce [22][23][24][25][26][27][28][29][30][31] . Most works report on the growth and properties of BBO on SrTiO 3 (STO) substrates 22,23,25,[28][29][30][31] , with both materials presenting a rather large lattice mismatch of ≈ 12 %. It was shown that this misfit is accommodated through the formation of an ultrathin interfacial δ-Bi 2 O 3 layer with a fluorite structure, which forms a coherent interface with the substrate and a semi-coherent interface with BBO 28,31 . Recent reports indicate that the electronic structure of BBO on STO is thickness dependent and is affected by chemical and structural effects imposed by the substrate/film interface 29,32 .

These results indicate that substrate effects on the BBO electronic structure might not be disregarded in thin epitaxial films.

With the aim of studying the effect of (B-site) Bi deficiency on the electronic and transport properties of BBO, we have developed a fabrication procedure allowing for the introduction of a large amount of Bi vacancies in stochiometric BBO thin films on silicon. From a combination of experiments and ab inito calculations, it is proposed that Bi deficiency leads to the formation of polarons that dominate the electrical transport. Our work shows the control of cation deficiency as an emerging degree of freedom that allows tuning the BBO electronic structure and its functional properties.

II. EXPERIMENTAL SECTION

BBO thin films were deposited by pulsed laser deposition using a Nd-YAG Spectra Physics laser with λ = 266 nm and 850 mJ per pulse. Structural characterization of the films was performed by means of X-ray diffraction (XRD) using two Panalytical Empyrean diffractometers with PIXcel ultrafast detectors, in standard Bragg-Brentano configuration (0.02º steps were used for 2θ scanning). Transmission electron microscopy (TEM) was performed with a Thermo Fisher Scientific Themis Z monochromated and double corrected microscope operated at 300 kV. In order to get Bi-deficient BBO, we followed this procedure: first, we deposited a stochiometric BBO film (≈ 165 nm thick) on commercial Si/SiO 2 substrate, at a temperature and oxygen pressure of 550 ºC and 0.01 mbar, repectively. We chose this substrate in order to obtain strain-free polycrystalline films, as will be shown below. Next, we deposited (in-situ), on top of BBO, a yttrium stabilized zirconia (YSZ, ≈ 17 nm thick) capping layer, at room temperature, at an optimized oxygen pressure of 0.1 mbar. YSZ is a well known oxygen conductor 33 but has also been reported to display high cationic conductivities 34 . In our case, under the proper deposition conditions, YSZ capping favours the diffusion of Bi from BBO to YSZ, leading to an interfacial BBO layer with a strong Bi deficiency. Figure 1(a) shows X-ray diffraction spectra recorded for both bare and capped BBO films. It is found in both cases the presence of peaks corresponding to the BBO perovskite structure while, for capped BBO, YSZ remains amorphous. BBO peaks were indexed considering a pseudo-cubic (PC) lattice. We extracted a cell parameter for bare BBO a P C = 4.339(5) Å, in excellent agreement with the value corresponding to bulk BBO (4.335 Å). Upon YSZ capping, there is a shift of X-ray peaks towards lower angles (see Figure 1(b)), reflecting an (average) BBO unit cell parameter expansion of ≈ 0.5%. This is a first indication of the formation of Bi vacancies at the BBO layer, as unit cell enlargement is a signature of B-site vacancies in perovskites [36][37][38] . We will come back later to this issue.

Figure 1(e) displays a HAADF-STEM cross-section of a capped BBO film. First we notice the existence of a thin BaO layer (≈ 11 nm, non detectable by XRD) at the substrate/BBO interface, confirmed by the EDS maps displayed at the Suppl. Information 35 . This interfacial BaO layer has been observed before in BBO thin films 27 and can be related to the poor sticking coefficient of Bi ions during the first growth stages 29,39 . Next, Figure 1(c) displays a BBO layer divided in two zones with different contrast, the darker one (≈ 30 nm thick) located at the interface with YSZ. As HAADF brightness is higher for species with high atomic number such as Bi, the darker BBO fringe is an indication of Bi-deficiency at this zone. We notice that it is difficult to perform a precise composition quantification, as during EDS mapping we observed some grain growth together with changes in the HAADF contrast between BBO and Bi-deficient BBO 35 , indicating atomic diffusion during the measurement. With this constraint in mind, we can roughly estimate a lower bound of ≈ 8-10% of Bi vacancies at the Bi-deficient BBO layer 35 . This statement is supported from the analysis of (200) P C X-ray peaks performed in the Suppl. Information 35 . The (200) P C reflection of bare BBO can be fitted with a single Gaussian component, centered at 2θ = 41.59 º, corresponding to a bulk-like cell parameter 4.339(5) Å. On the other hand, the capped BBO (200) P C reflection shows a broadened peak with a slight asymmetry. The peak can be deconvoluted by assuming two Gaussian contributions centered at 2θ = 41.59 º and 2θ = 41.25 º. These two components are attributed to BBO and Bi-deficient BBO, with cell parameters 4.339(5) Å and 4.374(5) Å, respectively. The latter cell parameter accounts for a ≈ 4 pm BBO cell expansion due to the presence of Bi-vacancies. We recall that in the case of Ti deficient SrTiO 3 such expansion corresponds to a 10-15% B-cation deficiency 36 , allowing an additional estimation of the amount of B-site vacancies in the Bi-deficient layer of capped BBO, which is consistent with the previous one obtained from microscopy analysis. Now we will focus on the electrical transport properties of our samples. BBO single crystals were reported to display an insulating, thermally activated behavior 40 . Transport properties in BBO thin films have been rarely published: some authors reported an "insulating behavior" 24 in BBO thin films on MgO but no temperature dependent resistivity was shown, while a negative resistance coeficent, in the range 200-300 K, was displayed for BBO films on silicon 26 . We recall that the transport properties in semiconducting oxide thin films can be highly affected by the presence of defects. For instance, BaSnO 3 thin films were shown to present dramatic changes in their resistivity values -and the corresponding temperature behavior-depending on the films' free carriers concentration 41 , strongly linked to the type and number of existing defects. Electrical transport has been measured on the deposited structures at low temperatures (50-200 K) in a 4-point configuration using a multipurpose Quantum Design Versalab cryostat. In order to gain electrical access to the capped BBO samples, we have covered some parts of the substrate by shadows masks during the deposition (Figure 2(b)). Electrical contacts were fixed with silver paste. Figure 2(c) displays the evolution of the 4points resistance versus temperature for both uncapped and capped BBO films. Bare BBO films show a resistance nearly independent of the temperature -in the measured temperature range-, indicating that the expected thermally activated behavior (negative resistance coefficient) is compensated with a reduction of carriers mobility as temperature increases. We note that the mentioned behavior of the mobility with temperature implies that mobility is limited by phonons rather than by impurity scattering, this balance being very sensitive to free carriers density 42 . On the other hand, capped BBO shows lower resistance values -a 6-fold decrease at 175 K-together with a semiconducting-like behavior characterized by a decrease of the resistance as the temperature increases. We checked in separated control experiments that both SiO 2 and YSZ present resistances several orders of magnitude higher than those shown in Figure 2(c) and this is also the case for the estimated resistance corresponding to the BaO layer (> 10 GΩ), given its reported resistivity (≈ 10 4 Ωcm at ≈ 700 ºC) and semiconducting-like temperature behavior 43 . Therefore we can safely exclude the contribution of these three layers to the electrical transport and assume that, for capped BBO, the electrical transport flows trough two parallel channels related to BBO and Bi-deficient BBO layers (see the sketch of Figure 2(a)). This scenario was considered to fit the experimental data by proposing an equivalent circuit formed by a temperature independent resistance R BBO for the BBO layer in parallel with a polaronic transport channel, attributed to the Bi-deficient BBO layer, with a temperature dependent resistance R P ol (T) given by 44

R P ol (T ) = K 0 a g d e 2 ν 0 k B T exp E σ k B T , (1)
where k B is the Boltzmann constant, T is the absolute temperature, e is the electronic charge, a is the polaron hopping distance, ν 0 is the jump attempt frequencyrelated to the phonon vibrations coupled to the polaron-, E σ is the activation energy for polaron hopping, g d is an adimensional factor that depends on the hopping geometry [START_REF]d =N ca 3 , where N is the available sites for a jump, c is the defect concentration and a was defined in the text[END_REF] and K 0 =L/Wt, with L, W and t being the length, width and thickness of the Bi-deficient BBO layer, respectively.

The R(T) for the bilayer is thus given by

R(T ) = R BBO R P ol (T ) R BBO + R P ol (T ) . ( 2
)
The fitting of R(T) is shown with a red full line in Figure 2(c) and displays a reasonable agreement with the experimental data for E σ =0.12 eV and K 0 ak B g d e 2 ν0 = 5.0Ω/K. If we take ν 0 =1.710 13 Hz (the breathing mode vibration frequency 45 ) , a=0.437nm -corresponding to this perovskite cell parameter-, g d =3/8 59 and K 0 =0.067 nm -1 (L/W≈2 given by our sample geometry and t=30 nm, as extracted from the TEM analysis of Figure 1) we ob- tain a theoretical value K 0 ak B g d e 2 ν0 = 2.5Ω/K, in very good agreement with the fitted value, giving consistency to our analysis. We recall that it has been theoretically proposed that, upon hole doping, BBO remains insulating due to the coupling between holes and phonons 46 and the subsequent formation of polarons. This is likely the case of Bi-deficient BBO, where the holes for polaron formation are provided by Bi vacancies, as it will be analyzed in the next section by means of first principles calculations.

III. DFT CALCULATIONS

To study the different physical mechanisms that could trigger the transport properties of Bi-deficient BaBiO 3 , we perform ab initio calculations within the framework of Density Functional Theory (DFT) and the projector augmented wave (PAW) method 47 , as implemented in the Vienna ab initio package (VASP) 48,49 . We explicitly treat 10 valence electrons for Ba (5s 2 5p 6 6s 2 ), 15 for Bi (5d 10 6s 2 6p 3 ) and 6 for O (2s 2 2p 4 ). To properly describe the electronic structure of charge ordered and semiconducting BaBiO 3 -specially in the presence of localized defects-, it is necessary to consider long range exchange interactions that are poorly described by the standard local and semilocal functionals as the local density approximation (LDA) or the general gradient approximation (GGA). Therefore, we use the semilocal functionals developed by Perdew, Burke and Ernzerhof (PBE) 50 and we add the hybrid functional including a fraction of non-local Hartree-Fock exchange interactions through the Heyd-Scuseria-Ernzerhof (HSE06) functionals 51,52 in all the calculations. The charge-ordered phase of BaBiO 3 is represented formally by the Ba +2 2 Bi +3 Bi +5 O -2 6 formula unit. The two Bi species are ordered in an alternated way within a monoclinic structure 53 . To compute the electronic properties of perfect BaBiO 3 and Bi-deficient BaBi 0.75 O 3 , containing a 25% of Bi-vacancies, we use a supercell rotated 45 • with respect to the simple cubic single perovskite and which contains 4 Bi sites in total, 2 per layer and two planes in the z-direction. To explore lower concentrations of Bi-vacancies we also perform calculations in a 2x2x2 FCC supercell containing 80 atoms 60 , as the one used in Ref. 46. This larger supercell contains 16 Bi-atoms per unit cell and we account for defect concentrations of 12.5% (2 Bi-vacancies) and 18.75% (3 Bi-vacancies). The HSE calculations within the smaller unit cell are performed using a 450 eV energy cutoff in the plane waves basis and to evaluate the integrals within the Brillouin Zone (BZ), a 4×4×3 Monkhorst-Pack k-point grid is employed. The corresponding values for the energy cutoff and k-points grid in the larger FCC supercell are 400 eV and 2×2×2, respectively. The experimental lattice parameters 53 are used in all the calculations and internal structural relaxations are performed until the forces on each ion are less than 0.02 eV/ Å.

For perfect BaBiO 3 , our calculations predict the semiconducting state with a band gap of 0.90 eV, exhibiting the charge disproportionation between the two different species, Bi +3 and Bi +5 in agreement with previous calculations 54 . In Figure 3(a) we show the corresponding bandstructure where the full 6s bands of Bi +3 are depicted in red and the empty ones of Bi +5 are in dark green. It is important to mention that the Bi-O bonds can be described by spatially extended hybridized 6s-2p orbitals and all the bands present a strong s-p hybridization. The color coding used in Figure 3 and in the following ones is just a guide to the eye. Distinct kinds of defects appear in the presence of Bi vacancies. It is worth to mention that several molecular dynamics steps, followed by relaxation, were necessary to stabilize these defects in order to avoid spurious metallic solutions 46 . In the following, we describe the results obtained for the different Bi vacancies concentrations considered. The system with 25 % of Bi vacancy concentration has been simulated in the small cell where one out of four Bi is withdrawn. Independently of the site of the vacancy, whether it is at a Bi +5 or Bi +3 site, the resulting electronic structure is very similar. For instance, if the vacancy is at a Bi +5 site, in order to compensate the lacking five electrons, each of the remaining two Bi +3 liberates two electrons, so that the originally full 6s bands become empty. These two trapped holes at each Bi +3 form a bipolaron. Since the remaining Bi +5 site has its 6s band already empty, the other missing electron turns into a hole trapped at the surrounding oxygen cloud next to the Bi vacancy. On the other way around, if the vacancy is at a Bi +3 site, the lacking three electrons give rise to a bipolaron at the remaining Bi +3 site and a hole trapped at the oxygen p-bands surrounding the vacancy. In As mentioned before, the system with 12.5 % and 18.5 % Bi vacancy concentrations has been simulated with a supercell containing 16 Bi sites. In the case with an even amount of vacancies, the kind of electronic defects are only of the bipolaron type.

In Figure 3 (c), we show the calculated electronic bandstructure obtained for the 12.5 % concentration, where there are two vacancies at Bi +3 sites in the unit cell. The three occupied bands in red are the full 6s bands of the remaining Bi +3 , the three bands in blue are those of the bipolarons, and, at higher energies, the eight 6s bands of the Bi +5 are in dark-green. In the case of even number of vacancies, the system is non-spin polarized so that the bands are doubly degenerate. When the Bi vacancy concentration is 18.75 %, there are three vacancies in the supercell. For example, if two of the vacancies are at a Bi +3 site and the other one at a Bi +5 , there are eleven electrons missing. The six remaining Bi +3 sites can provide this charge. That is, five of these Bi +3 sites turn into Bi +5 by trapping two holes each giving rise to five bipolarons (blue bands in Figure 3 (d)). The remaining Bi +3 keeps one electron forming a polaron with a localized spin-polarized charge (violet band in Figure 3 (d)).

In order to analyse the role played by each type of defects into the transport properties we proceed to estimate the activated diffusion barrier E σ in each case. The hole polaron trapped at the p-bands of the oxygens that surround the Bi vacancy in the highly doped case (25% vacancy concentration), is a mixed type of defect where the Bi vacancy is bound to a hole electronic cloud. Nudged elastic bands (NEB) 55 method has been used to estimate the activation barrier, taking into account two different migration pathways. The first one considers the migration of the defect within the ab-plane, that is in the plane that contains the part of the hole shared by the 4 in-plane oxygen atoms. The second considered path was with the defect travelling along the (111) -( 111 ) direction, with the inflection point located at an optimized height. That is, the defect travelling from the (1/2, 1/2, 1/2) to the (0, 0, 0) position, through (1/4, 1/4, h 0 ), being h 0 the optimized distance in the z-direction from the ab-plane. We obtain E σ =9.35 eV and 4.9 eV, respectively. In view of these results, we conclude that the contribution to the electronic conductivity of these type of mixed defects is negligible.

To estimate the maximun of the energy barrier of the polaron and bipolaron trapped at Bi +3 sites (violet and blue bands of Figure 3, respectively), we perform a static calculation where the lattice distortion coupled to the defect is transferred to an intermediate situation that forces the trapped charge to be shared between neighboring Bi sites. For the case with 18.75 % Bi vacancy concentration, in the top panel of Figure 6 (a), we show the partial density of states (PDOS) projected onto the s-orbitals of two Bi neighboring sites, one of them trapping the polaron (labeled as Bi 12 ) and the other one belongs to a Bi +5 , named Bi 5 . The corresponding charge density associated with the polaron is shown in Figure 6 (b). In the bottom panel of Figure 6 (a), we plot the PDOS projected on the same Bi sites but in the situation where the polaron charge is shared in equal parts between them. In Figure 6 (c) we show the related charge density. The energy difference between the initial and the intermediate situation provides an estimation of E σ , that is 0.54 eV. When this approach is done to estimate E σ for the bipolaron migration, where now two holes are hopping simultaneously, we obtain E σ =0.60 eV. Interestingly, it is not surprising that this energy scale is related to the electron-phonon coupling in BBO systems, specially with the breathing phonon modes that are strongly correlated to Bi +5 -Bi +3 charge ordering. The value around 0.5 eV for E σ is of the same order of magnitude as the energy shift induced in the electronic bands when frozen phonon pertubations for the breathing mode are performed for bulk BBO 56 . However, as it has been recently reported in Ref. 21, this value decreases to ∼ 0.2 eV in confined situations as the one described for the BBO/BaPbO 3 interface. So, it is likely that the value of ∼ 0.5 eV obtained for the bulk supercell described here is overestimated since it does not take into account the confinement effects present in the real samples. Therefore, we consider that the actual activation energy might be close to the 0.12 eV obtained from the transport experiments and that the migration of polarons trapped at Bi +3 sites plays an important role in the increment of the electrical conductivity with the temperature, experimentally observed in our Bi deficient films. The existence of polarons or bipolarons trapped at Bi sites in BBO was already theoretically described in Ref. 46. The difference is that in that work the doping is caused by K substitution at Ba sites. They show that local lattice relaxations are sufficiently large to screen the localized holes and that hole doped BBO remains semiconducting upon moderate hole doping. In our work, we show that the level of doping with Bi vacancies is 3 to 5 times larger than with K. Notwithstanding this, it is remarkable that the system remains semiconducting for all the concentrations studied.

IV. CONCLUSIONS

We have shown that it is possible to synthesize thin layers of BBO perovskite with a large amount of Bivacancies. Transport properties were modelled with a polaronic model and the existence of these defects was theoretically validated. Our DFT calculations performed for Bi-deficient BBO supercells show that Bi vacancies give rise to different type of polarons, depending on the vacancies concentration. After the Bi vacancy is formed, the missing electrons induce polaron and/or bipolarons that consist in one or two holes trapped at some of the remaining Bi sites, or even at the oxygen octahedron surrounding the vacancy for large vacancies concentrations. We conclude that hole doped BBO remains semiconducting despite the large hole concentration induced by the Bi deficiency. We finally estimate the activated diffusion barrier for all these type of defects to analyse their contribution to the transport properties. We show that the calculated energy barriers of the charges trapped at Bi sites are correlated to the electron-phonon coupling energy scale. The calculated activated barriers for the simulated bulk supercells are somewhat larger than the one obtained from the fit of the transport experiments. In view of previous theoretical results for other BBO interfaces, these larger values are expected to decrease in confined situations -such as the Bi-deficient BBO layer studied in this work-and match the experimentally extracted activation energies. We finally conclude that the polarons trapped at Bi sites strongly contribute to the electrical conductivity of Bi-deficient BBO and, in a more general framework, that the control of Bi vacancies can be used as a tool for tuning the electronic properties of this perovskite.

FIG. 1 :

 1 FIG. 1: (a) X-ray diffraction spectra corresponding to both bare and capped BBO. Peaks indexes correspond to the pseudocubic (PC) perovskite BBO lattice. Peaks marked with * correspond to the sample holder used for the measurement of the capped BBO sample 35 . Zero-shift corrections were made from the position of Si peaks; (b) Blow-up of the spectra displayed in (a), evidencing the shift of diffraction peaks to lower angles for capped BBO. Maxima are shifted by 0.24 º. Also, capped BBO displayed broader diffraction peaks (FWHM of ≈ 0.5 º and ≈ 0.7 º for BBO and capped BBO, respectively) with a slight asymmetry; (c) TEM HAADF image corrresponding to a BBO/YSZ bilayer on Si/SiO2, with the YSZ grown at an oxygen pressure of 0.1mbar. The Pt layer on top of the YSZ layer was deposited for the lamella preparation.

FIG. 2 :

 2 FIG. 2: Sketches showing (a) the Si2/BBO/Bi-def. BBO/YSZ stack and (b) a top view of the contacted device. Small parts of SiO2 (top layer of the substrate, blue) were left exposed in order to access the BBO and Bi-deficient BBO layers for 4-points transport measurements. Same colours are used for (a) and (b); (c) 4-points resistance as a function of the temperature recorded for a bare BBO film and a YSZcapped heterostructure.The red line shows a fitting assuming a polaronic conduction model for the Bi-deficient BBO layer (see text for details).

FIG. 3 :

 3 FIG. 3: Spin-polarized band structure of (a) perfect BaBiO3 and Bi-deficient BaBiO3 with (b) 25 %, (c) 12.5 % and (d) 18.75% of Bi-vacancies. Solid lines correspond to spin-up projection and dashed lines to spin-down. In (a) and (c) spin-up and spin-down bands are degenerate. Red bands correspond to Bi +3 , dark-green ones to original Bi +5 atoms, and blue to created bipolarons. The orange band in (b) represents the hole trapped at the oxygen cloud centered at the Bi-vacancy sketched in Figure 4 (a) and the violet bands in (d) correspond to the charge/hole polaron originated around Bi12 (see text and Figure 6 (b))

  Figure 3 (b) we show the electronic bandstructure, with the orange narrow band corresponding to holes trapped at the oxygen cloud around the vacancy, whose charge density is depicted in Figure 4(a). The bands in blue are the bipolarons with the corresponding charge density shown in Figure 4(b).

FIG. 4 :

 4 FIG. 4: Charge density plots of Bi-deficient BaBiO3 with a 25% of Bi-vacancies. Panel (a) shows the hole trapped at the oxygen cloud around the vacancy, correspoding to the orange band depicted in Figure 3 b). The charge density shown in (b) corresponds to the bipolaron at the remaining Bi +3 site (blue band in Figure 3 b).

FIG. 6 :

 6 FIG. 6: (a-top panel) PDOS projected onto the s-orbitals of two selected Bi atoms, the one with the polaron, i.e. Bi12 and one of its nearest Bi neighbors, i.e. Bi5; (a-bottom) PDOS projected onto the same Bi atoms for the intermediate situation where the polaronic charge is half-transfered between neighboring sites; (b) Charge-density plot corresponding to the localized state close to EF shown in the top panel of (a), i.e. the polaron centered at Bi12; (c) Idem (b) after the static perturbation where the polaronic charge is shared with Bi5.

V. ACKNOWLEDGMENTS

The authors received financial support from PICTs 2015-0869, 2016-0867, 2017-1836 and 2019-02128 of ANPCyT, Argentina. Also, funding from EU-H2020-RISE project "MELON" (SEP-2106565560) is acknowledged. DR thanks the Nederlandse Organisatie voor Wetenschappelijk Onderzoek for a visiting grant, project 040.11.735, "Memcapacitive elements for cognitive devices". We thank D. Vega, from CAC-CNEA, for some of the XRD measurements.