Matthieu Roy-Barman

François Thil

Louise Bordier

Arnaud Dapoigny

Lorna Foliot

Sophie Ayrault

François Lacan

Catherine Jeandel

Catherine Pradoux

Ester Garcia-Solsona

Thorium isotopes in the Southeast Atlantic Ocean: Tracking scavenging during water mass mixing along neutral density surfaces

Keywords: thorium isotopes, Southern Ocean, Atlantic, marine particles, isopycnal mixing Highlights

The distributions of dissolved and particulate thorium isotopes ( 230 Th and 232 Th) were established in samples from the BONUS GoodHope (BGH) IPY-GEOTRACES cruise in the SE Atlantic sector of the Southern Ocean (36°S-13°E to 57°S-0°, Feb.-Mar. 2008). The distribution of total (dissolved+particulate) 232 Th is dominated by the inputs from continental margins. The non-linear profiles of dissolved 230 Th are interpreted as due to the southward upwelling of the isopycnal surfaces. However, total 230 Th and 232 Th versus salinity plots illustrate departures from binary mixing and provides evidence for non-conservative behavior of both isotopes along the section. We propose a model for total 230 Th and 232 Th scavenging and mixing along isopycnal surfaces. We use this model to estimate particle settling speeds and isopycnal eddy diffusion coefficients along the BGH section. Data-model comparison suggests particle settling velocities in the range of 400-700 m/y and isopycnal eddy diffusivity of the order of 2000 m 2 /s.

Introduction

Thorium isotopes are recognized tracers of particle dynamics in the ocean that provide key information on ocean biogeochemical cycles [START_REF] Bacon | Distribution of Thorium isotopes between dissolved and particulate forms in the deep sea[END_REF][START_REF] Nozaki | The water column distribution of thorium isotopes in the western North Pacific[END_REF][START_REF] Roy-Barman | The sources and the fates of thorium[END_REF]. In particular, 230 Th is produced in situ by radioactive decay of 234 U and rapidly scavenged on marine particles, thus providing a chronometer for estimating the settling speed for marine particles. 232 Th is introduced in the ocean by continental inputs and hence can bring complementary information on marine particles compared to 230 Th (Roy- [START_REF] Roy-Barman | Thorium isotopes in the western Mediterranean Sea: An insight into the marine particle dynamics[END_REF].

In early studies, particle settling speeds were derived from 230 Th with a simple 1D (vertical) production-scavenging model, which was well adapted to oceanic areas with weak currents [START_REF] Bacon | Distribution of Thorium isotopes between dissolved and particulate forms in the deep sea[END_REF][START_REF] Nozaki | The water column distribution of thorium isotopes in the western North Pacific[END_REF][START_REF] Roy-Barman | The sources and the fates of thorium[END_REF][START_REF] Roy-Barman | Modelling the effect of boundary scavenging on Thorium and Protactinium profiles in the ocean[END_REF]. However, it is now clear that strong currents characterizing oceanic regions such as the Southern Ocean may significantly affect the distribution of 230 Th in the water column and hence potentially bias the estimate of the settling speed of marine particles obtained by the simple production-scavenging model. Therefore, a simple productionscavenging-advection model was proposed to constrain simultaneously the renewal rate of the deep waters and the settling speed of marine particles in the Atlantic sector of the Southern Ocean (Rutgers van der Loeff et al., 1996). However, it relies on rather crude approximations, such as uniform values for ventilation rate and particle settling speed over the whole water column. Another limitation of this model is that it does not consider the mixing of water masses with different origins and Th contents.

In the real ocean, the distribution of geochemical tracers depends on both nonconservative processes and mixing of water masses of different origins, in addition to advection by currents. This mixing is thought to occur preferentially along surfaces of constant density called isopycnal surfaces. [START_REF] Chase | Scavenging of 230 Th, 231 Pa and 10 Be in the Southern Ocean (SW Pacific sector): The importance of particle flux, particle composition and advection[END_REF] suggested that the small gradients in 230 Th concentration observed along neutral density surfaces in the deep waters of the Southern Ocean argue for small eddy diffusive fluxes and little scavenging of 230 Th across the Antarctic Polar Front (APF). However, they also advocated for further work to better constrain both nuclide concentrations along isopycnal surfaces and diffusion coefficients in the Southern Ocean.

In the present work, we take advantage of high precision 230 Th data collected during the Bonus GoodHope cruise to revisit the behavior of 230 Th across the Southern Ocean and to try to unravel the respective roles of water mass mixing and particle scavenging along neutral density surfaces. Our goal is to explore if Th data can help constraining both particle settling speeds and water mixing rates using a model that is more complete than with the simple mixing-scavenging model.

Sampling and methods

Hydrographic settings

The oceanographic structure determined from the Bonus GoodHope (BGH) expedition has already been described [START_REF] Chever | Physical speciation of iron in the Atlantic sector of the Southern Ocean along a transect from the subtropical domain to the Weddell Sea Gyre[END_REF][START_REF] Bown | The biogeochemical cycle of dissolved cobalt in the Atlantic and the Southern Ocean south off the coast of South Africa[END_REF][START_REF] Garcia-Solsona | Rare earth elements and Nd isotopes tracing water mass mixing and particleseawater interactions in the SE Atlantic[END_REF][START_REF] Abadie | Iron isotopes reveal distinct dissolved iron sources and pathways in the intermediate versus deep Southern Ocean[END_REF], so we only summarize the main features and focus on the 5 super stations (S1-S5 stations, Fig. 1) sampled for Th isotopes. Three main hydrological and circulation regimes characterize the studied area: the Subtropical domain (Station S1), the Antarctic Circumpolar Current (ACC; S2 to S4) and the Weddell Sea Gyre (S5). The subtropical domain is bounded to the south by the subtropical front (STF), which was located around 41°S-42°S, north of S2. The ACC domain, incorporating the Subantarctic Front (SAF; ~45°S) and the Polar Front (PF; ~50°S), is bounded to the south by the southern ACC boundary (Sby). From the Sby to the Antarctic continent, water circulation is dominated by the cyclonic Weddell Gyre (WG; S5).

The subtropical regime is influenced from 0 down to ~800 m by surface and intermediate waters of Indian Ocean origin ( ≥ 10 °C, S ≥ 34.5), conveyed south of Africa by the Agulhas Current and detached eddies [START_REF] Gladyshev | A hydrographic section from South Africa to the southern limit of the Antarctic Circumpolar Current at the Greenwich meridian[END_REF]. At intermediate depths (~300-900 m), Antarctic Intermediate Water (AAIW) originates from cooled surface waters of the ACC that are subsequently subducted northward at the PF. Two varieties of AAIW have been recognized from BGH cruise data: one originating from the Indian Ocean (i-AAIW) close to the African continent (S > 34.3; S1) and one formed in the Atlantic Ocean (a-AAIW; S < 34.3; S2-S4). The Upper Circumpolar Deep Water (UCDW; ranging from 500 to 2000 m depth) is identified below the AAIW and upwells southwards. Below UCDW, the eastward flowing (and southward drifting) NADW is detected by its salinity maximum within the 2000-3000 m depth range. Two NADW branches are identified in the BONUS section: North of the Agulhas Ridge, higher salinities indicate advection by an eastern boundary current referred to as SE-NADW (S1 and S2 stations). South of the Agulhas Ridge, the main (western) component of NADW (i.e., SW-NADW; S2-S4 stations) is incorporated into the ACC after flowing southward along the Argentine-Brazilian margin (Stramma and England, 1999). The other water mass encountered at southern ACC latitudes is high salinity Lower Circumpolar Deep Water (LCDW). Cold and fresher Antarctic Bottom Water (AABW) is also recognized, occupying the water column from LCDW/NADW depths down to the sea floor (Fig. 2). In the Weddell Sea Gyre domain (S5), we interpret deep oxygenated waters as due to input of recently ventilated AABW.

Water mass transport along the BGH section is mainly zonal. In the subtropical domain (station S1), it is dominated by the westward flowing Agulhas Current in the upper water column (z < 1500 m) and at depth by the eastward transport of SE-NADW as delineated by [START_REF] Sokolov | Circumpolar structure and distribution of the Antarctic Circumpolar fronts: 1. Mean circumpolar path[END_REF]. Stations S2, S3 and S4 are located within the Antarctic Circumpolar Current, flowing eastward throughout the water column. Finally station S5 is located in the Northern limb of the Weddell gyre, also flowing eastward. Therefore, the BGH section is approximately perpendicular to the climatological mean (dominant) flows, so that water masses should not flow directly from one BGH station to another. The different BGH stations would therefore not be directly connected by the mean flow. However, superimposed on this vigorous zonal circulation, slower meridional transports are clearly visible on meridional sections of water mass tracers such as salinity, temperature or nutrients, revealing notably the northward transport of AAIW and AABW, and southward transport of NADW [START_REF] Abadie | Iron isotopes reveal distinct dissolved iron sources and pathways in the intermediate versus deep Southern Ocean[END_REF]. This meridional transport is presumably a residual circulation in which lateral eddy fluxes largely balance the wind-driven circulation (Marshall and Speer, 2016). While the circulation in the BGH area is therefore clearly 3D and dominated by zonal flows (preventing a direct link by mean advection between our 5 stations), some of the BGH samples can still reflect the effect of these meridional transports of these water masses upstream of the BGH section. In the following, transport models will be used to link tracer properties measured along the BGH section. We do not assume that the transports actually act along the BGH section (pure meridional transport is not realistic). However, we use the tracer properties measured along the BGH section to elucidate and quantify the meridional component (advection, diffusion) embedded within a dominant zonal component and therefore acts upstream of the BGH section, for instance between the Drake Passage and the BGH section.

Sampling

The BONUS GoodHope cruise took place in the framework of the international polar year (IPY) from February 8 th 2008 till March 17 th 2008 on board the Research Vessel Marion-Dufresne II, from 34°S 18°E to 51°S 00°W and from 51°S to 58°S along the Greenwich meridian. The expedition took place during the late austral summer, when biological activity is high and ice cover is reduced. Five depth profiles (S1-S5 stations) were sampled for both dissolved and particulate Th with 10 or 11 depths per profile with 1 duplicated sample at stations S3, S4 and S5 (Fig. 1). An intercalibration (hereafter referred to as 'intercal') station was analyzed for 3 depths which were all duplicated.

Seawater samples were collected using Niskin bottles equipped with Teflon coated strings and mounted on a rosette frame with a CTD sensor. Challenger Oceanic in situ pumps (ISP) were used to filter large seawater volumes through SUPOR filters (pore size: 0.45 m, filter diameter: 293 mm). Typically, volumes ranging from 200 L to 1000 L were filtered.

Analytical procedure

Filtered seawater

Twenty liters of filtered seawater (Nuclepore, 90 mm diameter, 0.4 m pore size) were acidified, spiked with 229 Th and pre-concentrated on board using the ferric hydroxide coprecipitation technique as described in Piepgras et al. (1979). The precipitate was then recovered on filters (Supor, 90 mm, 0.4 m), which were kept in petri dishes and stored in the dark. Back in the laboratory, the filters were submerged in 40 mL of 6 M HCl for 24 h in 300 mL Polytetrafluoroethylene beakers, enabling the dissolution of precipitate. The filters were then carefully rinsed with deionized water to recover all the dissolved precipitate. After evaporation, the residue was dissolved in 0.5 mL 7M HNO3 and loaded on an anion-exchange column (AG1-X8 resin, 200-400 mesh) to separate the thorium fraction from the major elements, uranium and the rare earth elements (REE) (adapted from [START_REF] Jeandel | Single column sequential extraction of Ra, Nd, Th, Pa and U from a natural sample[END_REF]. The total 232 Th chemical blank is typically 48 ± 14 pg and represents generally less than 15% of the sample although peaks with 50% of the sample occur for very 232 Th-poor samples. Total 230 Th chemical blank is typically 2.3 ± 0.1 fg and represents generally less than 5 % of the sample although peaks with 30% of the sample occur for very shallow samples.

Particles collected from Niskin bottle filtrations

From the twenty liter seawater samples used for dissolved Th analyses (described above), particulate samples were also available on the 90mm diameter 0.4µm porosity Nuclepore membranes used for filtration (hereafter referred to as Niskin filtration). Membranes were transferred from their petri dishes to Teflon beaker containing 15mL 6M HCl and 2.5mL 15M HNO3. They were heated 90 minutes at 130°C. After leaving them cool, 0.5mL 23M HF was added, and the mixture was heated again as previously. The membranes were removed from the leaching solutions and rinsed over the beakers with milliQ water. After evaporation, the samples were dissolved in 3mL 0.32M HNO3 for analysis by Quadrupole Induced Coupled Plasma Mass Spectrometry.

Particles collected from in situ pump filtrations

Filters from the in situ pumps were cut into pieces using ceramic scissors. Several leaching steps were necessary to recover Th without attacking the filters. Filter pieces were first leached with 200 mL of 6N HCl and 1.5 mL of 24-25N HF in a 300 mL Teflon beaker for 2 days at 75°C. The filter pieces were removed from the leaching solution and rinsed with a MQ water squeeze bottle over the leaching beaker. The filter pieces were saved apart. The leaching and rinsing solutions were evaporated down to ~10 mL and transferred into a 30 mL Teflon beaker. Then, the filter pieces were leached a second time with 150 mL of 7N HNO3 and 0.05 mL of 29N HF for 2 days at 75°C. The filter pieces were removed from the leaching solution and rinsed with MQ water that was again recovered into the leaching beaker. After evaporation to a few mL, this solution was added to the first leaching/rinsing solution. The 300 mL Teflon beaker was rinsed in warm diluted HNO3 to remove any particle sticking on the beaker walls and the resulting solution was also added to the 30 mL beaker. The resulting solution was then spiked with 229 Th (and 233 Pa), dried and dissolved again in a solution with 4 mL of 14N

HNO3and 1 mL of 12N HCl. After 1 night on a hot plate, the solution was dried and the residue was dissolved again in 10 ml of 1N HNO3. Since filters were not rinsed immediately after filtration with distilled water on board, it was preferable to remove the salt before the column chemistry. Therefore, 40µL of a Fe solution (60 mg/g) were added. After 1h-heating, Fe was precipitated by raising the pH to about 8 with NH3. The Fe precipitate (that coprecipitates Th isotopes) was separated by centrifugation and rinsed several times. Finally, it was dissolved in 0.25 ml of 8N HNO3, ready for loading on an anionic column. Th isotopes were separated from

Fe by ion exchange chromatography on a small volume (0.5 ml) column of AG1X8 resin (200-400 mesh) with a procedure adapted from [START_REF] Jeandel | Single column sequential extraction of Ra, Nd, Th, Pa and U from a natural sample[END_REF][START_REF] Gdaniec | Thorium and protactinium isotopes as tracers of marine particle fluxes and deep water circulation in the Mediterranean Sea[END_REF].

The total 232 Th chemical blank was typically 42 ± 25 pg, although few blanks as high as 1400 pg were observed. Given the good agreement between the particulate 232 Th concentrations obtained with the Niskin filtrations (see Fig. ES1) and the ISP 232 Th concentrations corrected with the low blank values, the highest blanks were not taken into account to correct the ISP data. Total 230 Th chemical blank was typically 1 ± 1 fg. As for 232 Th, a small number of very large blanks (up to 125 fg) were obtained but assumed to be outliers.

Mass spectrometry

For dissolved and in situ pump samples, Th isotopes were measured on a MC-ICPMS Neptune + (Thermo-Fischer). The dissolved samples were analyzed with the initial configuration of the Neptune + with an Apex desolvator. The particulate samples were analyzed after the Neptune's upgrade with a Jet interface and an Aridus 2 desolvating nebulizer. The acquisition conditions are given in [START_REF] Gdaniec | Thorium and protactinium isotopes as tracers of marine particle fluxes and deep water circulation in the Mediterranean Sea[END_REF]. The accuracy of the MC-ICPMS measurement was checked by analysis of Th standard IRMM035 ( 230 Th/ 232 Th = (11.31±0.12)×10 -6 , n=6) and IRMM036 (3.059±0.033)×10 -6 , n=8), in good agreement with reference values: 230 Th/ 232 Th = (11.38±0.10)×10 -6 for IRMM035 and 230 Th/ 232 Th = (3.047±0.024)×10 -6 for IRMM036 [START_REF] Sims | An inter-laboratory assessment of the thorium isotopic composition of synthetic and rock reference materials[END_REF]. The quality of the analyses was checked through the GEOTRACES intercalibration exercise where we obtained precise and accurate dissolved Th concentrations well within the consensus values [START_REF] Anderson | GEOTRACES intercalibration of 230 Th, 232 Th, 231 Pa, and prospects for 10 Be[END_REF].

For Niskin filtration samples, 232 Th was analyzed by quadrupolar ICPMS (Agilent 7500CE of the Observatoire Midi Pyrénées, Toulouse). Indium and rhenium (100 ppt, added to the samples) were used as internal standards to monitor instrument sensitivity variations.

Reproducibility calculated from 3 duplicated samples (2 different Niskins) was on average 12%

(2 standard deviations). Detection limit was 0.5 pg/kg. Along the section, the lowest 232 Thd are found at station S3 and S4, far from the coasts of Africa and Antarctica ( 232 Thd = 4.7-30 pg/kg). All the 232 Thd duplicates agree within analytical uncertainties. for a given depth, 230 Thd is higher at stations S4 and S5 than at S1, S2 and S3. All 230 Thd duplicates agree within 9% (3 of them agreeing within the analytical uncertainties).

Results

3

Following GEOTRACES good practices, a comparison with 2 stations of the Zero-Drake cruise (Rutgers van de Loeff et al., 2016) is provided (Fig. ES3). Profiles agree within uncertainties for 230 Th and more qualitatively for 232 Th.

Particulate thorium

Particulate 232 Th concentrations ( 232 Thp) have been already reported elsewhere [START_REF] Garcia-Solsona | Rare earth elements and Nd isotopes tracing water mass mixing and particleseawater interactions in the SE Atlantic[END_REF]. 232 Thp from ISP filtrations range from 0.08 pg/kg to 371 pg/kg (Fig. 4a, Tab. ES2).

There is generally a good agreement between 232 Thp profiles obtained with in situ pumps and profiles obtained with Niskin filtrations. The most notable exception is above the seafloor at station S3, where Niskin bottle filtration gives 232 Thp = 27 pg/kg, whereas ISP filtration gives 232 Thp = 2.5 pg/kg. Apart from this particular feature recorded at S3, the 232 Thp profiles generally increase with depth with a sharp increase above the seafloor. Along the section, the highest 232 Thp concentrations are found at station S1 and S2, which are the stations closest to the African coast.

All the 232 Thp duplicates agree within analytical uncertainties. The variation in 232 Thp observed between in situ pumps and Niskin filtrations (Fig. ES1) is of the same order than observed for other insoluble elements [START_REF] Planquette | Sampling for particulate trace element determination using water sampling bottles: methodology and comparison to in situ pumps[END_REF].

Particulate 230 Th concentrations ( 230 Thp) range from 0.02 fg/kg to 28 fg/kg (Fig. 4b). From the surface down to 3000 m depth, the 230 Th profiles of S2 and S3 are similar to S1. In the deepest waters, S2 has a strong increase of its dissolved and particulate 230 Th (although less than S1) also suggesting inputs from the sediments or lateral transport from the African margin. By contrast, at S3, the low 230 Thd concentrations at a depth of 4500 m could also be attributed to bottom scavenging (a low REE content on the same sample was also noted by [START_REF] Garcia-Solsona | Rare earth elements and Nd isotopes tracing water mass mixing and particleseawater interactions in the SE Atlantic[END_REF].

The strongly concave shapes of the 230 Thd profiles at S4 and S5 strongly suggests these profiles are impacted by the upwelling of LCDW (Rutgers van der Loeff and Berger, 1993).

The 230 Thd profile at S3 has also a concave shape that can be regarded as a result of the LCDW upwelling. However, above 3000 m depth, it is not so different from the S1 profile where no LCDW upwelling is suspected. Also, the low 230 Thd concentrations at a depth of 4500 m, that accentuate the nonlinearity of the profile, could also be attributed to bottom scavenging.

At stations S1, S2 and S3, the 230 Thp profiles have very similar slopes, particularly from the surface down to 1500 m depth. We consider the conservation equation of the 1D reversible scavenging model at steady state (and neglecting 230 Th radioactive decay):

𝜕𝐶 𝑡 𝜕𝑡 = 𝑃 𝑑 -𝑤 𝑝 𝜕𝐶 𝑝 𝜕𝑧 = 0 (2)
Where Ct is the total (dissolved + particulate) 230 Th concentration, Pd is the in situ production rate of 230 Th (Pd = 0.56 fg/kg/y), wp is the particle settling speed and Cp the particulate 230 Th concentration. Hence, wp is given by: wp = P/(𝜕 230 Thp/𝜕z). Using an average slope of 0.001 fg/kg/m, we obtain that wp = 560 m/d, in agreement with Rutgers van der Loeff and Berger (1993).

In the southern part of the section, we first use the mixing-scavenging model introduced by Rutgers van der Loeff and Berger (1993) to evaluate the impact of lower Circumpolar Deep Water upwelling and input into the Weddell Gyre. This "scavengingmixing" model relies on several assumptions:

(i) The representation of mixing with a constant restoring term at all depths (eq. 3). This is largely ad hoc, so that the physical interpretation of Cit and w is unclear.

(ii) the omission of the effects of vertical transport by advection and diffusion, (iii) steady state, (iv) the partition coefficient K is vertically uniform. It implies that the 230 Thxs particulate fraction is constant over the depth profile. Clear variations of the 230 Thxs particulate fraction (Tab. ES2) occur in the deepest samples either due to particle inputs from the margins (stations S1, S2 and S5, see [START_REF] Gdaniec | 231 Pa and 230 Th in the Arctic Ocean: implications for Boundary Scavenging and 231 Pa-230 Th fractionation in the Eurasian Basin[END_REF] or following hydrothermal scavenging (station S3).

(v) K is very small compared to 1, so that Cd could be approximated by Ct. Apart for the deepest samples at stations S1 and S2, it introduces an underestimation of wp by 13 % at most.

(vi) the particle settling velocity wp is vertically uniform. This assumption is deduced the from linear dissolved and particulate 230 Th profiles observed in regions unaffected by strong deep currents (Roy- [START_REF] Roy-Barman | The sources and the fates of thorium[END_REF].

(vii) Cd is negligible near the surface.

Using the above assumptions, the tracer conservation equation of this model is:

𝜕𝐶 𝑡 𝜕𝑡 = 𝑃 𝑑 + (𝐶 𝑖_𝑡 -𝐶 𝑡 ) 𝜏 𝑤 -𝑤 𝑝 𝜕𝐶 𝑝 𝜕𝑧 = 0 (3)
Where w is restoring time scale taken equivalent to the water residence time and Ci_t is the total 230 Th concentration of the inflowing water. Introducing a partition coefficient between dissolved and particulate Th, K= Cp /Cd and integrating equation 3 vertically over the water column, the 230 Thd concentration at any depth is given by:

𝐶 𝑑 = (𝐶 𝑖_𝑡 + 𝑃 𝑑 𝜏 𝑤 ) × (1 -exp (- 𝑧 𝑤 𝑝 𝐾𝜏 𝑠 )) (4)
For station 5, we use the 230 Thi = 24 fg/kg for the LCDW (Rutgers van der Loeff and Berger, 1993), K=0.04 is the average ratio measured at this station. Then, we adjust wp = 875 m/y and w= 20 y to obtain a good agreement between the modelled curve and the data (Fig. 5). The low w obtained here compared to 45 y used by Rutgers van der [START_REF] Van Der Loeff | Scavenging of 230 Th and 231 Pa near the antarctic polar front in the South Atlantic[END_REF] derives presumably from the relatively peripheral position of station S5 in the Weddell Sea.

We confirm that the scavenging residence time in the Weddell Sea (s=d/(2Kwp) = 55 y for a water column of a depth of d = 4000 m) is larger than in the subtropical gyre (w = 25 y for a water column of a depth of 4000 m). On the whole, the model explains why the 230 Thd concentrations are high and do not increase linearly with depth in the Weddell Gyre.

We apply the same model to station S4. Here 230 Thi is also adjusted to obtain a good data-model agreement. We note that 2 modeled curves based on relatively different sets of w andwp values (ACC high and ACC low, Fig. 5) agree with the S4 data. This illustrates a limit of this method to estimate particle settling velocities. This advection-scavenging model has been largely used to describe the vertical distribution of 230 Th in the Southern Ocean [START_REF] Coppola | Thorium isotopes as tracers of particles dynamics in the Indian sector of the Southern Ocean[END_REF], Venchiarutti et al., 2011[START_REF] Van Der Loeff | Meridional circulation across the Antarctic Circumpolar Current serves as a double 231 Pa and 230 Th trap[END_REF]. Despite crude assumptions (e.g., uniform influence throughout the water column by a single water mass), this mixing-scavenging model has the great merit of highlighting the impact of upwelling on 230 Thd profiles and to provide an estimate of the particle settling speed when the 1D production-scavenging model cannot be used. Nevertheless, this model cannot fully describe the upwelling and mixing of 2 different water masses along an isopycnal surface as observed during BGH and described in the next section. Salinity data help to illustrate this point. For a conservative tracer (Pd = 0, wp = 0 and Cp = 0) and assuming steady state, equation (3) reduces to:

(𝑆 𝑡_𝑖 -𝑆 𝑡 ) = 0 (5)
This expresses that salinity is constant and equal to the salinity of the inflowing water (St = St_i).

In this model, LCDW is advected at all levels in the Weddell Gyre. It implies that the Weddell Gyre must be filled only with LCDW. Nevertheless, the 230 Tht concentration in the Weddell Tht concentration gradients are small on a given neutral density surface (Fig. 6). Figure 6 highlights that distribution of 230 Thd is mainly driven by transport from upwelling and mixing along isopycnal surfaces. There is more variability of 232 Thd or 232 Tht (dissolved + particulate 232 Th) on a given neutral density value, presumably because large lithogenic inputs occur at stations S1 and S2 (Fig. 6). In the following, we discuss concentration changes along isopycnal surfaces. As these changes are small, we restrict this discussion to the Bonus GoodHope data, all generated in the same laboratory, in order to avoid the systematic inter-laboratory differences identified during the GEOTRACES intercalibration [START_REF] Anderson | GEOTRACES intercalibration of 230 Th, 232 Th, 231 Pa, and prospects for 10 Be[END_REF].

We focus on 2 neutral density surfaces:  n = 27.865 kg/m 3 and  n = 28.094 kg/m 3 (Fig ES2). The neutral density surface  n = 27.865 kg/m 3 can be followed from z=1500 m at station S1 to z=139 m at station S5. It crosses the following water masses: I-UCDW (S1), A-UCDW (S2 and S3), DP-UCDW (S4) and AASW/WW (S5). The neutral density surface:  n = 28.094 kg/m 3 can be followed from z = 3000 m at station S1 to z = 200 m at station S5. It is close to the neutral density of the deep salinity maximum observed in 2004 along the same section [START_REF] Gladyshev | A hydrographic section from South Africa to the southern limit of the Antarctic Circumpolar Current at the Greenwich meridian[END_REF]. It crosses the following water masses: SE-NADW (S1 and S2), SW-NADW/AABW (S3), LCDW (S4) and AASW/WW (S5).

We choose these two neutral density surfaces because: (1) they are found along the whole section, (2) samples along these surfaces tend to form a linear trend on a -S diagram, indicating a binary mixing (Fig. 7), and (3) they are not interrupted by the Agulhas Ridge. When, at a given station, no sample was collected exactly on the selected neutral density surface, concentrations were linearly interpolated between the samples collected just above and below this surface at the same station. When the particulate Th was not analyzed with in situ pump samples at the target depth or in its vicinity (e.g., intercal station), we used the 232 Thp concentration obtained with the Niskin filtration profiles and the 230 Thp was calculated with the ( 230 Thp/ 232 Thp) ratio linearly interpolated between the ratios measured on particles just above and below or at the nearby station:

230 Thp-estimates = 232 Thp-Niskin filtration * ( 230 Thp/ 232 Thp)interpolated ( 6)

 n =27.865 kg/m 3
On a -S diagram, S1, S2, S3 and S4 plot along a straight line, suggesting a binary mixing between S1 and S4 (Fig. 7a). S5 does not fall on this line, possibly due to ocean-atmosphere interactions [START_REF] Abadie | Iron isotopes reveal distinct dissolved iron sources and pathways in the intermediate versus deep Southern Ocean[END_REF], so that this station does not correspond to the same binary mixing along this neutral density surface. Therefore, in the following, we will apply the mixing model mostly between stations S1 and S4 for this isopycnal.

On a S-232 Tht diagram, the stations S2 and S3 stand below the conservative mixing line between S1 and S4 (Fig. 7c). It corresponds to a depletion (net removal of 232 Th on sinking particles) of 232 Tht compared to a conservative mixing between S1 and S4:  232 Tht<0. The strong 232 Tht inputs from the South African margin to the ACC are marked by a strong 232 Tht concentration gradient ( 232 Tht decrease by a factor ~4 between stations S1 and S4). A small 232 Tht gradient between Station S5 and S4 suggests smaller inputs from the Weddell Sea.

On a S-230 Tht diagrams, there is much less change in 230 Tht concentration (~ 30% between stations S1 and S4) compared to 232 Tht presumably because the margin is not a strong source of 230 Tht for the ACC. Nevertheless, the 230 Tht concentrations generally decrease southward.

However, stations S2 and S3 stand above the conservative mixing line between S1 and S4 (Fig. 7e). This trend would correspond to an excess of 230 Tht compared to a conservative mixing between S1 and S4:  230 Tht>0. There is an accumulation of 230 Th due to in situ production despite the net removal of 232 Tht during the same mixing process.

 n =28.094 kg/m 3
On a -S diagram, S3, S4, intercal and S5 plot along a straight line, suggesting a binary mixing between S3 and S5 (Fig. 7b). S1 and S2 do not fall on this line, so that these 2 stations do not correspond to the same binary mixing along this neutral density surface. Therefore, in the following, we will restrict our analysis to stations S3 to S5 and, in the next section, we will apply the mixing model only between stations S3 and S5.

On a S-232 Tht diagram, we observe the same trend as on the previous isopycnal: a strong 232 Tht decrease from S1 to S3 and a weaker decrease from S5 to S4 (Fig 7d). This illustrates the role of the margins (and particularly the South African margin) as sources of 232 Th for the ACC.

We argued in the previous paragraph that waters from stations S1 and S2 do not seem directly involved in the mixing observed further south. However, high 232 Tht concentrations are also observed in the SW-NADW sampled upstream of the BGH section. These high values are acquired as this water flows along the Namibian margin [START_REF] Scholten | Advection and scavenging: Effects on 230Th and 231Pa distribution off Southwest Africa[END_REF] so that it could be a suitable 232 Tht source fueling the southern part of the BGH section. S4 and intercal stations stand below the S3-S5 binary mixing line, reflecting 232 Tht depletion at these 2 stations compared to conservative mixing, expressed as  232 Tht<0 (Fig 7). More generally, paired 232 Tht-S data indicates that there is a net Th removal by the settling particles while water masses mix between stations S3 and S5.

As for  n =27.865 kg/m 3 , the relative variations of 230 Tht are small (Fig. 7f). On a S-230 Tht diagram, the S4 and intercal stations stand above the S3-S5 binary conservative mixing trend (S1 and S2 are not considered because they are not involved in the mixing on this isopycnal). This trend would correspond to an excess of 230 Tht compared to a conservative mixing:  230 Tht>0. Like on the previous isopycnal, 230 Th accumulates due to in situ production despite the net removal of 232 Tht during water mass mixing. See table ES3 for the data calculated on isopycnals.

4.3. Modeling advection-diffusion-scavenging along a neutral density surface.

In the deep ocean, water mass mixing is thought to occur preferentially along isopycnal surfaces. There is both advection and eddy diffusion along the steep isopycnal surfaces of the Southern Ocean [START_REF] Zika | Diagnosing the Southern Ocean overturning from tracer fields[END_REF][START_REF] Garabato | Short-circuiting of the overturning circulation in the Antarctic Circumpolar Current[END_REF]. While diapycnal mixing is important in the Drake Passage due to the rough topography, it appears to be much weaker out of the passage [START_REF] Watson | Rapid cross-density ocean mixing at mid depths in Drake Passage measured by tracer release[END_REF], so we will neglect it for the present analysis. We approximate the combined effect of advection and diffusion along an isopycnal by extending the conservation equation ( 2) with the explicit introduction of advection and diffusion terms in an infinitesimal volume of water at steady state [START_REF] Jenkins | 3 H and 3 He in the Beta Triangle: Observations of gyre ventilation and oxygen utilization rates[END_REF]:

𝐾 𝑖 𝜕 2 𝐶 𝑡 𝜕𝑥 2 -𝑢 𝑖 𝜕𝐶 𝑡 𝜕𝑥 + 𝑃 𝑑 -𝑤 𝑝 ( 𝜕𝐶 𝑝 𝜕𝑧 ) = 0 (7)
Where Ct is the total (dissolved + particulate) concentration of the tracer, Cp the particulate concentration of the tracer, Ki is the cross-stream isopycnal eddy diffusion coefficient and ui is the cross-stream isopycnal advection velocity. The x-axis is parallel to the isopycnal surface in the meridional direction and oriented northward. The z axis is oriented perpendicular to the isopycnal surface and hence is almost vertical. Just like in the advection-scavenging model (equ. 3), keeping wp out of the partial derivative with respect to z (writing 𝑤 𝑝 ( 𝜕𝐶 𝑝

𝜕𝑧

) instead of ( 𝜕(𝑤 𝑝 𝐶 𝑝 )

𝜕𝑧

)) requires to assume that wp does not vary with z. We integrate Eq. 7 along isopycnal surfaces of the BGH section between the stations corresponding to the end members of water mass mixing (stations for S1 and S4 for  n =27.865 kg/m 3 and S3 and S5 for  n =28.094 kg/m 3 ). Unlike most Th models that are built by vertically integrating the conservation equation, here we integrate the conservation equation along surfaces that are only slightly tilted relative to the horizontal. This is a simplified view of water transport because while water moves across the ACC, it is also rapidly advected eastward by the ACC (S2-S4) and the northern limb of the Weddell Gyre (S5). We assume that Ki, ui, wp and (dCp/dz) are all constant along x in order to obtain analytical solutions. We note Ct = Ct-conc + Ct, where Ct-conc is the concentration of the tracer if it had a conservative behavior and Ct ))

(
Here A and B represent stations S1 and S4 or S3 and S5 along the neutral density surface. For a conservative tracer such as salinity, the concentration is given by:

𝑆 = 𝑒 𝑢 𝑖 𝐾 𝑖 (𝑥-𝑥 𝐴 ) -1 𝑒 𝑢 𝑖 𝐾 𝑖 (𝑥 𝐵 -𝑥 𝐴 ) -1 (𝑆 𝐵 -𝑆 𝐴 ) + 𝑆 𝐴 (10)
Figure 8: Advection-diffusion modelling of salinity. The distance considered is the meridian distance from the equator (d = 0 at the equator). Red curves are obtained with equation 10. ui/Ki is adjusted by least squares just based on the stations corresponding to the 2-water mass mixing and by forcing the curves to pass through the 2 water mass end members. The values of ui/Ki are given in Table 1. Dotted curves are extrapolated to the rest of the section but the corresponding stations are not used to calculate the curve parameters. Equation 10 allows determining the ui/Ki ratio using salinity data along the BGH 616 section. For  n = 27.865 kg/m 3 , we obtain ui/Ki = -1.57×10 -6 m -1 (Fig. 8, Tab. 1). The minus 617 sign corresponds to a southward advection (ui <0 and Ki is always positive). This is consistent 618 with the southward upwelling expected based on the average circulation of UCDW in the 619 Southern Ocean within this density range (section 2.1). This feature is observed for the 2 620 isopycnals studied here (see below for  n = 28.094 kg/m 3 ). 621 below the depth of the neutral density surface at each station (Fig. 9). For  n = 27.865 kg/m 3 , d 232 Thp/dz varies strongly between S1 and S4. Interestingly, we obtain low d 232 Thp/dz values when the 232 Tht horizontal gradient is low (e.g., d 232 Thp/dz = -0.001 pg/kg/m at S3) and higher d 232 Thp/dz values when the 232 Tht horizontal gradient is high (e.g., d 232 Thp/dz = + 0.006 pg/kg/m at S1). In other words, the increase of 232 Tht with depth is due to the horizontal transport of 232 Th by eddy diffusion from the margin, whereas without inputs from the margins the dissolved and particulate 232 Th concentration would be expected to be constant with depth (Roy- [START_REF] Roy-Barman | The sources and the fates of thorium[END_REF]. Ideally, it would be best to integrate equation 7 with d 232 Thp/dz and d 230 Thp/dz varying with x. In practice, it is difficult to determine the evolution of these gradients along the isopycnal or even to determine the exact average gradients between stations S1 and S4. As a first approach, we will just use the range of the gradients measured between S1 and S4 as input values to determine wp, ui and Ki. For 232 Th, we use a range d 232 Thp/dz from -0.001 to + 0.012 pg/kg/m (Fig. 9, Tab. 1). For 230 Th, most 230 Thp data plot along a slope ranging between 3.8×10 -4 fg/kg/m and 6.2 ×10 -4 fg/kg/m (Fig. 9, Tab. 1).

We then estimate wp and ui with a least square approach. Using Monte Carlo simulations, we maximize the coefficient of determination (R 2 ) of the 230 Tht versus S and The 2 Ki values estimated here are on the upper end of the range (300 to 1800 m 2 /s) reported in the Southern Ocean, whereas ui is on the lower but of the same order of magnitude than an estimate of the cross-stream horizontal velocity associated with the upwelling in the ACC (~ -0.010 ± 0.003 m/s, Naveira [START_REF] Garabato | Short-circuiting of the overturning circulation in the Antarctic Circumpolar Current[END_REF][START_REF] Zika | Diagnosing the Southern Ocean overturning from tracer fields[END_REF].

We have estimated wp at relatively close depths on the 2 isopycnals (1200-1500 m for  n = 27.865 kg/m 3 at stations S3 and S2 and 1000-1200 m for  n = 28.094 kg/m 3 at stations intercal and S4, Fig. ES2). Taking the corresponding wp at face value (but keeping in mind that they strongly rely on the particulate concentration gradient estimate), it seems that there is a southward decrease of wp (for  n = 27.865, wp = 600 m/y and for  n = 28.094 kg/m 3 , wp = 300 m/y). This southward decrease of settling speed seems consistent with the 238 U-234 Th disequilibrium in the surface waters that decreased from the north to the south of the polar front zone (PFZ) during the BGH cruise [START_REF] Planchon | Late summer particulate organic carbon export and twilight zone remineralisation in the Atlantic sector of the Southern Ocean[END_REF] and 2 weeks earlier on the same transect during the ANTXXIV/3 cruise (Rutgers van der Loeff et al., 2011). This gradient of the 238 U-234 Th disequilibrium indicates a higher removal rate of 234 Th and suggests a higher particle settling velocity in the north compared the south of the PFZ.

It is rewarding that the particle settling speeds obtained with the isopycnal mixing model Consequently, the wp does not depend on sharp changes of (d 230 Thp/dz) that could result from random errors on 230 Thp at intermediate depths. In the 1D mixing-scavenging model some advection is introduced but the model is still integrated from the surface to the seafloor, so that random errors on 230 Thp tend to be smoothed out (Rutgers van der Loeff and Berger, 1993). In addition, 234 U in situ decay is not the only source of dissolved 230 Th, because advection also brings dissolved and particulate 230 Th, which adds up uncertainty on the input rate of 230 Th.

In the isopycnal mixing model presented here, the uncertainties cited above can have much larger effects on wp because i) (d 230 Thp/dz) is calculated over much thinner depth ranges and ii) combined effects of high water flow rates and small 230 Th concentration gradients implies that the net effect of circulation on 230 Th transport is subject to large relative uncertainties.

Nevertheless, the isopycnal mixing model opens the possibility to calculate the vertical evolution of wp in highly dynamical environments where the hypotheses of the 1D model are clearly not fulfilled.

In a recent series of paper, non-linear dissolved 230 Th profiles are modelled with a purely 1D approach by adjusting adsorption/desorption/remineralization rates and particle settling velocities that are assumed to vary with depth [START_REF] Lerner | Testing models of thorium and particle cycling in the ocean using data from station GT11-22 of the US GEOTRACES North Atlantic Section[END_REF][START_REF] Lerner | Kinetics of thorium and particle cycling along the US GEOTRACES North Atlantic Transect[END_REF][START_REF] Lerner | Effects of particle composition on thorium scavenging in the North Atlantic[END_REF].

These numerous parameters are adjusted by using different Th isotopes ( 228 Th, 230 Th, 234 Th).

This approach is different but not inconsistent with our approach because: (1) Unlike the present paper, Lerner et al focus on areas where advection and eddy diffusion are expected to be negligible;

(2) In the present paper, we make no hypothesis on the processes driving dissolved-particulate transfers (Eq. 7). Ultimately, it should be possible to determine adsorption/desorption rates in areas with weak currents and then to use these values in regions with a vigorous circulation such as the ACC. Lerner's approach treats explicitly the case of wp varying with depth, whereas our approach to leave the possibility to determine the particle settling speed at different depths/density. To do so rigorously, it will be necessary to treat wp like a parameter varying with depth (see section 4.3).

Conclusion

While we still used crude approximations, the model proposed here opens the possibility to determine the particle settling speed and to bring constraints on the isopycnal mixing rate at different vertical levels and with a more realistic description of the vertical structure of the water column. The improvement that will come out of the method proposed here will require sampling with higher horizontal and vertical sampling resolutions, a goal partly achieved in the framework of the GEOTRACES program. It also requires progresses for the particulate Th analysis and constant efforts for the intercalibration of the different cruises.

Fig. 1 .Fig. 2 :

 12 Fig. 1. Sampling stations during the BONUS-GoodHope section and main fronts (following Sokolov and Rintoul, 2009; mean path position): STF = SubTropical Front; SAF = SubAntarctic Front (middle branch); PF = Polar Front (middle branch); N-SACCF = northern branch of the Southern ACC front; N-SB = northern branch of the Southern Boundary. A.R.: Aghulas Ridge.

  .1. Dissolved thorium Dissolved 232 Th concentrations ( 232 Thd) range from 4.7 pg/kg to 288 pg/kg (Fig. 2a, Tab. ES 1) and generally increase with depth. At stations S1 and S2, the 232 Thd profiles have a local maximum around 200-400 m depth and a sharp increase at bottom depths. Conversely, at station S3, 232 Thd for the deepest sample is relatively low. This feature was also noted for dissolved REE (Garcia-Solsona et al. 2014) and for 230 Thd (see below) analyzed on the same water sample.

Fig. 3 .

 3 Fig. 3. Dissolved thorium profiles. a: dissolved 232 Th. b: dissolved 230 Th.

Fig. 4 .

 4 Fig. 4. Particulate Thorium profiles collected with in situ pumps. (a) Particulate 232 Th. (b) Particulate 230 Th. Insets: expanded concentration scales to display particulate concentration maximum at S1.

Fig. 5 .

 5 Fig. 5. Modelling of the Bonus Good Hope profiles. The dotted and dashed curves are calculated with Eq. 4. Weddell: Ci = 24 fg/kg, wp = 875 m/y, w = 20 y and K = 0.04. ACC low: Ci = 12 fg/kg, wp = 803 m/y, w = 12 y and K = 0.07. ACC high: Ci = 15 fg/kg, wp = 1785 m/y, w = 4 y and K = 0.125. linear: Ci = 0 fg/kg, wp = 800 m/y, w = + ∞ and K = 0.1.

  Gyre changes with depth due to the scavenging transport of230 Tht by sinking particles. Hence, this mixing-scavenging model cannot represent the mixing of water masses with different origins and labeled with different  and S signatures as they are observed along the BGH section and elsewhere in the Southern Ocean. In principle, general circulation models, associated or coupled with biogeochemical model, have the potential to reproduce the coupled effect of scavenging and mixing, but the present generation of coupled models do not accurately describe the vertical distribution of scavenging intensity (e.g.:[START_REF] Dutay | The influence of particle size and type on 231 Pa and 230 Th simulation with a global coupled biogeochemical-ocean general circulation model: a first approach[END_REF][START_REF] Gu | 231 Pa and 230 Th in the ocean model of the Community Earth System Model (CESM1.3)[END_REF][START_REF] Van Hulten | A global scavenging and circulation ocean model of thorium-230 and protactinium-231 with an improved particle dynamics (NEMO-ProThorP 0.1)[END_REF]. 4.2. Th isotopes along neutral density surfaces While there are large horizontal concentration gradients of 230 Thd and total (dissolved + ¨particulate) 230 Th ( 230 Tht) along the BGH section (mostly between station S3 to S5, Fig 3), 230

Fig. 6 .

 6 Fig. 6. Neutral density versus 232 Tht and 230 Tht concentration. Grey dots: S1 and S2 stations with high 232 Tht contents. White dots: S3, S4 and S5. Dashed line:  n = 27.865 kg/m 3 . Dotted line:  n = 28.094 kg/m 3 .

Figure 7 :

 7 Figure 7: Potential temperature, 232 Tht and 230 Tht versus salinity on different neutral density surfaces. Left:  n =27.865 kg/m 3 . Right:  n =28.094 kg/m 3 . a-b: Potential temperature versus salinity. c-d: 232 Tht

Figure 9 :

 9 Figure 9: 232 Thp and 230 Thp profiles. For each isopycnal, we also plot the sample analyzed just above

Figure 10 :

 10 Figure 10: Advection-diffusion modelling of S and Tht isotopes compared to data. Staked red curves: modelled curves with the highest R 2 values (Tab. 1). Thin black dotted line: conservative mixing between S3 and S5.

232

  Tht versus S curves simultaneously. In equation 9, wp and the dThp/dz terms are present as products so that their opposite variations cancel each other. As a consequence, a relatively large range of wp, d 230 Thp/dz and d 232 Thp/dz values yield similar R 2 values. Therefore, we report wp and ui obtained for a range of high R 2 values. We obtain wp = 674 ± 250 m/y and ui = -0.0033 ± 0.0014 m/s. Using the ui/Ki ratio estimated previously, we deduce Ki = 2000 ± 840 m/s 2 .   n = 28.094 kg/m 3For  n = 28.094 kg/m 3 , most of the 230 Thp data plot along a slope of 5×10 -4 fg/kg/m, whereas232 Thp data plot along a slope of 2×10 -4 pg/kg/m (Fig.10). In both diagrams, the slope of the S4 station is poorly defined and potentially lower than the average. Using the same method, we obtain ui/Ki = -1.65×10 -6 m -1 , wp = 435 ± 324 m/y, ui = -0.0036 ± 0.0008 m/s and Ki = 2180 ± 480 m/s 2 (Fig 8-9, Tab. 1).

(

  which does not integrate the conservation equation over the water column) are in the range of the speeds generally obtained with the vertical 1D model (see section 4.1; Rutgers van der[START_REF] Van Der Loeff | Scavenging of 230 Th and 231 Pa near the antarctic polar front in the South Atlantic[END_REF]. In the 1D production-scavenging model, wp is given by (see section 4.1.): wp = P/(d 230 Thp/dz). However, wp is usually obtained by averaging (d 230 Thp/dz) over the whole water column or a large part of the water column as we did in section 4.1.

Fig. ES1 :Figure ES3 :

 ES1ES3 Fig. ES1: comparison of particulate 232 Th collected by in situ pump and Niskin filtrations. 947

Table ES1 :

 ES1 Dissolved Thorium data. .808 27.7752 40.2 ± 1.4 11.64 ± 0.18 11.46 ± 0.18 SE-NADW 2.132 34.838 27.8316 42.0 ± 1.3 14.71 ± 0.19 14.52 ±0.19 SE-NADW 1.945 34.836 27.8444 42.2 ± 1.4 14.94 ± 0.17 14.76 ± 0.17 SE-NADW 0.892 34.749 27.8506 70.0 ± 1.4 17.55 ± 0.23 17.25 ± 0.23 AABW 0.716 34.733 27.8493 69.1 ± 1.6 21.70 ± 0.27 21.40 ±0.27 AABW 0.568 34.72 27.8481 288.5 ± 3.5 37.10 ± 0.41 35.84 ± 0.41 AABW S2: 42.47°S, 08.93°E, m bottom depth 20 12.951 34.521 26.0547 11.1 ± 1.4 0.74 ± 0.02 0.70 ± 0.03 ICW 124 9.083 34.482 26.7815 19.8 ± 1.4 1.84 ± 0.03 1.75 ± 0.04 ICW 292 6.654 34.267 27.0053 36.2 ± 1.3 3.04 ± 0.04 2.89 ± 0.04 ICW/a-AAIW 398 6.01 34.292 27.1231 38.8 ± 1.5 4.13 ± 0.06 3.97 ± 0.06 ICW/a-AAIW 595 4.138 34.185 27.2829 30.9 ± 1.4 5.14 ± 0.08 5.01 ± 0.08 a-AAIW 940 3.075 34.313 27.5162 25.7 ± 1.3 7.02 ± 0.09 6.90 ± 0.09 a-AAIW 2.626 34.585 27.7945 26.1 ± 1.3 9.19 ± 0.11 9.07 ± 0.11 UCDW 2.346 34.79 27.994 31.1 ± 1.3 13.36 ± 0.16 13.22 ± 0.16 SE-NADW 1.967 34.827 28.0789 35.8 ± 1.3 14.43 ± 0.16 14.28 ± 0.16 SE-NADW 1.427 34.792 28.1359 36.9 ± 1.6 16.13 ± 0.24 15.96 ± 0.24 SE-NADW 0.81 34.741 28.1961 63.9 ± 1.5 17.96 ± 0.22 17.68 ± 0.22 AABW .407 27.6574 16.8 ± 1.3 6.46 ± 0.11 6.39 ± 0.11 a-AAIW/UCDW 2.447 34.575 27.8099 18.2 ± 1.3 8.26 ± 0.10 8.18 ± 0.10 A_UCDW 2.334 34.716 27.936 22.2 ± 1.5 10.93 ± 0.17 10.83 ± 0.17 SW-NADW 1.971 34.771 28.0376 23.2 ± 1.2 12.68 ± 0.15 12.58 ± 0.15 SW-NADW 23.2 ± 1.5 13.45 ± 0.16 13.34 ± 0.16 0.955 34.731 28.1653 29.2 ± 1.6 14.94 ± 0.19 14.81 ± 0.19 SW-NADW 0.415 34.695 28.2262 18.5 ± 1.3 8.83 ± 0.12 8.75 ± 0.12 AABW Table ES3: Total thorium isotopes recalculated on neutral density surfaces

	918
	Station S1: 36.50°S, 13.10°E, m bottom depth Depth (m)  (°C) 29 20.785 35.617 25.0222 34.7 ± 1.3 1.21 ± 0.03 1.06 ± 0.03 ICW Salinity  n (kg/m 3 ) 232 Thd 1 (pg/kg) 230 Thp 1 (pg/kg) 230 Thp-xs 2 (pg/kg) 198 12.364 35.014 26.53 52.1 ± 1.3 2.48 ± 0.04 2.25 ± 0.04 ICW 397 9.824 34.757 26.7918 50.1 ± 1.5 4.52 ± 0.06 4.30 ± 0.06 ICW 753 4.973 34.346 27.1602 37.2 ± 1.3 5.45 ± 0.08 5.29 ± 0.08 i-AAIW Water mass 3.279 34.551 27.5004 38.0 ± 1.3 7.38 ± 0.12 7.22 ± 0.12 i-AAIW 2.515 34S3: 47.55°S, 04.37°E, m bottom depth 40 6.32 33.733 26.5685 7.5 ± 1.4 1.19 ± 0.04 1.15 ± 0.04 AASW 124 4.029 33.867 27.0072 8.3 ± 1.3 1.82 ± 0.04 1.79 ± 0.04 AASW 248 3.633 34.101 27.2697 15.3 ± 1.4 3.55 ± 0.07 3.48 ± 0.07 AASW/a-AAIW 495 2.739 34.233 27.4888 17.1 ± 1.2 5.07 ± 0.07 5.00 ± 0.07 a-AAIW 742 (°C) Salinity 232 Th (pg/kg) error 230 Th (fg/kg) error  n = 27.865 kg/m 3 S1 -1550 2.732 34.695 10.67 0.15 72.4 1.5 S2 -1631 2.564 34.667 11.51 0.12 36.3 1.3 S3 -1227 2.44 34.644 10.07 0.13 25.1 2.9 S4 -345 1.838 34.563 7.59 0.11 16.8 1.8 S5 -98 -0.665 34.348 8.52 0.11 21.8 0.1  n = 28.094 kg/m 3 S1 3073 1.86 34.827 15.1 0.2 60.7 1.5 S2 3082 1.849 34.823 17.3 0.8 52.9 0.2 S3 2440 1.722 34.756 15.2 0.2 32.1 1.3 S4 1205 1.322 34.726 15.9 0.2 29.3 1.3 intercal -999.6 1.205 34.71 16.7 0.2 30.1 1.2 intercal -999.6 1.205 34.71 16.7 0.2 30.2 1.3 2.591 34Depth (m)  S5 247 0.28616 34.619 16.1 0.2 31.9 1.6
	919
	920
	921
	922
	923
	924
	925
	926

Th d (pg/kg) Depth (m) (a)

Acknowledgements:

The authors would like to thank S. Speich and M. Boyé, the chief scientists of the BGH cruise, and the crew of the research vessel Marion Dufresne II (IPEV). We acknowledge Christophe Guillerm, Lionel Scouarnec, Thomas Arsouze and Amandine Radic for the sample collection on board. We are grateful to the 2 anonymous reviewers who improved the manuscript with their very helpful comments. Support provided by the ANR (Programme Blanc 2007: ANR-07-BLAN-0146-02) and the INSU/CNRS (programme INSU/LEFE-CYBER) is gratefully acknowledged.

, where the ( 230 Th/ 232 Th)Litogenic = 4.0 × 10 -6 (mol/mol), close to average crust 230 Th/ 232 Th = 4.1 x 10 -6 ).

The lithogenic correction is small: generally less than 2% of 230 Thd and occasionally reaching 6 to 11% in the surface waters of S2 and S1 respectively. Therefore the 230 Thd and 230 Thd-xs profiles do not differ significantly.


M L L V V Gallani

Rastei
email: rastei@ipcms.unistra.fr

C Iacovita

J Hurst

G Manfredi

P A Hervieux

B Donnio

J L Gallani

M V Rastei

Magnetic force fields of isolated small nanoparticles clusters †

Magnetic force fields of isolated small nanoparticle clusters

Introduction

Stable and well characterized magnetic nanoparticles (NPs) are much sought after in many research fields. In particular, bio-related applications making use of iron-oxide NPs (magnetite -Fe 3 O 4 and maghemite -Fe 2 O 3 ) are currently one of the pharmaceutical markets with the fastest growth 1,2 . Such nanoparticles are now intensively used in biomedical imaging, diagnostic and therapeutic applications [3][4][5][6][7] . Their rather easy synthesis [8][9][10] , along with their peculiar magnetic properties 11 and relatively low toxicity 12,13 also make them ideal candidates as contrast agents for magnetic resonance imaging (MRI) [14][15][16][17] , heat generators for cancer treatment by magnetic hyperthermia [18][19][20][21] , as well as carriers in drug delivery systems 22,23 . As always the case for any compound which is meant to be used in human health, there are stringent requirements on the properties of these materials. From the nanometer-scale interaction point of view, the magnetic properties of NPs is of utmost importance for the entire spectrum of applications 24 . For example, the local magnetic field surrounding the NPs is crucial for self-organization of NPs in solutions and on surfaces 25,26 . Or, the magnetic field of NPs is key for the magnetic relaxation-based detection, which is not limited to MRI imaging 27 . Magnetic relaxation-based detection exploits the local magnetic field generated by NPs by altering the spin-spin relaxation time of nearby molecules.

The local magnetic field is nonetheless strongly modified by aggregation, an effect which has been recognized quite early to further enhance the rate of transverse spin relaxation of nearby water molecules for instance [27][28][29][30] . An efficient implementation of NPs in applications is consequently strongly related not only to their saturation magnetization -an active research field in synthetic inorganic chemistry -but also to the magnetic properties of NPs, which has been shown to be impacted by the degree of clustering 22,[31][32][33][34][35] , or the formation of controlled multicore particulate architectures [36][37][38] . The characterization of the local magnetic field generated by either individual or few interacting NPs (hereafter referred to as clusters) is, therefore, also of great interest.

Here, we present an experimental and theoretical study on a quantitative evaluation of the magnetic force fields created by clusters consisting of 1 to 5 NPs. A homogeneous external magnetic field able to saturate the magnetization of NPs has been used. It was found that the local magnetic force fields do not increase proportionally with the number of NPs in the clusters. Moreover, minute modifications of the spatial arrangement of NPs in the cluster were found to drastically modify the magnetic stray field, which in turn has a great impact on the quantitative determination of the magnetic moments of the clusters. In particular, it is shown that at distances smaller than the size of the clusters, the NPs do not equally participate to the measured signal, yielding to an overestimation of the magnetic moment of clusters when using a two-dipole model. This model is nevertheless shown to be suited for distances exceeding some critical values, which can be roughly approximated as twice the diameter of NPs forming the clusters. For measurements performed at larger distances, the application of the two-dipole model leads to a good estimation of the magnetic moments, improving further for clusters having a compact configuration.

Our numerical analyses were carried out by virtually splitting both the tip and the clusters into a large number of microscopic elements, permitting a full pairwise integration of magnetostatic potentials. This method allows taking into account the exact geometry of both tip and sample, which is crucial in the present study. Using the numerical results given by our model in combination with the experimental findings, we bring unprecedented insights into the quantitative evaluation of the magnetic field generated by clusters of NPs. The importance of the tip-sample distance and of the spatial organization of the particulate clusters is clearly emphasized. The results shed light on the magnetic force field distribution near clusters of NPs, which is relevant for a broad range of applications.

Experimental and modeling considerations

The investigation of the local field distribution above magnetic nanostructures of dimensions down to a few nanometers can be experimentally achieved by means of either magnetic force microscopy (MFM) 39 or electron holography 40 . The former technique has been extensively used in the last two decades and is now a mature approach, providing unique information about nanoscale magnetism of on-surface [START_REF] Passeri | Magnetic Force Microscopy, Magnetic Characterization Techniques for Nanomaterials[END_REF][START_REF] Ares | [END_REF] , immersed 43 or embedded nanostructures 44 . MFM uses the attractive and/or repulsive magnetostatic interactions established between the magnetic probe of the microscope and the near-field magnetic stray field emerging from the nanostructure, yielding a force gradient image. The signal is generally measured by converting the local magnetostatic interaction to an amplitude or phase change of the vibrating cantilevered MFM tip 39,[START_REF] Passeri | Magnetic Force Microscopy, Magnetic Characterization Techniques for Nanomaterials[END_REF] . In a first approximation, the phase shift of the cantilever vibration can be written as a function of the gradient force along the z-direction:

∆ϕ = - Q k ∂ F z ∂ z , ( 1
)
where Q is the oscillation quality factor of the cantilever, and k is the cantilever spring constant. For quantitative information, the phase shift needs to be probed at various tip-sample distances. This is in order to allow a 3D imaging of the magnetic field distribution generated by the magnetic nanostructure 39,45 . Moreover, a quantitative interpretation of the phase shift is also complex, since it depends on the detailed magnetic properties of the tip 46,47 . A few experimental studies have nevertheless proven feasible the extraction of quantitative information, such as magnetization or coercive fields, providing that the exact spatial distribution of the tip magnetization is known [48][49][50][51] . Nonetheless, in many studies the two-dipole model 52 was used to interpret the experimental findings [53][54][55] .

The two-dipole model, also called point-probe model, replaces the magnetic nanostructure by a unique magnetic dipole moment, and idealizes the tip magnetization by a magnetic dipole moment which is virtually located in the tip at a certain distance from the apex 52 . Within this model, the phase shift of the cantilever is directly proportional to the nanostructure magnetic dipole moment and to the magnetic dipole moment of the tip 52 . This two-dipole model has been successfully applied for extracting quantitative magnetic information, particularly on nanostructures able to generate field geometries similar to those used for the tip calibration [56][57][START_REF] Sievers | Small[END_REF] . Accordingly, spherical NPs with known magnetization have been employed to calibrate the magnetic tips, which were subsequently used to determine the magnetic moment of other spherical NPs or clusters [START_REF] Sievers | Small[END_REF][START_REF] Angeloni | [END_REF]63 .

Our experiments have been performed in ambient conditions using an ICON-AFM apparatus controlled by a Nanoscope V electronics. The topographic images have been acquired in tapping mode, while the lift mode was used for the magnetic (MFM) imaging. The tip-sample distance has been determined by adding at the mean topographic scan height (30 nm) to the tip-lift values used in the MFM scan. The spring constant of the used cantilevers was between 2 and 3 N/m and their resonance quality factor ranged between 500 and 1000. The probes were silicon tips covered with a hard magnetic coating presenting a coercive field of about 250 Oe. A magnetic field of 0.5 T was permanently applied in all experiments reported in this work. This field is high enough to saturate both the magnetic MFM probe and the NPs. Other details about the experiments are described in Sec. SI3. A characterization of the tip magnetic properties in the framework of two-dipole model is presented in Sec. SI4. The magnetic moments of the NPs evaluated with the two-dipole model are in Sec. SI5.

The numerical simulations were conducted by considering the real geometry of the tip and of the NPs number and organization. The configuration used in simulations is sketched in Fig. SI9. The advantage of our simulation model is that it considers a tip and a sample composed of small magnetic elements (meshes), which better account for the geometry of the systems. A magnetic moment corresponding to the respective volume is associated to each mesh element. Various meshes sizes have been checked, until a convergence was obtained. This typically corresponds to square meshes of a lateral sizes below 2 nm (see Sec. SI6). The experimental height of the tip H was also considered to get a convergence of the calculated force. The elementary force between two elementary magnetic elements of magnetic moments ⃗ M 1 and ⃗ M 2 separated by ⃗ r, reads:

∆ ⃗ F ( ⃗ M 1 , ⃗ M 2 ,⃗ r ) = 3 µ 0 4π r 5 [ ⃗ M 2 ( ⃗ M 1 •⃗ r ) + ⃗ M 1 ( ⃗ M 2 •⃗ r ) + + ⃗ r ( ⃗ M 1 • ⃗ M 2 ) - 5( ⃗ M 1 •⃗ r)( ⃗ M 2 •⃗ r) r 2 ⃗ r ] . ( 2
)
For a saturated magnetization along z-direction, the total magnetic force is then obtained by integrating over the whole system:

⃗ F = ∫ d⃗ r ′ d⃗ r ′′ d ⃗ F ( --→ M NPs ( ⃗ r ′ ) , --→ M Tip ( ⃗ r ′′ ) ,⃗ r ′′ -⃗ r ′ ) , ( 3
)
where m represents the magnetization density. This expression can be numerically evaluated as a function of tip-sample distance z, yielding as well to the first derivative of force versus distance, which can then be compared with the experimental phase shift. Section SI6 is a complete compilation of all parameters used in simulations.

Results and discussions

The high-temperature decomposition of magnetic precursors in organic solvents enables the formation of highly crystalline NPs of various sizes, shapes, compositions and magnetic properties 10,64 . Spherical Fe 3 O 4 NPs with diameters of 30 nm ± 5 nm with high saturation magnetization values have been used (Sec. SI1 and SI2). Prior to the deposition, the NPs were dispersed in chloroform. The residual traces of water in chloroform solvent lead to the assembly of NPs into ring structures of various diameters, as reported elsewhere for gold nanoparticles 65 . In between the ring structures, clusters comprising a different number of NPs can be found (Fig. 1, also Sec.SI1).

Magnetic phase-signal and two-dipole model

Phase shift profiles as a function of tip-sample distance are shown in Fig. 1. The total magnetic moment of the clusters is expected to proportionally increase with the number of NPs composing the cluster, inducing a proportional increase of the phase shift as well. As seen in Fig. 1, this is not the case, since the variation of the phase shift signal does not scale with the number of NPs in the cluster. It is important to notice that for large distances (80 -95 nm: light-brownish curves), the phase signal tends, however, to vary proportionally with the number of NPs. We will see below that the phase shift also depend on the spatial arrangement of NPs in the cluster. In any case, the experimental results presented in Fig. 1 already justify a deeper analysis of how the magnetic field generated by different clusters evolves with the number and organization of NPs. Moreover, experimental observations also reveal the need to develop an accurate numerical model which can be subsequently used for quantitative characterizations of the magnetic fields generated by clusters of NPs.

Stray fields of single NPs

We have pursued our analysis by extracting quantitative information on the magnetic moments of various NPs clusters. This was first realized by characterizing the MFM tip on reference samples consisting in six different isolated NPs having diameters ranging from 26 nm to 34 nm (Fig. 2 and Fig. SI5). This analysis is a necessary step towards obtaining a quantitative information about the magnetic moment of clusters [START_REF] Sievers | Small[END_REF][START_REF] Angeloni | [END_REF] . It is worth noting that the diameters of the NPs have been estimated from height profiles and not from lateral sizes, as lateral dimensions are usually overestimated by tip convolution effects [START_REF] Angeloni | [END_REF] . The convolution effects not being present at the top of a nanoparticle, the extracted diameters of NPs are then in agreement with size histograms obtained from electron microscopy images (Sec. SI2).

Another important issue in non-contact phase imaging is the electrostatic contributions which can interfere with "magnetic" phase shifts 60 . However, since all the measurements presented here are done in the saturation state of the magnetization, the electrostatic contribution is negligible as compared to the magnetic signal [START_REF] Angeloni | AIP Conf. Proc. 2016[END_REF] . A negligible contribution of the electrostatic field is also indicated by the asymptotic variations of the phase signal with the tip-sample distance, which can only be explained by an insignificant gradient of the electrostatic field in the respective tip-distance interval [START_REF] Angeloni | [END_REF] . The magnetic field produced by homogenously magnetized spherical NPs (especially at saturation) is hence well approximated to the field of a point dipole positioned at the center of the spherical NP. This is the core of the two-dipole model as used in MFM, and is quite accurate for individual NPs [START_REF] Angeloni | [END_REF] . The strength of the magnetic field is expected to be proportional to the magnetic moment of the nanoparticle, which in turn increases with its diameter (Sec. SI4). Within this description, a slight deviation of the magnetic moment induces a significant modification of the field above the nanoparticles. Consequently, different volumes of the magnetic tip are involved in collecting the magnetic signal. Evolution of the MFM phase shift with the tip-sample distance is perfectly fitted with the two-dipole model, as seen in Fig. 2(a), provided that the magnetic dipole moment of the tip (m tip ) as well as its virtual location (δ tip ) within the tip volume are free fitting parameters. Note the significant variations with the NP diameter. Since the magnetic moment of each NP is well known, the two fitting parameters can be determined as a function of NP diameter. The resulting values are summarized in Sec. SI4. Those values are helpful for the subsequent analyses of clusters.

Each pair of parameters (m tip , δ tip ) reflects the magnetic force field generated by an isolated NP. Both, m tip and δ tip , are found to increase as the diameter of the particle increases. The values of m tip are of the order of 10 -17 Am 2 , which is in good agreement with known values of magnetic moment of the tips. Plotting the m tip as a function of δ tip results in a linear dependence (Fig. 2). In our case, this dependence might be triggered by the reduced variation in the NPs size, since for fields generated by nanoscale current-carrying parallel lithographed wires, a power-law behavior has been found 57 .

Dimers and role of planarity

The calibration of the magnetic tip enables the extraction of magnetic properties of clusters presenting different spatial organizations. For instance, clusters D1 and D2 in Fig. 3 were formed by assembling two NPs (D1: 30 nm and 28.5 nm; D2: 31 nm and 27.5 nm, respectively, also see Fig. SI6). From those sizes, we can consider that both dimers have an identical volume of magnetic material even if their spatial arrangement is different. The main difference between the two dimers is that, the NPs in D2 do not lie in the same surface plane. They partially overlap (Sec. SI3), appearing in the image in a more compact form than for D1 where both NPs are in the surface plane (Fig. 3).

It is now interesting to realize that the field generated by D2 emerges from a smaller area displaying a larger decay length. This is a simple and clear example of how the NPs in the dimer can impact the MFM signal along a particular direction. Moreover, the phase signal as a function of tip-sample distance (proportionally to the tip lift in MFM) of D2 exhibits a more sloped profile, which means that D2 generates a magnetic force field with a greater gradient magnitude in the normal direction.

In order to quantify the stray field generated by the clusters and to observe at which distances from the cluster surface the two-dipole model increases its accuracy we calculate the ratio between the phase signals of the clusters and of a single NP. As an example, the ratio between the experimental MFM phase shifts recorded for the two dimers (D1 and D2) and a single NP (S4) are displayed in Fig. 3. The NP has been chosen in order to approach as close as possible the mean diameter of the NPs composing the cluster. As seen, for distances above 80 nm, the D2/S4 ratio is indeed almost 2, in agreement with the two-dipole model (which stipulates that if the magnetic moment of the investigated cluster is doubled, the phase shift should also double). However, at lower distances, the ratio decreases down to about 1.5. The D1/S4 ratio starts at 1.7 at the largest distance and smoothly decreases towards 1.2 for the smallest distances [Fig. 3(b)], remaining therefore significantly lower than the ideal value of 2. These findings suggest in both cases that the tip does not perceive the entire cluster magnetic moment, which is a significant deviation from theory.

From the fits, δ tip is 79 and 82 nm for D1 and D2, respec- tively. The δ tip values are larger in comparison to the ones obtained for isolated NPs. By applying the linear dependence from Fig. 2(b), the values of magnetic moment of the tip are 2.95 × 10 -17 Am 2 and 3.23 × 10 -17 Am 2 , when using D1 and D2, respectively. The higher value of the moment observed for D2 is due to the compact geometry, as the stray field interacts with a larger tip volume than in the case of D1. Finally, the extracted values of the magnetic dipole moment of D1 and D2 are 6.91 × 10 -18 Am 2 and 7.69 × 10 -18 Am 2 , respectively. These values represent a percentage of 68.5% and 75.5% of the sum of magnetic moments of the two NPs forming the D1 and D2. Consequently, when individuals NPs are getting closer to each other and start clustering, only a fraction of the entire magnetic moment is detected by the MFM tip, undermining the two-dipole model. Interestingly, in compact geometries, where NPs partially overlaid, the tip starts to perceive the entire magnetic moment for tip-sample distances exceeding roughly twice the diameter of a single NP.

Trimers, tetramers, pentamers and compactness

It is pertinent to verify this trend and to further investigate the effect of the arrangement of NPs in clusters of higher order. The association of a third NP to a dimer leads to a variety of geometrical "trimer" configurations. We analyzed 5 different 3-fold clusters (trimers) (Fig. 4). In Tr1, a NP partially superimposes over the other two, giving rise to a slightly out-of-plane and compact triangular configuration when looking from above; Tr2 has an almost perfect isosceles in-plane triangular configuration, while Tr3 presents a slightly more open triangular shape. Tr4 is less bent, resembling a linear chain, whereas Tr5 has the two outer NPs in part overlaying the central one, forming a more compact geometry from atop perspective. This has consequences for the field gradient generated in the out-of-plane direction, as discussed below in relation with the phase ratio comparison of the trimers. As highlighted in Fig. 4 (b), the phase shift for 3-fold clusters is even more drastically influenced by the NPs arrangement. The evolution of the phase shift as a function of distance is very abrupt in the case of compact Tr1, i.e. the generated magnetic force field has the steepest gradient. The slope of the signal gradually decreases for Tr2 and Tr3, becoming the lowest for Tr4. In others words, as the trimers evolve from a compact triangular configuration to a linear-like chain, the magnitude of the generated magnetic force field decreases. In turn, the force field for Tr5 becomes again comparable with the one of Tr2 and Tr3, effect which is due to the partial overlapping of outer NPs. Hence, we can conclude that the assemblage of NPs in linear-like chains generates reduced magnetic force fields above the center of the chain, whereas compact or overlapping NPs geometries yield gradient fields of greater magnitudes.

The phase shift signals on the trimers have also been normalized by the phase signal of a single NP. Again, the latter were estimated by averaging the diameters of the NPs constituting the respective trimers (Fig. SI7). The magnetic moments of the trimers were evaluated by applying the same procedure as previously used for the dimers. Assuming for instance that the compact Tr1 trimer has a magnetic moment three times larger than the one of a single NP, it is expected that the phase shifts ratio reaches a value of 3 over the entire range of distance. As shown in Fig. 4(b), the ratio is largely below 3 and only tends to this value for distances above 85 nm. For Tr1, the extracted magnetic moment almost equals the sum of the magnetic moments of the individual monomers -90% recovery (Sec. SI5). In the case of Tr2 and Tr3 -where the apparent surface area increases as a result of the larger spacing of the NPs -the recovery percentages are 85.5% and 80%, respectively. The phase shift ratios for Tr2/S3 and Tr3/S2 are close to 3 only for large distances (85 -95 nm), and steeply fall to 1.7 when the distance decreases to 50 nm [Fig. 4(b)]. The recovering percentage reduces significantly to 66% and 70% for Tr4 and Tr5 clusters, respectively. Accordingly, the Tr4/S2 ratio varies very rapidly from 1.3 to 2.6, while Tr5/S5 ratio has a smooth linear trend between 1.5 and 2.

The magnitude of the detected magnetic moment of a cluster therefore significantly depends on the degree of compactness and on the tip-sample distance. For compact triangular clusters the tip detects a full magnetic moment only at distances exceeding 80 nm. As the separation between NPs increases, the recorded magnetic force field decreases, affecting the overall magnitude of the magnetic moment of clusters.

We have expanded our study to larger clusters formed by the assembly of four and five NPs. In the case of the "tetramers" (Te), the NPs are spatially arranged in such way that they are located in the corner of a slightly distorted square, while NPs in "pentamers" (Pe) form a almost regular pentagon (Fig. 5 & Fig. SI7). The phase shift of the pentamer has a slightly larger magnitude than that given by the tetramer (Fig. SI8). In general, we find that the less compact geometries have a great impact on the magnitude of the extracted magnetic moment. The recovery percentage is 64% and 67% for tetramer and pentamer, respectively (Table SI2). This indicates that a larger 2D spatial extension of NPs results in a magnetic force field of a smaller magnitude and a shorter decay length. As can be seen in Fig. 5, where phase shift ratios are plotted for clusters with increasing connectivity, the Te/S4 and Pe/S4 phase shift ratios only approach values of 3.5 and 4.6, respectively, at the largest distances. Only the smaller compact clusters have the tendency to reach a phase shift ratio identical to the number of NPs. These findings again cannot be understood in the framework of the two-dipole model where the clusters are described by point-like magnetic dipole moments.

Numerical simulations

Calculations of the magnetostatic interactions between the microscope tip and the NPs were modeled by using the configurations found experimentally. A detailed description of the method can be found in Sec. SI6. The interaction between the magnetic tip and the sample has been simulated by placing the tip apex above the center of mass of the studied clusters. In Fig. 6 (a), we show the computed magnetic forces as a function of distance for three isolated NPs of different diameters: D = 26 nm, 30 nm and 34 nm. As expected, the magnetic force increases with the size of the NPs while the total magnetic moment scales with the volume. In Fig. 6 (b) are plotted the theoretical force gradient values (which are proportional to the experimental phase shift) with respect to the tip-sample distance. Comparison with the experimental results presented in Fig. 2 reveals that the simulated force gradients exhibit the same trend as in experiments. Indeed, in both simulated and experimental cases, the force gradient increases for larger NPs. The difference between the force gradients however shrinks as the distance between the tip and the NP increases.

Due to the finite element modeling used here, the numeri- gives rise to triangular clusters labeled: θ = π /8 and θ = π /4, respectively. The latter cluster has a more compact structure with the three NPs forming an equilateral triangle θ = π /3. From Fig. 7 (a), it can be seen that the magnetic force exerted on the tip gradually increases as the cluster changes from a linear to a triangular configuration. The differences between the four configurations are as expected more visible at small distances. The more compact the trimer, the stronger the magnetic force. The magnitude of the magnetic force field therefore depends on the geometrical configuration of the trimer. The slope of the force gradient vs. distance is also significantly reduced for low compactness [Fig. 7(b)]. The computed force gradients are consequently in very good agreement with the experimental phase shifts plotted in Fig. 4. Moreover, at large distances (≫ D), the force gradients converge to the same value because all four computed configurations have the same total magnetic moment. In the case of the most compact trimer (θ = π /3), the force gradient shows the greatest variation, as also found experimentally (Tr1 in Fig. 4).

The force gradients shown in Fig. 7 (b) are divided by the force gradient obtained for a single NP. The resulting curves are shown with solid lines in Fig. 8 (a). It is seen that the ratios display a smoother variation as the trimers adopt a less compact geometry. It is worth recalling that these ratios indicate how the magnetic fields generated by a cluster deviates from the fields generated by a single NP. As a matter of fact, the ratio cannot be larger than three for a trimer.

When the Tr/S ratio reaches the maximal value of three, it means that the magnetic force field is not any more dependent on the spatial arrangement in the trimer. In our numerical analyses, we found a maximum value of 2.8 at 95 nm in the case of the most compact trimer (θ =π /3). This maximum value decreases considerably at the same distance for less compact trimers. It means that, even at distances larger than 2-3 times the NPs diameter, the tip is still sensitive to the 2D spatial arrangement of the NPs in the trimer.

For the sake of comparison, we also calculated the force gradient ratios for θ = 0 and θ = π /3 trimers using the twodipole model. The results are shown with dashed lines in Fig. 8(a). A full analysis performed with the two-dipole model is presented in Sec. SI7. At large distances > 2 D, the ratio for the compact trimer is greater than the ratio of linear trimer, which is qualitatively in agreement with both previous simulations and experiments. For small distances (< 2D), the ratios show a completely different variation as compared with the respective curves obtained with our finite element model and with the experiments. For instane, at distances of the order of D, the two dipole model predicts a ratio for the linear-like chain cluster higher than the one for the compact trimer, aspect which is not observed in our finite-element modeling nor in the experiments. The ratio even becomes negative at very short distances (< D), meaning that the force between tip and clusters is attractive. This happens because the tip dipole moment is so close to the dipole moment of NPs that there is an angle where the sign of the magnetic force changes. This is not the case in our finite element modeling neither in the experiments, because of the finite size of the NPs and of the tip. Altogether, these results already show that two-dipole model cannot capture a reliable magnetostatic interaction at separation distances < 2D.

The force gradients of several clusters ranging from two to five NPs divided by the force gradient of a single NP are shown in Fig. 8(b). These theoretical ratios can be compared with the experimental curves presented in Fig. 5. As seen, the numerical results support the experimental findings. Depending on the number of NPs composing the clusters, the ratios progressively increase as the distance increases, approaching a saturation value. Above approximately 70 nm, the ratios are almost 2 and 3 for dimers and trimers, respectively. Instead, at the same 70 nm separation distance, the ratios for tetramer and pentamer are only 3.6 and 4.5, saturating at 3.9 and 4.8 only for very large distances. This is due to the fact that the tip senses differently the magnetic force fields depending on the 2D organization of the NPs in the cluster, as also observed experimentally.

Summary and conclusions

We presented a combined experimental and theoretical study on a quantitative evaluation of the magnetic force fields generated by clusters of NPs formed by associations of single magnetic NPs on a solid non-magnetic substrate. The experimental results obtained under a saturation magnetic field of 0.5 T have shown that the magnetic force fields do not increase proportionally with the number of NPs in the clusters. Furthermore, for dimeric and trimeric clusters, minute modifications of the 2D spatial arrangement of NPs in the clusters drastically modify the surrounding magnetic field. For these clusters, the magnitude of the generated force field follows the prediction of the two-dipole model, but only at distances exceeding some critical values (about twice the diameter of a NP). For tetramers and pentamers the two-dipole model is found to overestimate the value of the magnetic moments in the whole range of the experimental distances. The largest the in-plane lateral expansion of the cluster the less accurate is the two-dipole model in extracting quantitative information. This is because the principal condition for the use of a twodipole model is related to the lateral extension the clusters which should be smaller than the distance at which the magnetic force field is evaluated. A good agreement with the experimental results, in the entire range of tip-sample distance, was obtained by considering the real geometry of the system, namely by splitting both tip and sample in a finite number of microscopic elements and performing a full pairwise integration of magnetostatic potentials. The numerical calculations conducted with this finiteelement model brought unprecedented insights into the quantitative evaluation of the magnetic field generated by assemblies of NPs, and put forward the importance of the spatial organization. The magnetic force field in the proximity of clusters has a complex distance-dependent structure, being driven significantly by the exact number of NPs. Compact clusters have a general tendency to generate an enhanced force field able to extend at larger distances, which is the result of interparticle magnetic interactions and cooperative effects. Our results may have an impact on the chemical design of NPs-based clusters for various applications, including biomedical and therapeutic strategies or magnetic information technologies.

Fig. 1

 1 Fig. 1 Top row: AFM topographic images (200 × 200 nm 2 ) of a single NP and of various clusters made from two to five iron-oxide nanoparticles. The vertical scale bar is 34 nm. Middle row: MFM phase shift profiles (horizontal direction) as a function of tip-sample distance [tip lifts from 20 nm (dark blue) to 65 nm (light brown), values measured with respect to the topographic scan which is at about 30 nm from the sample surface]. Lower row: Phase shift profiles (vertical direction). Profiles are vertically offset for clarity.

2 )Fig. 2

 22 Fig. 2 (a) Phase shift of six single NPs of various diameters. Solid lines are fits obtained with the two-dipole model. Insets are topographic images of the investigated NPs (150 × 80 nm 2 ). (b) Plot of the magnetic dipole moment of the tip as a function of its virtual location in the tip volume. Blue line is a linear fit.

Fig. 3 (

 3 Fig. 3 (a) Phase shift of two dimers (D1 and D2) of different organization. Solid lines are fits obtained with the two-dipole model. The insets are topographic images of D1 and D2 (150 × 100 nm 2 ). (b) Distance dependence of the phase shift for the two dimers divided by the phase shift of a single NP.

Fig. 4

 4 Fig. 4 (a) Phase shift of five 3-fold clusters (trimers) of various compactness. Solid lines are fits using the two-dipole model. The insets are topographic images of the trimers Tr1 (scale bar: the large side of images is 170 nm). (b) Phase shifts for the five trimers divided by the phase shift of a single NP, reported as a function of tip-sample distance.

Fig. 5

 5 Fig. 5 Phase shift of several clusters divided by the phase shift of a single NP. The NP has been chosen with a diameter close to the mean diameter of those composing the cluster. Insets are topographic images of the analyzed clusters.

Fig. 6 ( 4 . 1 7 .

 6417 Fig.6(a) Variation of the magnetic force between the tip and three single NPs of different diameters (D = 26 nm, 30 nm and 34 nm) as a function of tip-NP distance. (b) Variation of the force gradient as a function of distance for the three different single NPs. Inset exhibits the variation of the gradient force ratio between the larger singles NPs (D = 30 and 34 nm) and the smallest single NP (D = 26 nm) as a function of distance by using our simulation method (full curve) and an equivalent two-dipole model (dashed curve).

Fig. 7 (

 7 Fig. 7 (a) Force and (b) force gradient between the tip and the four clusters depicted in (c). (c) Top view of four trimers of different compactness. The red crosses represent the lateral position of the tip.

Fig. 8 (

 8 Fig.8(a) Force gradients for the four trimers divided by the force gradient of a single NP (30 nm in diameter). For comparison, the results for θ = 0 and θ = π /3 as computed with the two-dipole model are shown with dashed lines. (b) Force gradients for clusters composed from two to five NPs divided by the force gradient of a NP.

Acknowledgments

We thank N. Beyer for technical assistance, O. Ersen for TEM measurements, S. Colis for magnetic hysteresis measurements, and R. Dudric for XRD data. This work is supported by ANR METABIP ( 12BS10 003 01), ANR OH-RISQUE SMARAGD (14 OHRI 0008 01), EOARD (FA8655-13-1-3001) and partially by ANR-11-LABX-0058-NIE within the Investissement d'Avenir program ANR-10-IDEX-0002-02. C. Iacovita acknowledges financial support from Romanian National Authority for Scientific Research, CNCSIS-UEFISCDI, through researcher mobility project No. PN-III-P1-1.1-MC-2017-0981 and research project to stimulate young independent teams No. PN-III-P1-1.1-TE-2016-0967. We thank CNRS & University of Strasbourg for support.