Maria Tsekhmistrenko
email: maria.tsekhmistrenko@gmail.com

Karin Sigloch

Kasra Hosseini

Keywords: Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information

 could also point to more rising heat than admitted by the surface tallies based on narrow hotspot tracks 2,3 .

 , tomographic illumination across all mantle depths is particularly useful in providing a set of self-consistent geometric constraints that mantle modelling should reproduce. This is the aim of our study.

The Réunion Hotspot and Upper Mantle-Réunions Unterer Mantel experiment

). In first using RHUM-RUM body-wave data for tomography, our hemispheric-scale study is able to trace at least one upwelling from the African LLVP to its surface location in high detail.

 .

P-wave tomography from surface to core

 and detailed in the Methods. Seismic velocity anomalies are caused primarily by A tree of Indo-African mantle plumes imaged by seismic tomography

mantle temperature variations, with slower-than-average regions probably hot and therefore buoyant [START_REF] Masters | Earth's Deep Interior: Mineral Physics and Tomography from the Atomic to the Global Scale[END_REF] .

The teleseismic model in Fig. 2a,b shows a voluminous body of slow P velocities resides in the lower mantle beneath southern Africa-the African LLVP with its sharp western edge, tilting eastward towards the surface [START_REF] Ritsema | Complex shear wave velocity structure imaged beneath Africa and Iceland[END_REF][START_REF] Ni | Sharp sides to the African superplume[END_REF] . The upper mantle is slow beneath the Réunion hotspot and Madagascar. The three-dimensional (3D) hulls in Fig. 2b emphasize the smoothness of this conventional global tomography, characteristically lacking upper-mantle detail beneath the hotspot because seismic stations in the region were limited to a few widely spaced island locations.

The base of the mantle is constrained only marginally by existing P-wave tomographies [START_REF] Montelli | Finite-frequency tomography reveals a variety of plumes in the mantle[END_REF][START_REF] Grand | Global seismic tomography: a snapshot of convection in the Earth[END_REF] , but it is sampled extensively by our P diff waveform measurements. Figure 2c,d shows the effect of adding these data. The lower mantle becomes slower and more structured than in Fig. 2a,b. On the CMB, a narrow, slow corridor extends southeastwards from the LLVP's main doming, towards the viewer in Fig. 2c, rising 500-900 km above the CMB. We term this the Southern Indian Ocean (SIO) corridor.

The model resulting from the joint inversion of all three datasets, including RHUM-RUM, is named RROx-19 (Fig. 2e,f). Slow-velocity anomalies at all depths are further intensified and concentrated spatially, compared not just with Fig. 2a-d but also with published models (Extended Data Fig. 1). Crucially, the eastward-tilting mid-mantle structure sharpens into an alignment of distinctive blobs, each 700-800 km in diameter, which are robust in numerous resolution tests (Methods). Beneath the area instrumented by RHUM-RUM, the upper mantle acquires the expected ). Yellow circles and cyan hexagons represent broadband stations of partner projects SELASOMA [START_REF] Tilmann | SELASOMA Project, Madagascar 2012-2014[END_REF] and MACOMO [START_REF] Wysession | Investigation of sources of intraplate volcanism using PASSCAL broadband instruments in Madagascar, the Comores, and Mozambique[END_REF] , respectively. Lower panel shows different 3D views of usable teleseismic rays incident on OBS stations during their 13-month deployment. Red/blue rays represent seismically slower-/faster-than-average wave paths (modelled by finite-frequency kernels, not rays). Overall station event distribution can be viewed in Supplementary Fig. 1. R, La Réunion Island; M, Mauritius; SWIR, Southwest Indian Ridge.

layering of fast lithosphere underlain by slow asthenosphere and reveals small-scale, lateral subdivisions. The upwellings have a strongly 3D character. Figure 3 is an annotated zoom into the preferred joint inversion of Fig. 2f, additionally showing a northeast-southwest cross section through Afar and Bouvet hotspots and the East African Rift region. Beneath Réunion, a 200-km thin vertical conduit appears in the upper mantle (right dashed arrow in Fig. 2f) and is robustly resolved (Extended Data Fig. 2). The conduit is anchored in a slow blob (marked by a bullseye), the distal part of the 700-to 800-km-thick mid-branch that tilts up eastwards from the LLVP cusp. If a vertical slow conduit beneath Réunion continued below 1,000 km, it would be much thinner than this tilting structure because resolution in this entire vertical column is excellent according to the resolution tests (Supplementary Figs. 4 and5). Despite its disrupted downward connection, the hotspot sits vertically above slow lowermost mantle, on the SIO corridor's northern flank (Extended Data Fig. 3).

Another vertical, upper-mantle conduit sits beneath southernmost Madagascar (Fig. 3; ~100 km in diameter, left dashed arrow) and is well resolved (Extended Data Fig. 2). This conduit could relate to the unexplained, late Cenozoic onset of volcanism on Madagascar 27 , present-day low velocities in the Madagascan asthenosphere [START_REF] Pratt | Shear velocity structure of the crust and upper mantle of Madagascar derived from surface wave tomography[END_REF] and 1-2 km of dynamically supported uplift since 15 Ma 29 . However, no time-progressive volcanic track is known, and the geochemistry of central Madagascan volcanism requires no substantial excess temperatures [START_REF] Cucciniello | From olivine nephelinite, basanite and basalt to peralkaline trachyphonolite and comendite in the Ankaratra volcanic complex, Madagascar: 40 Ar/ 39 Ar ages, phase compositions and bulk-rock geochemical and isotopic evolution[END_REF] . Farther west, LLVP-connected anomalies extend to asthenospheric depths of ~250 km but appear contained beneath old Jurassic lithosphere of the Mozambique Channel [START_REF] Mazzullo | Anisotropic tomography around La Réunion Island from Rayleigh waves[END_REF] .

Under the Central Indian Ridge, the uppermost mantle is slow (Fig. 2e), as expected for active spreading [START_REF] Mazzullo | Anisotropic tomography around La Réunion Island from Rayleigh waves[END_REF][START_REF] Barruol | Large-scale flow of Indian Ocean asthenosphere driven by Réunion plume[END_REF] . The westward-down extension towards the mid-mantle blob is probably a smearing artefact due to lack of instrumentation because surface-wave tomography confirms slow velocities beneath the ridge to be limited to the upper 100 km [START_REF] Mazzullo | Anisotropic tomography around La Réunion Island from Rayleigh waves[END_REF][START_REF] Barruol | Large-scale flow of Indian Ocean asthenosphere driven by Réunion plume[END_REF] .

Structural elements of the LLVP

The 3D connectivity of slow structure (presumably upwellings) is rendered in Fig. 4, centred on southern Africa. Figure 4a renders all low-velocity anomalies between 1,000 and 2,900 km (colour codes for depth); Fig. 4b-e separates the same structure into narrower depth bins.

The lowermost mantle (2,900-2,400 km, Fig. 4b) features a single, hemisphere-spanning structure, the African LLVP. It consists of a broadened central part under southern Africa, where three narrowly elongated corridors converge. The longest is the SIO corridor (~7000 km long, ~800 km wide on the CMB). The Kerguelen hotspot sits above its distal end, and Réunion and Marion hotspots near its northern and southern flanks, respectively. The East African (EA) corridor extends northwards under Kenya, atop which the Afar hotspot was sitting ~30 Ma [START_REF] O'neill | On the uncertainties in hot spot reconstructions and the significance of moving hot spot reference frames[END_REF][START_REF] Matthews | Global plate boundary evolution and kinematics since the late Paleozoic[END_REF] when its flood basalts erupted [START_REF] Ring | The east African rift system[END_REF] . (Africa has drifted north since.) A third, shorter corridor is incident from the west (South Atlantic). These CMB corridors are discernible in some previous global tomographies (Extended Data Fig. 4) but are rendered more narrowly elongated by our P diff data, hence more distinct from the LLVP's central region under South Africa.

Between 2,300 and 1,900 km (Fig. 4c), slow velocities fade above the corridors and gather in a narrow dome under South Africa, forming the LLVP's cusp (Fig. 3, dashed purple line). A smaller dome is present beneath the Kerguelen area. At 1,600-1,300 km depth (Fig. 4d), the Kerguelen dome has faded, and the central dome has split into three diverging branches, one running southeast (above the SIO CMB corridor), another running north (above the EA corridor), and the third running southwest towards the Bouvet hotspot.

The three branches continue to diverge and rise between 1,100 and 900 km depth (Fig. 4e). South of Madagascar, the SIO splits into two sub-branches running towards Réunion and Kerguelen.

From deep to shallow levels, slow anomalies become increasingly small scale and scattered. In composite Fig. 4a, slow anomalies are notably absent above the LLVP's central cusp, instead clustering above the red-brown SIO and EA corridors. Into the upper mantle, we lose resolution outside instrumented areas, but the general increase in structural complexity is real (further confirmed by Supplementary Fig. 4a,b, above 800 km) and presumably reflects the lower viscosities of the upper mantle [START_REF] Hager | Lower mantle heterogeneity, dynamic topography and the geoid[END_REF] .

The key observation of Fig. 4a, of the two major tilting branches (EA, SIO) vertically overlying the two major LLVP corridors, means that the essence of 3D structure can be captured by the 2D sections of Fig. 3 and cartoon Fig. 4f. Figure 3b slices through the EA and Bouvet mid-mantle branches, the former rising from the LLVP cusp northwards, toward Malawi and Kenya, further connecting to Afar through anomalously slow upper mantle.

The mid-mantle branch towards Bouvet rises southwestwards from the cusp, similar to the SIO and EA branches in overall appearance, though less well resolved (Supplementary Fig. 2). A slow CMB corridor underlying the Bouvet branch is currently not resolvable (see targeted test, Extended Data Fig. 5) but may be present because most global models in Extended Data Fig. 4 show the resolvable vestiges expected for the case that a corridor were present.

The tilting East African Rift branch was well known [START_REF] Ritsema | Complex shear wave velocity structure imaged beneath Africa and Iceland[END_REF][START_REF] Hansen | Mantle structure beneath Africa and Arabia from adaptively parameterized P-wave tomography: implications for the origin of Cenozoic Afro-Arabian tectonism[END_REF] , and the tilting SIO branch is discernible in global tomographies (Extended Data Fig. 1). The Bouvet branch is newly resolved, as is the emergence of all three branches from one specific area and depth under the southern Mozambique Channel (the cusp), which thus becomes recognizable as the upper limit of the LLVP.

A family tree of Indo-African plumes

None of the existing models for LLVP-plume interaction [START_REF] Garnero | Continent-sized anomalous zones with low seismic velocity at the base of Earth's mantle[END_REF] anticipated or explain the observed, multiple occurrences of tilting mid-mantle branches, vertically underlain by slow CMB corridors. Three instances of this structural element converge under southern Africa, giving the LLVP a tree-like appearance (Fig. 4f). Narrow, root-like CMB corridors converge on a central doming region (the trunk). Three branches split off at the top of this trunk (~1,500 km depth), diverging upwards and further splitting into a canopy of splotchy, small-scale structure above 1,000 km.

Figure 5 proposes a new model for how the African LLVP generates plumes and hotspot volcanism. It accounts for our seismic observations and hindcasts the spatiotemporal patterns of plume-head eruptions in detail (large igneous provinces, LIPs). Each mid-mantle branch consists of several blobs or proto-plumes (~800 km in diameter, Fig. 3a), which rose vertically from the CMB corridor underneath. The blobs are diagonally aligned because they detached in temporally staggered sequence (not because plume material flows diagonally up within the branch).

Blob 1 detaches from the distal end of the CMB corridor (Fig. 5a). This promontory region is characterized by strong lateral temperature gradients (dT/dx) to its non-LLVP surroundings in three directions (to the right, and in/out of the viewing plane). Settings of strong dT/dx are known to be preferential formation regions for thermal instabilities [START_REF] Turner | Buoyancy Effects in Fluids[END_REF][START_REF] Namiki | The influence of boundary heterogeneity in experimental models of mantle convection with internal heat sources[END_REF] . The upward release of Blob 1 shifts the location of maximum dT/dx left along the corridor (Fig. 5b), where the next instability is generated (Blob 2). Departure of Blob 2 causes instability 3 to its left, and so on. Thus, a self-triggering cascade of blob-shaped detachments propagates along a hot CMB corridor, towards the centre of the LLVP. Upon reaching the boundary to the upper mantle with its drop in ambient viscosity, each proto-plume blob spawns a plume-head/plume-tail sequence (Fig. 5d). This explains why the thin plume-tail conduit under present-day Réunion is anchored in a mid-mantle blob rather than the CMB.

Tank and numerical models show a strong tendency for plume material to rise vertically [START_REF] Davaille | Simultaneous generation of hotspots and superwells by convection in a heterogeneous planetary mantle[END_REF][START_REF] Kumagai | Mantle plumes: thin, fat, successful, or failing? Constraints to explain hot spot volcanism through time and space[END_REF][START_REF] Steinberger | Conduit diameter and buoyant rising speed of mantle plumes: implications for the motion of hot spots and shape of plume conduits[END_REF] ; notable tilts are produced only when plume detachment regions sweep laterally across the CMB over time [START_REF] Steinberger | Conduit diameter and buoyant rising speed of mantle plumes: implications for the motion of hot spots and shape of plume conduits[END_REF] , as in Fig. 5. Moreover, diagonal upwards flow in the branches would not require or explain CMB corridors vertically underneath. Compared with earlier S-wave models, our P diff observations narrow the CMB corridors to <1,000 km, no wider than the detaching blobs-a line-like LLVP geometry removes or supersedes the debate of whether plumes rise from the centre versus the interior of the LLVP [START_REF] Burke | Plume generation zones at the margins of large low shear velocity provinces on the core-mantle boundary[END_REF][START_REF] Austermann | A statistical analysis of the correlation between large igneous provinces and lower mantle seismic structure[END_REF] .

Our model hindcasts a linearly staggered sequence of plume-head arrivals at the surface (LIP eruptions), which are indeed observed for the EA and SIO branches. Located distally above the EA corridor, the Afar region experienced flood basalt eruptions ~30 Ma, followed by transition from continental rifting to seafloor spreading. Continental rifting is increasingly recent and less advanced towards the southern EA rift [START_REF] Ring | The east African rift system[END_REF] , and the associated magmatism does not (yet) have LIP character, although slow-velocity anomalies are already flooding the upper mantle (Extended Data Fig. 6). As predicted, no rifting is observed south of Malawi/Mozambique, where the tilting Afar branch is still at mid-mantle depths.

Similarly, the Kerguelen flood basalts occurred ~120 Ma 19,41,42 above the distal SIO corridor and moved inwards to Marion ~90 Ma [START_REF] Morgan | The Oceanic Lithosphere[END_REF][START_REF] Gibbons | The breakup of East Gondwana: assimilating constraints from Cretaceous ocean basins around India into a best-fit tectonic model[END_REF][START_REF] Coffin | Kerguelen hotspot magma output since 130[END_REF] (Madagascar Rise LIP) and to Réunion ~65 Ma (Deccan Traps LIP) [START_REF] Courtillot | Three distinct types of hotspots in the Earth's mantle[END_REF][START_REF] Morgan | The Oceanic Lithosphere[END_REF][START_REF] Gibbons | The breakup of East Gondwana: assimilating constraints from Cretaceous ocean basins around India into a best-fit tectonic model[END_REF] , but no flood basalts have occurred closer to the cusp, where the upwelling SIO branch is still at mid-mantle depths. The three LIP arrivals above the SIO corridor are thought to have created the Indian Ocean as they unzipped India from its Gondwana supercontinent surroundings, in the sequence India from Antarctica [START_REF] Morgan | The Oceanic Lithosphere[END_REF][START_REF] Gibbons | The breakup of East Gondwana: assimilating constraints from Cretaceous ocean basins around India into a best-fit tectonic model[END_REF][START_REF] Coffin | Kerguelen hotspot magma output since 130[END_REF] ~120 Ma (Kerguelen); India from Madagascar [START_REF] Morgan | The Oceanic Lithosphere[END_REF][START_REF] Gibbons | The breakup of East Gondwana: assimilating constraints from Cretaceous ocean basins around India into a best-fit tectonic model[END_REF][START_REF] Storey | Timing of hot spot-related volcanism and the breakup of Madagascar and India[END_REF] ~90 Ma (Marion); and India from Seychelles [START_REF] Morgan | The Oceanic Lithosphere[END_REF][START_REF] Gibbons | The breakup of East Gondwana: assimilating constraints from Cretaceous ocean basins around India into a best-fit tectonic model[END_REF][START_REF] Coffin | Kerguelen hotspot magma output since 130[END_REF] ~65 Ma (La Réunion). Similarly, the plume heat emanating from the EA corridor may be enabling the splitting of Africa, potentially creating a new ocean basin [START_REF] Ring | The east African rift system[END_REF] . Thus, SIO and EA corridor geometries foreshadowed the strikes of major future spreading systems. The linear staggering of proto-plume detachments (Fig. 5) explains the effectiveness of flood basalts in unzipping continents and creating new spreading ridges [START_REF] Morgan | The Oceanic Lithosphere[END_REF][START_REF] Gibbons | The breakup of East Gondwana: assimilating constraints from Cretaceous ocean basins around India into a best-fit tectonic model[END_REF][START_REF] Nürnberg | The tectonic evolution of the South Atlantic from Late Jurassic to present[END_REF] (compared with randomly sequenced LIP arrivals).

This newly predictable sequencing of LIPs informs the question of whether flood basalts precede continental rifting or vice versa [START_REF] Courtillot | On causal links between flood basalts and continental breakup[END_REF][START_REF] White | Magmatism at rift zones: the generation of volcanic continental margins and flood basalts[END_REF] . Under Kenya, for example, a future LIP seems to be rising, but rifting and hot asthenosphere inflow may have long been encroaching from the north (Afar). Thus, Kenyan rifting precedes its LIP but is driven (or facilitated) by the adjacent, preceding LIP within the EA corridor sequence.

We predict cataclysmic flood basalts, some tens of millions of years in the future, for the surface overlying the LLVP cusp, where three inward-sweeping mid-mantle branches converge-possibly analogous to widespread, Cretaceous plateau eruptions above the Pacific LLVP [START_REF] Burke | Plume generation zones at the margins of large low shear velocity provinces on the core-mantle boundary[END_REF][START_REF] Suetsugu | South Pacific mantle plumes imaged by seismic observation on islands and seafloor[END_REF][START_REF] Morgan | The Oceanic Lithosphere[END_REF][START_REF] Courtillot | On the ages of flood basalt events[END_REF][START_REF] Madrigal | Record of massive upwellings from the Pacific large low shear velocity province[END_REF] . Upwelling of the cusp probably already shapes the large-scale topography anomaly of southern Africa, which stands 500-1,000 m higher than expected [START_REF] Jones | Geochemical and geophysical constrains on the dynamic topography of the Southern African Plateau[END_REF] .

Our model of blobby proto-plumes (Fig. 5) predicts a Réunion-like situation under every active hotspot of the region: a classical, thin and vertical plume conduit, but only down to ~700 km depth. There it is anchored in a mid-mantle blob without a tail to the CMB (for example, ref. [START_REF] Davaille | Convective patterns under the Indo-Atlantic «box[END_REF] ). Currently resolved structure under Afar [START_REF] Ritsema | Complex shear wave velocity structure imaged beneath Africa and Iceland[END_REF][START_REF] Hansen | Mantle structure beneath Africa and Arabia from adaptively parameterized P-wave tomography: implications for the origin of Cenozoic Afro-Arabian tectonism[END_REF] , Bouvet (Extended Data Fig. 6) and Kerguelen is consistent with this prediction. Hence, the observed, near-stationarity of Indo-Atlantic hotspots (relative motions of <10 mm yr -1 ), long harnessed for constructing absolute reference frames [START_REF] Morgan | The Oceanic Lithosphere[END_REF][START_REF] O'neill | On the uncertainties in hot spot reconstructions and the significance of moving hot spot reference frames[END_REF] , may measure the advection of mid-mantle reservoirs at 700-1,500 km depth rather than in the lowermost mantle.

The thin, upper-mantle conduits under Réunion and Madagascar may represent mostly thermal plumes released from the proto-plume, whose denser components tend to stay behind in the broad mid-mantle blob (Fig. 5d). Tailless blobs of large diameters (700-800 km under Réunion) point towards a dense chemical component in the LLVP and in the proto-plumes [START_REF] Kumagai | Mantle plumes: thin, fat, successful, or failing? Constraints to explain hot spot volcanism through time and space[END_REF] . Dense blobs need to become hotter to rise, perhaps partly explaining their robust tomographic visibility. Purely thermal plumes require smaller excess temperatures to rise, should be thin and should retain their stems to the CMB [START_REF] Kumagai | Mantle plumes: thin, fat, successful, or failing? Constraints to explain hot spot volcanism through time and space[END_REF] .

The hotspots that form the branches of the southern African LLVP tree essentially coincide with the ensemble of hotspots defining the geochemical DUPAL anomaly [START_REF] Dupre | Pb-Sr isotope variation in Indian Ocean basalts and mixing phenomena[END_REF][START_REF] Jackson | Evidence for a deep mantle source for EM and HIMU domains from integrated geochemical and geophysical constraints[END_REF] , which distinguishes hotspot basalts of the Indo-Southern Ocean from others around the world. Hence, the source region of the DUPAL anomaly may be identified with the LLVP tree structure below its branching depth of ~1,500 km. Indian and Southern Ocean hotspot eruptions that were previously not known to be specifically related beyond their shared DUPAL signature are aligned into predictable spatiotemporal sequences by our observation of narrow LLVP corridors paired with vertically detached upwellings in the mid-and upper mantle. This inference has been enabled by the very detailed imaging, from lithosphere to core, of one complete 'unzipping corridor' under the On reaching the viscosity boundary between lower and upper mantle, proto-plume 1 spawns the 'classical' sequence of a short-lived plume head (LIP), followed by one or more thin, long-lived plume tails, which remain rooted in their mid-mantle blob. Arrival of the plume head dynamically uplifted the lithosphere [START_REF] Cande | Indian and African plate motions driven by the push force of the Réunion plume head[END_REF] (vertically exaggerated). Proto-plume 1 corresponds with the Réunion plume in Fig. 3 (top) and Fig. 4f; 3 corresponds with the LLVP cusp.

Articles

Nature GeoscieNce

Methods RHUM-RUM. The RHUM-RUM experiment instrumented 2,000 × 2,000 km of the western Indian Ocean (Fig. 1). Its central component was a 13-month deployment of 57 German and French broadband OBSs in 2,300-to 5,600-m-deep water, supplemented by 2-3 yearlong deployments of 37 island stations on La Réunion, Mauritius, Rodrigues, the southern Seychelles, the Iles Eparses and southern Madagascar. Two complementary partner projects across Madagascar [START_REF] Wysession | Investigation of sources of intraplate volcanism using PASSCAL broadband instruments in Madagascar, the Comores, and Mozambique[END_REF] and Mozambique 52 extended the deployment westwards with another 33 stations. Combined with a cluster of OBSs that stretched eastwards to the Central Indian Ridge, the overall array was elongated in an east-west direction, making optimal use of crossing wave paths from the two most prolific incidence directions of teleseismic earthquakes: the southwest Pacific subduction zones to the east and the South Sandwich subduction to the west. A line of distant southerly OBSs recorded seismicity from the Alpine-Himalayan mountain belts, and Seychelles stations recorded sparser events from the south.

Pre-processing. The RHUM-RUM dataset [START_REF] Barruol | HUM-RUM experiment, 2011-2015, code YV (Réunion Hotspot and Upper Mantle -Réunion's Unterer Mantel) funded by ANR, DFG[END_REF] has been assigned the FDSN network code YV and is freely accessible at the French RESIF data centre (http://seismology. resif.fr). Large volumes of seismic waveforms were retrieved and pre-processed using obspyDMT [START_REF] Hosseini | ObspyDMT: a Python toolbox for retrieving and processing large seismological data sets[END_REF] , a program that automatically applies instrument correction and filters and cuts the waveforms around possible earthquakes according to a selected earthquake catalogue. A general description and quality assessment of the waveforms recorded by OBS stations was presented by Stähler et al. [START_REF] Stähler | Performance report of the RHUM-RUM ocean bottom seismometer network around La Réunion, western Indian Ocean[END_REF] ; the main challenge in working with the ocean-bottom data was ambient (microseismic) noise and high self-noise (only in the German DEPAS stations). Dataset 3: teleseismic, multi-frequency measurements on Indian Ocean stations. We followed the approach of Sigloch and Nolet [START_REF] Sigloch | Measuring finite-frequency body-wave amplitudes and traveltimes[END_REF] , briefly outlined in the following. Observed broadband waveforms from RHUM-RUM OBS, island stations and other land stations were Gabor-filtered in eight frequency passbands between 30 and 2.7 s centre period. To measure travel-time delays, these observations were cross correlated with bandpassed synthetic waveforms [START_REF] Hosseini | Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves[END_REF][START_REF] Sigloch | Measuring finite-frequency body-wave amplitudes and traveltimes[END_REF] (not with neighbouring observed seismograms), yielding travel-time delays relative to spherically symmetric Earth model IASP91. A measurement was accepted for input to the tomography if the Pearson cross-correlation coefficient exceeded 0.8. This procedure yielded the 55,657 teleseismic P-wave travel times in the eight frequency passbands, from 25,130 unique source-receiver paths, generated by 635 earthquakes that occurred 2011-2016 (Supplementary Fig. 1). Multi-frequency measurements enable selectively rejecting data in noisy frequency bands, while often retaining at least one acceptable band per wave path. Since OBS seismograms are more affected by microseismic noise than are land stations, the method salvages a maximum of wave paths for these expensive and valuable measurements.

Synthetic waveforms are calculated as the convolution of a (broadband) Green's function with a (broadband) source time function (STF), which are subsequently Gabor-filtered to the same passbands as the observations. Fully numerical Green's functions are calculated in the broadband range of 0.2 MHz to 1 Hz, using the axisymmetric spectral element method AxiSEM [START_REF] Nissen-Meyer | AxiSEM: broadband 3-D seismic wavefields in axisymmetric media[END_REF] (via its look-up database Instaseis [START_REF] Van Driel | Instaseis: instant global seismograms based on a broadband waveform database[END_REF] ). We propagated waves through the spherically symmetric Earth model IASP91 [START_REF] Kennett | IASPEI 1991 seismological tables[END_REF] with density and attenuation taken from PREM [START_REF] Dziewonski | Preliminary reference Earth model[END_REF] .

For each of the 635 earthquakes, a broadband STF is estimated (jointly deconvolved) from a suitable, spatially evenly distributed subset of broadband seismograms recorded for the event in question and from their broadband synthetics. For STF deconvolution, we use all stations at teleseismic distances that are part of the permanent, international networks (network codes II, IU, G or GE), but no RHUM-RUM stations or other temporary, regionally clustered networks. The computation of STFs is required for making correlation measurements because our study is almost global in scale. We cannot calculate finite-frequency travel times by correlating data with data (as in regional studies) because the waveforms may look too different to correlate well (due to the effects of radiation pattern and depth phases pP and sP). We have to correlate data with synthetics, and our frequencies are so high, and most events so shallow (depth phases pP and sP spaced by only a few seconds), that approximating the STF as Dirac pulse does not yield satisfactory correlations. STF shape is also used as a quality check; that is, events are not used for subsequent travel-time measurements when their STF estimates look non-physical [START_REF] Hosseini | Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves[END_REF] (see Sigloch and Nolet 58 for details).

Dataset 2: multi-frequency travel-time measurements of P diff . These deepest-diving body waves are abundantly observed on broadband seismograms, but they have hardly been used in tomography [START_REF] Wysession | Large-scale structure at the core-mantle boundary from diffracted waves[END_REF][START_REF] Kárason | Tomographic imaging of the lowermost mantle with differential times of refracted and diffracted core phases (PKP, Pdiff)[END_REF] because computational modelling of their diffraction by Earth's core is challenging. In a decade-long method-development effort, we have measured multi-frequency cross-correlation travel times on broadband, core-diffracted P waveform [START_REF] Hosseini | Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves[END_REF][START_REF] Hosseini | Multifrequency measurements of core-diffracted P waves (Pdiff) for global waveform tomography[END_REF] , modelled by the full-waveform propagation method AxiSEM [START_REF] Nissen-Meyer | AxiSEM: broadband 3-D seismic wavefields in axisymmetric media[END_REF] . Figure 2c,d illustrates the benefits of jointly inverting P diff travel times with the teleseismic ISC-EHB picks. P diff waves dive so deep that they sense Earth's core and are diffracted around it. Their abundance and relatively short wavelengths make them the most highly resolving wave type that extensively samples the lower third of the mantle. These longest-distance P waves (epicentral distances of ~97° or greater) have been made usable for tomography by Hosseini and Sigloch [START_REF] Hosseini | Multifrequency measurements of core-diffracted P waves (Pdiff) for global waveform tomography[END_REF] and Hosseini et al. [START_REF] Hosseini | Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves[END_REF] . This study uses a regionally optimized subset of 189,439 cross-correlation measurements of P diff travel times in five frequency bands (periods 30 to 7.5 s), taken from the global dataset of ~1,400,000 travel times inverted by Hosseini et al. [START_REF] Hosseini | Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves[END_REF] . Full-waveform synthetics were calculated using the AxiSEM 59 software, as for dataset 3.

Dataset 1: teleseismic ISC-EHB travel-time picks. The backbone dataset of over 2 million analyst-picked, teleseismic P-wave arrivals, was chosen because these data have been recorded around the sparsely instrumented Indian Ocean for three decades longer than modern broadband networks. To maximize the data coverage in the mid-and upper mantle beneath the western Indian Ocean, and to embed our new multi-frequency measurements in a global context, we complemented the inversion by a subset of analyst-picked P-wave arrival times from the EHB 21 catalogue (picks reported between 1970 and 1999) and from the new ISC-EHB 22 catalogue (2000-2014). Our aim was to select the smallest useful subset of teleseismic P-wave picks that would fill as many gaps as possible in our multi-frequency data coverage. We considered P-wave arrival picks at epicentral distances 30°-90° if the precision of their arrival-time reading was better (smaller) than 0.1 s. We also used the method of Hosseini et al. [START_REF] Hosseini | Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves[END_REF] to homogenize the directional density coverage of measurements at source-receiver locations. The final selection includes 2,342,621 teleseismic P-wave phase picks.

Tomography. All P-wave travel times were corrected for crustal structure, topography and ellipticity. We correct for crustal structure using the model CRUST2.0 66 in a ray-theoretical framework [START_REF] Nolet | A Breviary of Seismic Tomography[END_REF] , that is, not accounting for frequency dependence [START_REF] Ritsema | Long-period body wave traveltimes through the crust: implication for crustal corrections and seismic tomography[END_REF] . Finite-frequency sensitivity kernels were computed for the P diff travel times (as in ref. [START_REF] Hosseini | Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves[END_REF] ) and for the teleseismic RHUM-RUM travel times (following the method of Dahlen et al. [START_REF] Dahlen | Frechet kernels for finite-frequency traveltimes-I[END_REF] ). The ISC data were modelled with ray-like sensitivities. Each sensitivity kernel is numerically integrated on the tetrahedral inversion grid [START_REF] Hosseini | Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves[END_REF][START_REF] Montelli | A catalogue of deep mantle plumes: new results from finite-frequency tomography[END_REF][START_REF] Sigloch | Mantle provinces under North America from multifrequency P wave tomography[END_REF] . The earthquakes were relocated during the inversion stage. A combination of norm damping and Laplacian smoothing operators was used to regularize the ill-posed linear system of equations. By changing the intensity of the two regularization parameters, the trade-off between residual data misfit ||Gm-d|| 2 (where G is the matrix containing the sensitivity kernels discretized over our model m, and d is the data vector containing the travel-time measurements) in terms of χ 2 red (absolute χ 2 divided by the number of data N) can be written as:

χ 2 red = 1 N N ∑ i=1 ( d obs i -d pred i σi ) 2
and the model norm ||m|| 2 can be analysed. In this equation, d i obs is the observed travel time, d i pred is the predicted travel time and σ is the a priori uncertainty assigned to each measurement according to the cross-correlation coefficient x c . The values of x c translated to assigned measurement uncertainties σ as follows: 0.8 ≤ x c < 0.85 to σ = 1.0 s; 0.85 ≤ x c < 0.90 to σ = 0.6 s; 0.90 ≤ x c < 0.95 to σ = 0.4 s; 0.95 ≤ x c < 1.0 to σ = 0.3 s. All ISC measurements are assigned an uncertainty of σ = 1.0 s.

Each L curve is derived by assuming a constant ratio of the damping smoothing parameters but with different absolute amplitudes. χ 2 red = 1 is the optimal value to target if all σ are perfectly known. In practice, imperfect knowledge about the measurement uncertainties means some experimentation and subjective model appraisal in a relatively narrow range around this theoretically optimal χ 2 value (usually models close to the bending of the L curve). Our preferred model has χ 2 = 1.02.

The resulting tomography models can be examined in Fig. 2. In the ISC-only model (Fig. 2a), the uppermost mantle between the Madagascar and the Réunion areas is laterally uniformly red (slow) and lacks any vertical layering. The shallow, fast (blue) blob just west of Réunion in Fig. 2c (versus Fig. 2a) is introduced by paths incident on the Réunion station 'RER' from core-diffracted distances from the west (various azimuths, Andean earthquakes and Mid-Atlantic Ridge). By contrast, western paths from teleseismic distances in Fig. 2a are essentially limited to the South Sandwich arc. These teleseismic observations are evidently not sufficiently numerous and/or not sufficiently localizable along their ray paths to 'override' the tomography's smoothing constraint, which acts to blur the red anomalies to either side of the blue anomaly (Fig. 2c) into a single, laterally extended red anomaly in Fig. 2a. All of these structures in the ISC-only and in the ISC + P diff models are essentially constrained by only one station on Réunion, one on Mauritius, and one or two on Madagascar. In Fig. 2e, the RHUM-RUM data confirm and intensify the intervening blue anomaly. In addition, the laterally extended RHUM-RUM network introduces a vertical layering of lithosphere, asthenosphere and transition zone in the upper mantle under this region.

The three tomography models compared in Fig. 2 are chosen such that they feature most nearly the same global χ 2 misfit statistics (of close to 1.0 for the global dataset and ~1.4 for data recorded in the Indian Ocean and Africa). The data fit for the ISC data does not deteriorate when the P diff data get introduced in addition (transition from Fig. 2a,b to Fig. 2c,d), and the fits of ISC and P diff datasets do not deteriorate when the RHUM-RUM data get introduced (transition from model Fig. 2c,d to Fig. 2e,f). This combined model fits all three datasets without conflicts or trade-offs.

Visualization. All 3D visualizations (Figs. 2, 3 and 4) were done using ParaView [START_REF] Ahrens | The Visualization Handbook[END_REF] (http://www.paraview.org), an open-source data analysis and visualization application. Low velocities are rendered at several isovalues (dV p /V p = -0.25%, -0.5%, -1%, -2% and -3%); this leads to a 3D rendered effect, where seismically slow structures are enclosed by the nested isosurfaces. In Figs. 2 and3, the corresponding, nested cylindrical colour bar clarifies this visualization concept. Each cylinder represents one velocity isovalue, several of which are stacked into each other. Cut obliquely, we can identify the different cylinders and their isovalues (darker red means slower dV/V). The same idea of nested hypersurfaces applies when assessing the tomographic model.

Figures that present nested isosurfaces in thick cross-sectional slices always render volumes that are ring shaped (segments of a great circle) and 1,000 km wide; that is, the rendered region extends 500 km towards the viewer and 500 km away from the viewer, with a conventional, 2D, red-blue cross section forming the dividing plane. In addition to this, Figs. [START_REF] Hosseini | Global mantle structure from multifrequency tomography using P, PP and P-diffracted waves[END_REF] , GAP-P4 [START_REF] Obayashi | Finite frequency whole mantle P wave tomography: improvement of subducted slab images[END_REF] , HMSL-P06 [START_REF] Houser | Shear and compressional velocity models of the mantle from cluster analysis of long-period waveforms[END_REF] , PRI-S05 [START_REF] Montelli | A catalogue of deep mantle plumes: new results from finite-frequency tomography[END_REF] , S363ANI + M 76 , S40RTS [START_REF] Ritsema | RTS: a degree-40 shear-velocity model for the mantle from new Rayleigh wave dispersion, teleseismic traveltime and normal-mode splitting function measurements[END_REF] , savani [START_REF] Auer | a variable resolution whole-mantle model of anisotropic shear velocity variations based on multiple data sets[END_REF] , SEMUCB-WM1 79 and TX2015 [START_REF] Lu | The effect of subducting slabs in global shear wave tomography[END_REF] . It is generally observed that dV/V anomalies in P-wave models have lower amplitudes than in S-wave models. This difference is driven by the relatively larger changes in shear modulus versus bulk modulus, (mainly) as a function of temperature heterogeneities in the mantle.

Extended Data Fig. 1 compares models in the 1,000-km-thick slice containing the Réunion mid-mantle branch and the SIO CMB corridor beneath it. DETOX-P3 and GAP-P4 capture the general structure of the LLVP and cusp area but lack resolution beneath Réunion. PRI-S05, S40RTS and TX2015 have the highest overall resemblance to RROx-19. All models capture the sharp, tilting western edge of the LLVP. None of the models shows comparable details to RROx-19 in the upper mantle, and most of them fail to image a connection of the LLVP with the upper mantle. Notably, DETOX-P3 does not resolve nearly as much structure under this study area, despite incorporating RHUM-RUM data in the inversion. This is no surprise since the data coverage in a global model is very heterogeneous and a dataset like RHUM-RUM would not fall into full focus.

Extended Data Fig. 6 compares models in the 1,000-km-thick slice containing the two mid-mantle branches towards Afar (East Africa) and Bouvet. The Afar branch has been well imaged in previous studies [START_REF] Ritsema | Complex shear wave velocity structure imaged beneath Africa and Iceland[END_REF][START_REF] Hansen | Mantle structure beneath Africa and Arabia from adaptively parameterized P-wave tomography: implications for the origin of Cenozoic Afro-Arabian tectonism[END_REF] . DETOX-P3 and GAP-P4 show the general shape of the LLVP and the connection to Afar. On the CMB vertically beneath Afar, a low velocity is visible (Supplementary Fig. 4, best seen in models S362ANI + M, SEMUCB-WM1 and TX2015), which does not seem connected to low-velocity anomalies vertically above it, in the upper third of the mantle. Only PRI-S05, S40RTS and TX2915 capture similar structures as RROx-19, including the mid-branch towards Bouvet.

Extended Data Fig. 4 compares models within ~500 km above the CMB. Amplitudes in P models are lower but are in good agreement with S models about the extent of the LLVP. The southeast-striking SIO corridor is prominent in all models but is more narrowly focused in RROx-19 than in other models (except perhaps S40RTS). All models (except for GAP-P4) show signs of a north-south corridor along the CMB extending from South Africa towards Afar. No model shows a clear, narrow CMB corridor towards Bouvet, but there are hints (a southwestward bulge of the slowest isosurfaces) in RROx-19, DETOX-P3 and the S-wave models (PRI-S05, SEMUCB-WM1, S40RTS and S362ANI + M).

Resolution tests. For resolution tests, synthetic travel-time data vectors are computed on the basis of synthetic, user-defined input models and a noise term (with realistic statistical properties). Solutions are then obtained through the same inversion procedure and parameter settings as for the observed data. The comparison between input and output models gives insights into how faithfully and completely the tomography can recover the hypothetical Earth structures.

In Extended Data Fig. 2, we assume two narrow (r = 75 km and 700 km deep) vertical plume conduits beneath the southernmost tip of Madagascar and southwest of Réunion. Both conduits are resolved, but the Madagascar conduit is better resolved (less smearing) than the Réunion conduit due to OBS station coverage.

In Extended Data Fig. 5, we assume ~700-km-broad and ~400-km-high LLVP corridors to investigate the resolvability of these corridors. EA and SIO are well resolved whereas the Bouvet corridor vanishes. This is due to the low coverage of P diff phases in this undersampled region. Supplementary Fig. 3 is whole-mantle resolution tests using a Gaussian chequerboard. The chequers are Gaussian spheres where dV p /V p peaks at their centre (±3%) and tapers exponentially to 1/e of this value at 500 km radius (for the near-surface spheres; sphere radius decreases proportionally with depth.) Recovery in the upper mantle is good where stations are present. We observe excellent recovery of the chequerboard pattern between 800 and 2,000 km depth. Moderate

Fig. 1 |

 1 Fig. 1 | The RHUM-RUM experiment instrumented the southwestern Indian ocean with 57 broadband oBSs and 35 broadband land seismometers. Red dots represent broadband OBSs (October 2012-November 2013). Pink diamonds represent broadband land seismometers (2011-2016). Yellow circles and cyan hexagons represent broadband stations of partner projects SELASOMA[START_REF] Tilmann | SELASOMA Project, Madagascar 2012-2014[END_REF] and MACOMO[START_REF] Wysession | Investigation of sources of intraplate volcanism using PASSCAL broadband instruments in Madagascar, the Comores, and Mozambique[END_REF] , respectively. Lower panel shows different 3D views of usable teleseismic rays incident on OBS stations during their 13-month deployment. Red/blue rays represent seismically slower-/faster-than-average wave paths (modelled by finite-frequency kernels, not rays). Overall station event distribution can be viewed in Supplementary Fig.1. R, La Réunion Island; M, Mauritius; SWIR, Southwest Indian Ridge.

Fig. 2 |

 2 Fig. 2 | east-west cross sections through La Réunion and South Africa, showing the contributions of three tomography datasets to resolving whole-mantle structure. a, Global inversion using ISC-EHB/EHB teleseismic travel-time picks. The location of this 2D section is marked by the black great-circle line on the 3D globe in the top left. Red and blue shades indicate slower-and faster-than-average P velocities (V p ), respectively. Vertical black column extends from La Réunion down to the CMB. Md, Madagascar; CIR, Central Indian Ridge; PGZ, plume-generation zone. b, The same model with slow anomalies rendered as five velocity isosurfaces (dV p /V p = -0.25%, -0.5%, -1%, -2% and -3%) within a 1,000-km-thick, cross-sectional slice (shaded blue on the globe) and everywhere in the lowermost 500 km. c,d, 2D cross section (c) and 2D cross section with 3D rendering (d) of joint inversion of data used in a and b plus a global dataset of 189,439 multi-frequency P diff travel times (in five frequency bands; periods 30-7.5 s). e,f, 2D cross section (e) and 2D cross section with 3D rendering (f) of joint inversion of data used in c and d plus 55,657 multi-frequency teleseismic travel times measured on the RHUM-RUM array.

Fig. 3 |

 3 Fig. 3 | 3D rendering of slow P-velocity anomalies in model RRox-19. a, The same, 1,000-km-wide, cross-sectional band through La Réunion as Fig. 2f. Long-dashed purple line marks the top of the LLVP. The bullseye and cross mark the areas where Kerguelen and East African mid-mantle branches emerge out of and into the viewing plane. Structure outside the grey dashed lines is not 3D rendered for better visibility of the main features; the same is true for the grey area between points 3 and 4, which contains a smearing artefact emerging out of plane. b, An analogous cross section through the Bouvet mid-mantle branch, the LLVP cusp (purple dashed line), the East Africa mid-mantle branch and Afar hotspot (vertical line).

Fig. 4 |

 4 Fig.4| 3D rendering of seismically slow regions in the lower mantle, presumed to be upwelling. a, Isosurface of dV p /V p = -0.3% encloses slow material between 1,000 and 2,900 km, with colour indicating depth and changing every 100 km. Coastlines are filled light grey at the surface and outlined in dark grey at the deepest rendered depth (2,900 km). 3D rendered area is outlined in black, extending from 55º S, 20º W to 12º N, 50º E (at all rendered depths, in all panels). Hotspot locations are marked by black vertical columns from the surface to the deepest rendered depth: A, Afar; R, La Réunion; K, Kerguelen; M, Marion; B, Bouvet. b, Slow anomalies between 2,900 and 2,400 km depth, contoured at dV p /V p = -0.3%. A 2D isosurface is rendered every 100 km (five isosurfaces in total, the deepest filled with colour according to the colour bar). Other elements as described for panel a. c, Slow anomalies between 2,300 and 1,900 km depth. d, Slow anomalies between 1,600 and 1,300 km depth. Arrows highlight three diverging branches; they are not meant to imply trajectories of material flow. e, Slow anomalies between 1,100 and 900 km depth. As in d, arrows highlight diverging branches. f, Conceptual cartoon of the plume tree, its East Africa and SIO branches flattened into a plane.

Fig. 5 |

 5 Fig.5| Interpreted evolution of a tilted mid-mantle upwelling underlain by a low-velocity cMB corridor. a, A proto-plume, 1, ~800 km in diameter, develops at the distal end of an LLVP corridor, a region of maximum dT/dx. b, 1 detaches and rises vertically. The next instability, 2, develops left of 1, in the new instantaneous region of maximum dT/dx. c, Continued vertical rise of 1 and 2 and incipient detachment of 3 in the adjacent, new region of maximum dT/dx. d, Present-day snapshot shows tilted alignment of rising proto-plume blobs and a thinned low-velocity corridor left behind on the CMB. On reaching the viscosity boundary between lower and upper mantle, proto-plume 1 spawns the 'classical' sequence of a short-lived plume head (LIP), followed by one or more thin, long-lived plume tails, which remain rooted in their mid-mantle blob. Arrival of the plume head dynamically uplifted the lithosphere[START_REF] Cande | Indian and African plate motions driven by the push force of the Réunion plume head[END_REF] (vertically exaggerated). Proto-plume 1 corresponds with the Réunion plume in Fig.3(top) and Fig.4f; 3 corresponds with the LLVP cusp.

  2 and 3 render nested isosurfaces within 500 km of the CMB. The published tomographic models used for comparison in Extended Data Figs. 1, 3 and 5 were downloaded from our free SubMachine 73 web portal and visualized with ParaView using the same techniques. Comparison with other tomographic models. Extended Data Figs. 1, 3 and 5 compare our model RROx-19 with nine other global tomography models (from top to bottom and left to right): DETOX-P3

NATURe GeoScIeNce | www.nature.com/naturegeoscience

© The Author(s), under exclusive licence to Springer Nature Limited 2021NATURe GeoScIeNce | www.nature.com/naturegeoscience

Acknowledgements

RHUM-RUM was funded by Deutsche Forschungsgemeinschaft in Germany (grants SI1538/2-1, SI1538/3-1 and SI1538/4-1 and ship grant R/V Meteor cruise 101) and by Agence Nationale de la Recherche in France (project ANR-11-BS56-0013). Additional supports came from Centre National de la Recherche Scientifique-Institut National des Sciences de l'Univers (CNRS-INSU), Terres Australes et Antarctiques Françaises (TAAF), Institut Polaire Paul Emile Victor (IPEV), University of La Réunion and Alfred Wegener Institut (AWI). OBSs were provided by DEPAS (Deutsche Geräte-Pool für Amphibische Seismologie, Bremerhaven), GEOMAR Helmholtz-Zentrum für Ozeanforschung Kiel and INSU-IPGP (Institut National des Sciences de l'Univers -Institut de Physique du Globe de Paris). M.T. was funded by the NERC DTP in Environmental Research award NE/L002612/1 and a St. Edmund Hall RCUK Partnership award. K.S. acknowledges additional funding from the People Programme (Marie Curie Actions) of the European Union's Seventh Framework Programme FP7/2007-2013/ under REA grant agreement no. PCIG14-GA-2013-631104 'RHUM-RUM' and by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement 639003 'DEEP TIME'). K.H. acknowledges support by The Alan Turing Institute under the EPSRC grant EP/N510129/1. We thank M. Wysession for sharing MACOMO data. This is IPGP contribution 4097.

Data availability

All data used in this study are freely retrievable from international data centres that support the 'web service' protocol of the International Federation of Digital Seismograph Networks. We used the obspyDMT software of Hosseini and Sigloch [START_REF] Hosseini | ObspyDMT: a Python toolbox for retrieving and processing large seismological data sets[END_REF] to retrieve and pre-process all the waveforms used in this study. Our tomography model and the models compared with in this paper are freely available through the 'SubMachine' tomography web portal: http://submachine.earth.ox.ac.uk/ (where the user can select models, depths and other parameters).

code availability

The software used for tomographic inversion is available upon request from the corresponding author. It was customized from the software freely available at https://www.geoazur.fr/GLOBALSEIS/Soft.html.

Articles

Nature GeoscieNce SIO, including one of its currently active hotspots (La Réunion), and by much-improved resolution of the southern African LLVP. Surface relationships to geologically datable events (flood basalts, continental rift progression, dynamic topography, geochemistry) open new paths for investigating the dynamics of plumes, their source regions and their role in continental uplift and rifting.

Articles

Nature GeoscieNce

smearing occurs in the mid-mantle and near the CMB, where input anomalies are smaller than the inversion mesh (typically 400 km facet length for tetrahedra at the CMB). Supplementary Fig. 5 is a whole-mantle resolution test simulating vertical, whole-mantle plumes of the kind that we might have expected but did not observe under any Indian Ocean hotspot. This test shows that they would be well recovered if they were similarly thick as the mid-mantle branches that are observed instead. Note the good recovery at the CMB due to the inclusion of P diff measurements. Some lateral smearing appears below Kerguelen. In the South Atlantic, the model is not well resolved; most anomalies smear with low amplitudes or vanish completely under Tristan da Cunha and St. Helena. Supplementary Fig. 4 is a resolution test for upper-mantle structure under the RHUM-RUM instrumented area, using smaller chequers (500 km 2 diameter at surface). Recovery is good in areas where OBSs were deployed. Note how the resolution drops in areas with no stations, for example, in the northern part of cross-section 1 and the northeastern part of cross-section 2, where the recovered model smears more and dV/V amplitude drops. By contrast, in the east of cross-section 3 where OBS are present, we recover the chequers with little smearing. The tomography shows the bulk of the low-velocity anomalies not centred beneath La Réunion but slightly shifted to the southwest (coinciding with the expected progression direction of the hotspot track).

Supplementary Fig. 6 tests the extent of artificial upward smearing in the upper 400 km. The resolution input model is RROx-19 with all structure in the upper 400 km removed. The output shows upward smearing beneath the southern Mozambique Channel and the Central Indian Ridge but no vertical smearing beneath Réunion itself (the region most densely instrumented by OBS).

Supplementary Fig. 7 tests the extent of artificial downward smearing below 300 km. This is an ambitious test because body-wave tomographies tend to have poor ray crossing in the uppermost mantle, and the crustal model is never perfect. Both factors contribute to artificial downward smearing. input model is a surface-wave tomography model inverted from RHUM-RUM Rayleigh-wave measurements by Mazzullo et al. [START_REF] Mazzullo | Anisotropic tomography around La Réunion Island from Rayleigh waves[END_REF] , which extends from 0 to 300 km depth. At 100-200 km depth, our body-wave inversion recovers the overall structure of the surface-wave model well, but it does not resolve vertical transitions, such as the change from red (at 100 km) to blue (at 200 km) in the southeastern quadrangle of the model. Moderate artificial smearing is observed down to 400 km, and weak smearing down to 600 km. Smearing is more pronounced outside of the instrumented area.

Author contributions

M.T. analysed the RHUM-RUM body-wave data, modelled their sensitivities and developed the 3D visualization techniques. K.H. analysed the P-diffracted data and computed sensitivity kernels. M.T. and K.H. jointly integrated the datasets into a global-scale tomographic model. K.S. designed the tomography studies and acted as the primary academic adviser to M.T. and K.H. G.B. and K.S. designed the RHUM-RUM experiment and led the acquisition of its data. All authors contributed to writing this manuscript.

competing interests

The authors declare no competing interests.

Additional information

Extended data is available for this paper at https://doi.org/10.1038/s41561-021-00762-9.

Supplementary information

The online version contains supplementary material available at https://doi.org/10.1038/s41561-021-00762-9.

Correspondence and requests for materials should be addressed to M.T. [START_REF] Burke | Plume generation zones at the margins of large low shear velocity provinces on the core-mantle boundary[END_REF] . Southwestward away from the LLVP's core (reddest area), in direction of 7 o'clock, DETOX-P3 and all six S-wave models show a bulge towards Bouvet hotspot. This observation is consistent with the existence of a fully developed CMB corridor towards Bouvet, which would be recovered only as such a vestigial bulge (see targeted resolution test Extended Data Fig. 5) by our tomography, and presumably by others.

Peer review information

Articles

Nature GeoscieNce

Extended Data Fig. 5 | Seismic resolvability of three slow corridors at the cMB. Left: the resolution test input consists of three corridor-like, slow anomalies diverging from the LLVP centre under South Africa. Its geometries are idealised from lowermost mantle structure imaged by RROx-19. As in the tomography, corridors are ~700 km wide and 500 km 'high' above the CMB. The white column shows the location of Bouvet. Right: the results of this test confirm that the Southern Indian Ocean and East African corridors can be recovered well, while a Bouvet corridor would not be resolvable in our model, except for remnant southeastward bulging, away from the LLVP's core. Such a bulge is actually observed by all S-wave tomographies in Extended Data


Keywords: Autonomous Vehicles, Linear parametrically varying (LPV) methodologies, Robust control, Pole placement

This article presents the development of a control law for vehicle steering, acceleration and braking in the autonomous driving context. A self-scheduled robust feedback control design strategy is proposed for the stabilization of the vehicle along a predefined reference collision avoidance maneuver. A non-linear vehicle model with tire forces saturation is developed. Its linearization along the reference trajectory leads to a Linear Parameter-Varying (LPV) model. An LPV H ∞ control law with pole placement is computed for the LPV model. Simulation results illustrate the stabilization of the vehicle along the reference trajectory.

Self-scheduled H ∞ control of autonomous vehicle in collision avoidance maneuvers

Dario Penco * Joan Davins-Valldaura * * Emmanuel Godoy * Pedro Kvieska * * Giorgio Valmorbida *

INTRODUCTION

Autonomous driving is widely considered a possible solution to several private transportation challenges, as for example efficiency, energy consumption and safety, as described in [START_REF] Litman | Autonomous Vehicle Implementation Predictions: Implications for Transport Planning[END_REF]. Several Advanced Driver-Assistance Systems (ADAS) have been developed in the recent years, see [START_REF] Bengler | Three Decades of Driver Assistance Systems: Review and Future Perspectives[END_REF]. Some of these systems are able to control the vehicle under the driver supervision in specific environmental conditions. They constitute the first basic components of autonomous driving. To move towards increased levels of autonomous driving, defined in SAE (2018), less driver supervision is required. Hence safety functions have to be developed, to control the vehicle in emergency conditions. The development of methods for driving automation in collision avoidance maneuvers are of great interest for both academical and industrial communities.

The ADAS systems dealing with driving automation that are currently available on commercial vehicles mainly handle low dynamics maneuvers, i.e. maneuvers with a maximum vehicle acceleration of 0.3g. In these conditions, longitudinal and lateral vehicle dynamics can be considered uncoupled, and the tire forces behave linearly. The control of the longitudinal and lateral dynamics can hence be treated separately.

The longitudinal control architecture is usually designed hierarchically, with an higher level controller, generating the reference value for the acceleration, and a lower level controller that tracks the desired acceleration through throttle and brakes. A common solution for the higher level controller is the Linear Proportional Integral Derivative (PID) control law, as proposed in [START_REF] Persson | Stop and go controller for adaptive cruise control[END_REF], [START_REF] Naus | String-stable CACC design and experimental validation: A frequency-domain approach[END_REF] and [START_REF] Hima | Trajectory tracking for highly automated passenger vehicles[END_REF]. For the lower level controller the sliding mode method can be used for engine control to deliver the required torque, as proposed in [START_REF] Hedrick | Longitudinal Vehicle Controller Design for IVHS Systems[END_REF] and Sei-Bum [START_REF] Choi | Vehicle longitudinal control using an adaptive observer for automated highway systems[END_REF].

Lateral control is often based on the linear bicycle model, [START_REF] Rajamani | Vehicle dynamics and control[END_REF]. There exist several solutions that propose H 2 and H ∞ control laws, as in [START_REF] Shimakage | Design of lane-keeping control with steering torque input[END_REF] and [START_REF] Hima | Trajectory tracking for highly automated passenger vehicles[END_REF], to track lateral and yaw errors with respect to a reference trajectory. In high dynamics maneuvers, i.e. with accelerations greater than 0.3g, other phenomena occur. The longitudinal and lateral dynamics cannot be considered decoupled. Moreover the tire forces are subject to saturation and combined longitudinal and lateral slip. As a consequence in literature a large majority of the solutions developed to control the vehicle in high dynamics conditions consider vehicle models with coupled longitudinal and lateral dynamics.

There exist several solutions based on Model Predictive Control (MPC) for collision avoidance maneuvers. In [START_REF] Falcone | Predictive active steering control for autonomous vehicle systems[END_REF] an MPC controller is developed based on a non-linear bicycle model with both longitudinal and lateral dynamics. An alternative controller based on local linearizations of the non-linear model is also proposed to reduce the computational load. [START_REF] Falcone | MPC-based yaw and lateral stabilisation via active front steering and braking[END_REF] and [START_REF] Yin | Integration of motion planning and model-predictivecontrol-based control system for autonomous electric vehicles[END_REF] propose a similar controller based on a non-linear four-wheel vehicle model. In [START_REF] Gao | A tube-based robust nonlinear predictive control approach to semiautonomous ground vehicles[END_REF] tube-based MPC is used to address also robustness criteria.

The main advantage of MPC is the possibility to take into consideration the inputs and states constraints in the underlying optimization problem. It is then easy to represent actuators limits. However the framework requires a large computational load. This is a crucial point for the integration of the controller on a vehicle Electronic Control Unit (ECU), which for economic reasons has a limited computational capacity. On the contrary, the Linear Parameter-Varying (LPV) framework offers methods for the synthesis of controllers that do not require a large computation load; see [START_REF] Wu | Control of Linear Parameter Varying Systems[END_REF], [START_REF] Biannic | Commande robuste des systèmes a paramètres variables[END_REF], [START_REF] Blanchini | Set-Theoretic Methods in Control[END_REF]. Moreover it has already been proven effective in the vehicle dynamics control domain. In literature there exist several solutions that employ LPV controllers for Global Chassis Control (GCC) systems, which aim is either assist the driver by increasing vehicle stability in high dynamics conditions or improve comfort through several chassis actuators. In [START_REF] Fergani | Robust multivariable control for vehicle dynamics[END_REF], [START_REF] Gáspár | Towards global chassis control by integrating the brake and suspention system[END_REF] and [START_REF] Poussot-Vassal | Vehicle dynamic stability improvements through gain-scheduled steering and braking control[END_REF] different solutions for vehicle stabilisation through combinations of active steering, braking and active suspensions are proposed. They are all based on robust LPV H ∞ control law. This paper proposes a self-scheduled LPV H ∞ controller for vehicle stabilization along a reference collision avoidance maneuver. The controller is based on a LPV model for the vehicle dynamics along the maneuver, presented by the authors in [START_REF] Penco | Control for autonomous vehicles in high dynamics maneuvers : LPV modeling and static feedback controller[END_REF]. This LPV model is obtained from the linearization along the reference maneuver of a nonlinear model that considers the coupling between longitudinal and lateral dynamics and the tire nonlinear behavior. In particular, due to an algebraic loop introduced by the saturation function of the tires, we propose a strategy to obtain an uncertain parameter related to the slope of sector non-linearities.

The paper is organized as follows. In Section 2 the nonlinear vehicle model and the procedure to obtain the LPV model are discussed. The method used for the synthesis of the controller is illustrated in Section 3. In Section 4 the results of the application of the controller in a simulation environment are shown. The conclusions and future work are discussed in Section 5.

PRELIMINARY RESULTS

In this section we discuss the development of a non-linear vehicle model and the procedure to obtain an LPV model.

A more detailed description of the non-linear model is available in [START_REF] Penco | Control for autonomous vehicles in high dynamics maneuvers : LPV modeling and static feedback controller[END_REF]. The difference between the model presented in [START_REF] Penco | Control for autonomous vehicles in high dynamics maneuvers : LPV modeling and static feedback controller[END_REF] and the model studied in this paper is the influence of the road slope and camber in the dynamical model.

Non-linear vehicle model

To obtain a model of the vehicle we consider that the two wheels of each axle have been regrouped at the center of the axle track. The vehicle degrees of freedom allow it to move on the xy plane, so it can translate along the x and y axis and rotate around the z axis. The vehicle axis are shown in Figure 1.

In Table 1 there are the model state variables and inputs.

The road slope and camber angles are shown in Figure 2. Table 2 contains the forces shown in Figures 1 and3. The forces applied on the front wheel have each a component on the x and y axis of the vehicle, that depends on δ f . We denote the components of F xf on the xy plane:

F xxf = F xf cos δ f F yxf = F xf sin δ f .
The model equations of motion in the vehicle frame are:

m ( v -r u) = 2 F xxf + Fxr -R xxf -Rxr -F xyf -Faero -Px (1a)
m ( u + r v) = 2 F yxf -R yxf + F yyf + Fyr -Py (1b) Izz ṙ = 2 F yxf + F yyf -R yxf f -2 Fyr r (1c) Iw y ωwf = τ wf -2 F xf re (1d)
Iw y ωwr = τwr -2 Fxr re.

(1e) In high dynamics maneuvers, the tires are subject to saturation and combined slip [START_REF] Pacejka | Tyre and vehicle dynamics[END_REF]). It is important then to take into consideration these phenomena by using a non-linear model for the tire forces. The tire model in [START_REF] Kissai | A new linear tire model with varying parameters[END_REF], developed for control purposes, is a variation of the Dugoff 's model, proposed in [START_REF] Dugoff | Tire performance characteristics affecting vehicle response to steering and braking control input[END_REF]. It is simpler and requires less parameters than other commonly used models, as the Pacejka's magic formula, in [START_REF] Bakker | Tyre modelling for use in vehicle dynamics studies[END_REF], and the brush model, in [START_REF] Svendenius | Tire Modeling and Friction Estimation[END_REF]. We rely on this model since it allows to represent tire forces saturation and combined slip. The longitudinal and lateral forces are linear respectively on the longitudinal and lateral slip, κ and α as

Fx = c * κ (µ, N, α) κ (2a) Fy = c * α (µ, N, κ) α. ( 2b
)
Longitudinal and lateral slip κ and α are defined in [START_REF] Pacejka | Tyre and vehicle dynamics[END_REF]. The coefficients c * κ and c * α are non-linear with respect to the α and κ respectively, the wheel normal force N and the tire-ground friction coefficient µ.

For the sake of brevity we will not discuss the tire model in details. It is however important to discuss of the tire force saturation. Here we use a circle to approximate the saturation of the tire forces as in [START_REF] Schuring | A model for combined tire cornering and braking forces[END_REF], thus simplifying the more general ellipse model, discussed in [START_REF] Svendenius | Tire Modeling and Friction Estimation[END_REF]. The radius of the saturation circle is given by F max = µN . The logistic function, shown in Figure 4, is used to express the saturation

σ F , u, = u - 1 + e -4 u-( F -u- 2 ) + (3)
where u and are the upper and lower saturation bounds.

The unsaturated forces Fx and Fy of ( 2) are the inputs of (3), the outputs are the saturated tire forces F x and F y .

The force N on the tire is present both in the coefficients c * κ and c * α and the saturation bound. Considering the importance of N in the computation of the tire forces, a variable normal force model has then been developed. The procedure to obtain it has been discussed in [START_REF] Penco | Control for autonomous vehicles in high dynamics maneuvers : LPV modeling and static feedback controller[END_REF] and for brevity is not shown in this paper. The expressions of N f and N r depend on the longitudinal acceleration v. N f and N r are saturated between 0 and P z using the logistic function. To be brief the other forces are not discussed in this paper, if needed the reader can consult [START_REF] Rajamani | Vehicle dynamics and control[END_REF]. Finally we complete the vehicle model by adding the vehicle position and orientation. Their dynamics on the inertial frame are governed by:

ẋ = v cos ψ -u sin ψ (4a) ẏ = u cos ψ + v sin ψ (4b) ψ = r.
(4c)

Model linearization

The state and input vectors of the model in (1) are:

Fig. 4. Logistic function. x =     v u r ω wf ωwr     u =     δ f τ wf τwr θx θy     .
The model expressed by ( 1) is in implicit form f ( ẋ, x, u) = 0. Due to the saturation functions used for the tire forces, it is not possible to express the model in explicit form. We can write the model in a feedback form as in Figure 5: ẋ = g (x, u) + B σ (u) σ (h) (5) where: The vectors h ( ẋ, x, u, σ) and σ (h) contain respectively the inputs and the outputs of the saturation functions. The model non-linearities due to the saturation functions have then been isolated and they result as an input of the model. The expressions of g (x, u) and B σ (u) for brevity are not detailed here. The state dynamics expressed as in ( 5) is then linearized along a reference trajectory ( ẋ0 (t), x 0 (t), u 0 (t), σ 0 (t)), with σ 0 (t) = σ(h( ẋ0 (t), x 0 (t), u 0 (t))). For simplicity, in the following, the reference trajectory, containing both the states and the inputs, will be denoted by T (t).

σ (h) =       σ (h 1 ( ẋ, x, u)) σ (h 2 ( ẋ, x, u)) σ (h 3 ( ẋ, x, u, σ)) σ (h 4 ( ẋ, x, u, σ)) σ (h 5 ( ẋ, x, u, σ)) σ (h 6 ( ẋ, x, u, σ))       h ( ẋ, x, u, σ) =       h 1 ( ẋ, x, u) h 2 ( ẋ, x, u) h 3 ( ẋ, x, u, σ) h 4 ( ẋ, x, u, σ) h 5 ( ẋ, x, u, σ) h 6 ( ẋ, x, u, σ)       =        Nf Nr Fxf σ Nf Fxr σ Nr Fyf σ Nf Fyr σ Nr        .
∆ ẋ = A(t)∆x + B(t)∆u + Bσ(t)∆σ ∆h = C(t)∆x + D(t)∆u + Dσ(t)∆σ A(t) = ∂g (x, u) ∂x T (t) B(t) = ∂g (x, u) ∂u T (t) + ∂ (Bσ (u) σ (h)) ∂u T (t) Bσ(t) = Bσ (u)| T (t) C(t) = ∂h ( ẋ, x, u, σ) ∂ ẋ T (t) A(t) + ∂h ( ẋ, x, u, σ) ∂x T (t) D(t) = ∂h ( ẋ, x, u, σ) ∂ ẋ T (t) B(t) + ∂h ( ẋ, x, u, σ) ∂u T (t) Dσ(t) = ∂h ( ẋ, x, u, σ) ∂ ẋ T (t) Bσ(t) + ∂h ( ẋ, x, u, σ) ∂σ T (t)
In this paper we consider the reference trajectory of Figure 6. It corresponds to a two lanes change maneuver without braking, starting at the speed of 70km/h. The saturation functions can be approximated as sector- bounded non-linearities, whose slope varies according to the system working point. Figure 7 shows the slope sector for each element of σ (h) along reference trajectory. The wheels torques are equal to zero, hence the longitudinal tire forces and the load transfer are small and their sectors are then narrow. On the contrary the lateral tire forces magnitude in the reference trajectory includes a wider sector.

-1 ,5 0 0 -1 ,0 0 0 -5 0 0 0 5 0 0 1 ,0 0 0 1 ,5 0 0 We then approximate the non-linear function by a timevarying affine function as σ (h) ≈ K σ h. K σ is a timevarying diagonal matrix, with the slope of the saturation functions along the diagonal. Applying this approximation in the expression of ∆h in ( 6) we obtain:

-1,000 0 1,000 Nf σ Nf -1 ,5 0 0 -1 ,0 0 0 -5 0 0 0 5 0 0 1 ,0 0 0 1 ,5 0 0 -1,000 0 1,000 Nr σ Nr -1 ,5 0 0 -1 ,0 0 0 -5 0 0 0 5 0 0 1 ,0 0 0 1 ,5 0 0 -1,000 0 1,000 Fxf σ Fxf -1 ,5 0 0 -1 ,0 0 0 -5 0 0 0 5 0 0 1 ,0 0 0 1 ,5 0 0 -1,000 0
∆h = (I -D σ (t) K σ ) -1 C(t)∆x + (I -D σ (t) K σ ) -1 D(t)∆u.
Replacing this expression of ∆h in ( 6), the linearized state dynamics with the sector-bounded approximation of the saturation functions is:

∆ ẋ = A(t)∆x + B(t)∆u (7a) A(t) = A(t) + B σ (t) K σ (I -D σ (t) K σ ) -1 C(t) (7b) B(t) = B(t) + B σ (t) K σ (I -D σ (t) K σ ) -1 D(t). (7c)

LPV polytopic model for control synthesis

Given that the linearized model is obtained along a reference trajectory and the saturation functions are approximated by time-varying gains, the matrices A(t) and B(t) in ( 7b) and (7c) are time-varying.

In the following, the elements of the matrix K σ will be considered as the average slope coefficients of the slope sector shown in Figure 7. It is possible to consider the elements of A(t) and B(t) that vary with time as varying parameters. We call θ i , i = 1, ... , q these parameters. p is the vector containing all the parameters. We then express the linearized model in ( 7) as an LPV model, with the state matrix A and the input matrix B affinely dependent on the varying parameters.

Along the reference trajectory shown in Figure 6, a total of 23 elements of A(t) and B(t) vary with time. However several of these elements have a small variation. It is possible to consider just the 6 elements of A(t) and B(t) that vary the most to constitute the vector p of parameters for synthesis purposes. The range of each parameter θ i , i = 1, ... , 6 is between a minimum p i and a maximum p i . p i and p i are the minimum and maximum values of p i along the reference trajectory T (t).

For the control synthesis we add the longitudinal, lateral and yaw angle errors with respect to T (t) to the state of the model. They are defined as:

x L = x -x 0 y L = y -y 0 ∆ψ = ψ -ψ 0 .
x 0 , y 0 and ψ 0 are the vehicle's absolute position and attitude along the reference trajectory T (t). The vehicle's position dynamics in (4a) and (4b) can be simplified neglecting the vehicle's lateral speed u, which is considerably smaller than the longitudinal speed v ẋ = v cos ψ (9a) ẏ = v sin ψ.

(9b) Considering that v = v 0 +∆v, u = u 0 +∆u and ψ = ψ 0 +∆ψ and that the ∆ψ is small along the reference trajectory, the dynamics of the error variables are: ẋL = ∆v ẏL = ∆v ψ 0 + ∆ψ v 0 cos ψ 0 ∆ ψ = ∆r. The expressions ψ 0 and v 0 cos ψ 0 are added to p. It is finally possible to define the polytopic domain P as the hypercube of dimension 8, each dimension being the parameter p i , i = 1, ... , 8, between p i and p i . The polytopic model is defined at its N = 2 8 vertices:

S i = ( A i B i ) i = 1, ... , N.
(11) The model inside P is defined as:

S p = N i=1 α i ( A i B i )
where α i for i = 1, ... , N are the interpolation parameters and N i=1 α i = 1.

MAIN RESULTS

A dynamic LPV H ∞ state feedback controller has been designed starting from the polytopic LPV model described in Section 2.3. We consider the system of order n in the form:

ẋ = A i x + B 1i w + B 2 u (12a) z = C 1i x + D 11i w + D 12 u (12b) y = C 2 x + D 21 w (12c)
for i = 1, ... , N , where y is the controller input, z contains the controlled variables, w the exogenous inputs and u the control inputs. The system in the polytopic domain P is then expressed with respect to the N vertices:

S p = N i=1 α i A i B 1i B 2 C 1i D 11i D 12 C 2 D 21 0 (13)
where α i for i = 1, ... , N are the interpolation coefficients.

In this article, the solution proposed in [START_REF] Chilali | H ∞ Design with Pole Placement Constraints: an LMI Approach[END_REF] is used for the synthesis of the LPV H ∞ controller. This method allows to specify limits to the dynamics of the closed loop. This is obtained by computing a full-order LPV controller K p that forces the closed loop eigenvalues into some LMI region for all admissible values of (13) in P.

An LMI region D = z ∈ C : L + zM + zM T < 0 of order s is a subset of the complex plane, with L = L T ∈ R s×s and M ∈ R s×s . Given a polytopic plant ( 12), an region D and the H ∞ performance γ > 0, we search an LPV controller such that: (A) the closed loop poles for all admissible values of (13) in P are in D; (B) the L 2 -norm of the input/output mapping from w to z is bounded by γ. It exists if there exist two symmetric matrices R, S ∈ R n×n and matrices A K , B K , C K and D K such that:

R I I S > 0 Ψ 11 Ψ T 21 Ψ 21 Ψ 22 < 0 L ⊗ R I I S + M ⊗ Φ + M T ⊗ Φ T < 0
where, for i = 1, ... , N

Φ = A i R + B 2 C K A i + B 2 D K C 2 A K SA i + B K C 2 Ψ 11 = A i R + RA T i + B 2 C K + C T K B T 2 B 1i + B 2 D K D 21 (B 1i + B 2 D K D 21 ) T
-γI

Ψ 21 = A K + (A i + B 2 D K C 2 ) T SB 1i + B K D 21 C 1i R + D 12 C K D 11i + D 12 D K D 21 Ψ 22 = A T i S + SA i + B K C 2 + C T 2 B T K (C 1i + D 12 D K C 2 ) T C 1i + D 12 D K C 2 -γI
We compute any square matrices M and N such that M N T = I -RS. Solving the system of equation:

B K = N B Ki + SB 2 D K C K = C Ki M T + D K C 2 R A K = N A Ki M T + N B Ki C 2 R + SB 2 C Ki M T + S (A i + B 2 D K C 2 ) R for i = 1, ... , N , for B Ki , C Ki
and A Ki we find the controller matrices at the vertices of the polytopic domain P. The full-order LPV H ∞ controller is:

K p = N i=1 α i A Ki B Ki C Ki D K .
For the control synthesis an LMI region has been used. It guarantees that the closed loop eigenvalues for all admissible values in P are between λ = -8500 and λ = -1.8, to bound slow and fast dynamics, and θ = 1.12rad, to have a minimum damping ratio ζ > 0.9.

To bound the bandwidth of the control signals a diagonal weighting function W u as been used, where the elements on the diagonal are W u1 (s) = 100s+31.57 s+315.7 and W u2 (s) = W u3 (s) = 100s+315.7 s+3157 . The control architecture is represented in Figure 8. T , y = x and:

x =        ∆v ∆u ∆r ∆ω wf ∆ωwr x L y L ∆ψ        w = ∆θx ∆θx u = ∆δ f ∆τ wf ∆τwr z 1 = x L y L ∆ψ . (17)
The matrix B 1i and B 2 are obtained from B i in (11), given the definitions of w and u in (17).

NUMERICAL RESULTS

In this Section some numerical results of the controller are shown. The controller has been tested on MAT-LAB/Simulink on the non-linear bicycle model described in Section 2.1. The trajectory considered is the reference trajectory shown in Section 2.2, Figure 6. Several simulations have been carried out, with a set of different initial conditions for ∆v and ∆u, to establish empirically the region of initial conditions for which the controller is able to stabilize the system and converge to the reference trajectory. Figure 9 shows the empirical region of attraction of the controller for non-zero initial conditions for the variables ∆v and ∆u. The capacity of the controller to stabilize the system and bring back the vehicle along the reference trajectory is limited for -5.2m/s ∆v (0) 4.3m/s and -5m/s ∆u (0) 8.4m/s. A case with initial conditions of ∆v (0) = 3m/s, ∆u (0) = 2m/s and ∆r (0) = 15deg/s is simulated and shown in detail. Figure 10 shows the commands computed by the controller, the vehicle's position and attitude error and the vehicle overall longitudinal and lateral accelerations. Due to the large initial conditions, the position errors are initially large, but they eventually converge to zero. Indeed the vehicle acceleration belong to the high dynamics category.

CONCLUSION AND PERSPECTIVES

The LPV model discussed in Section 2 had already been used by the authors for the synthesis of a static state feedback controller in [START_REF] Penco | Control for autonomous vehicles in high dynamics maneuvers : LPV modeling and static feedback controller[END_REF]. The results of Section 4 show that the LPV H ∞ dynamic controller developed in Section 3, compared to the static controller developed in [START_REF] Penco | Control for autonomous vehicles in high dynamics maneuvers : LPV modeling and static feedback controller[END_REF], considerably enlarges the range of initial conditions for which the controller is able to stabilize the system along the reference trajectory. Moreover this has been obtained considering also the road inclination as system perturbation. A more thorough methodology to estimate the region of attraction of the closed-loop systems, without simulations, will be proposed by considering sector inequalities associated to the saturation function.

Fig. 1 .

 1 Fig. 1. Vehicle's xy plane, with the forces listed in Table2.

Fig. 2 .

 2 Fig. 2. Road inclination.

Fig. 5 .

 5 Fig. 5. Isolation of the implicit part of (1).

  Fig. 6. Reference trajectory.

Fig. 7 .

 7 Fig. 7. Slope sector for reference trajectory of Figure 6.

Fig. 8 .

 8 Fig. 8. Control architecture. We define z = [z 1 , z 2 ]T , y = x and:

Fig. 9 .

 9 Fig. 9. Empirical region of attraction for ∆v and ∆u.

  Fig. 10. Simulation results.

Table 1 .

 1 Model state variables and inputs

			State variables
	v	[m/s]	longitudinal speed
	u	[m/s]	speed
	r	[rad/s]	yaw rate
	ω wf	[rad/s]	front wheel rotational speed
	ωwr	[rad/s]	rear wheel rotational speed
	x	[m]	absolute longitudinal position
	y	[m]	absolute lateral position
	ψ	[rad]	yaw angle
			System's inputs
	δ f	[deg]	front wheel steering angle
	τ wf	[Nm]	front wheel torque
	τwr	[Nm]	rear wheel torque
	θx	[rad]	road camber angle
	θy	[rad]	road slope angle

Table 2 .

 2 Model's forces

	F xf	front tire longitudinal force	Fxr	rear tire longitudinal force
	F yf	front tire lateral force	Fyr	rear tire lateral force
	R xf	front wheel rolling friction	Rxr	rear wheel rolling friction
	Faero	aerodynamic drag	Pz	weight component in z axis
	Px	weight component in x axis	Py	weight component in y axis
	N f	front wheel normal force	Nr	rear wheel normal force