Lucas Hardouin

Christine Delire
email: christine.delire@meteo.fr

Bertrand Decharme

David M Lawrence

Julia E M S Nabel

Victor Brovkin

Nathan Collier

Rosie Fisher

Forrest M Hoffman

Charles D Koven

Roland Séférian

Tobias Stacke

Uncertainty in land carbon budget simulated by terrestrial biosphere models: the role of atmospheric forcing

Keywords: uncertainty, land carbon sink, atmospheric forcing dataset, terrestrial biosphere model

Global estimates of the land carbon sink are often based on simulations by terrestrial biosphere models (TBMs). The use of a large number of models that differ in their underlying hypotheses, structure and parameters is one way to assess the uncertainty in the historical land carbon sink.

Here we show that the atmospheric forcing datasets used to drive these TBMs represent a significant source of uncertainty that is currently not systematically accounted for in land carbon cycle evaluations. We present results from three TBMs each forced with three different historical atmospheric forcing reconstructions over the period 1850-2015. We perform an analysis of variance to quantify the relative uncertainty in carbon fluxes arising from the models themselves, atmospheric forcing, and model-forcing interactions. We find that atmospheric forcing in this set of simulations plays a dominant role on uncertainties in global gross primary productivity (GPP) (75% of variability) and autotrophic respiration (90%), and a significant but reduced role on net primary productivity and heterotrophic respiration (30%). Atmospheric forcing is the dominant driver (52%) of variability for the net ecosystem exchange flux, defined as the difference between GPP and respiration (both autotrophic and heterotrophic respiration). In contrast, for wildfire-driven carbon emissions model uncertainties dominate and, as a result, model uncertainties dominate for net ecosystem productivity. At regional scales, the contribution of atmospheric forcing to uncertainty shows a very heterogeneous pattern and is smaller on average than at the global scale. We find that this difference in the relative importance of forcing uncertainty between global and regional scales is related to large differences in regional model flux estimates, which partially offset each other when integrated globally, while the flux differences driven by forcing are mainly consistent across the world and therefore add up to a larger fractional contribution to global uncertainty.

Introduction

During the last decade, about 45% of anthropogenic carbon dioxide (CO 2 ) emissions remained and accumulated in the atmosphere. Model-and observation-based studies suggest that the remaining term was shared among ocean (24%) and land (32%) since 1960 [START_REF] Friedlingstein | Global carbon budget[END_REF]. Together, the contribution of the ocean and terrestrial CO 2 uptake nearly halves the increase in atmospheric CO 2 , damping the pace of climate change [START_REF] Canadell | Global carbon and other biogeochemical cycles and feedbacks Climate Change[END_REF]. These sinks significantly increased since the middle of the 20th century, mainly due to the acceleration of the increase in CO 2 concentrations, caused by fossil fuel emissions. However, carbon uptake processes depend strongly on climate variability, especially for the terrestrial biosphere (Le [START_REF] Quéré | Trends in the sources and sinks of carbon dioxide[END_REF][START_REF] Devries | Decadal trends in the ocean carbon sink[END_REF]. While the ocean sink grew from 1.1 ± 0.4 PgC•yr -1 in the 1960s to 2.8 ± 0.4 PgC•yr -1 in 2011-2020, with an inter-annual variability of a few tenths of PgC•yr -1 , the land sink rose from 1.2 ± 0.5 PgC•yr -1 to 3.1 ± 0.6 PgC•yr -1 , but with interannual variations up to 2 PgC•yr -1 [START_REF] Friedlingstein | Global carbon budget[END_REF].

Accurate evaluations of anthropogenic emissions, atmospheric CO 2 levels and carbon cycle perturbation are necessary to monitor, understand and predict climate change. The Global Carbon Project (GCP) has published an annual report since 2013 (Le Quéré et al 2013, 2014[START_REF] Friedlingstein | Global carbon budget[END_REF], 2021) that quantifies the magnitude and uncertainty of the five major components of the global carbon budget: fossil fuel emissions (E FOS ), emissions from land use change (E LUC ), the growth rate of atmospheric CO 2 (G ATM ), the ocean sink (S OCEAN ) and the land sink (S LAND ). In earlier analyses provided by the GCP, the terrestrial carbon cycle simulated by the terrestrial biosphere models (TBMs, called by the authors Land Biosphere Models or Dynamical Global Vegetation Models) was not considered reliable enough to be used in the report to estimate the land carbon sink (Le [START_REF] Quéré | Trends in the sources and sinks of carbon dioxide[END_REF][START_REF] Schaefer | A model-data comparison of gross primary productivity: results from the north american carbon program site synthesis[END_REF][START_REF] Todd-Brown | Causes of variation in soil carbon simulations from CMIP5 earth system models and comparison with observations[END_REF]. The land sink S LAND was hence diagnosed indirectly as the residual of the other terms. However, the methodology was updated in the 2017 global carbon budget [START_REF] Quéré | Global carbon budget[END_REF], partly due to improvements in carbon cycle representation [START_REF] Collier | The international land model benchmarking (ILAMB) system: design, theory and implementation[END_REF], Lawrence et al 2019[START_REF] Arora | Carbon-concentration and carbon-climate feedbacks in CMIP6 models and their comparison to CMIP5 models[END_REF][START_REF] Davies-Barnard | Nitrogen cycling in CMIP6 land surface models: progress and limitations[END_REF], but also because of evidence of underestimation and uncertainty in the ocean sink variability [START_REF] Landschützer | The reinvigoration of the southern ocean carbon sink[END_REF][START_REF] Devries | Recent increase in oceanic carbon uptake driven by weaker upper-ocean overturning[END_REF]. The global land sink is now estimated by the multi-model mean of the TBM simulations, with the budget imbalance term, B IM , representing residual uncertainty and/or unexplained aspects of the actual global carbon cycle.

To account for model construction, structural, and parametric uncertainty, as many as 17 TBMs were used to estimate S LAND by [START_REF] Friedlingstein | Global carbon budget[END_REF]. Simulations were run using a unified land use change data set (LUH2; Hurtt et al 2020), global atmospheric CO 2 trend, and atmospheric climate forcing (CRUJRA; Harris 2019). Use of a single forcing dataset for climatic forcing, however, means that any model spread due to climate forcing data uncertainty is not represented by the GCP ensemble. Yet, several studies have identified historical climate forcing as a large source of uncertainty in terrestrial carbon cycle modelling [START_REF] Hicke | NCEP and GISS solar radiation data sets available for ecosystem modeling: description, differences and impacts on net primary production[END_REF][START_REF] Jung | Uncertainties of modeling gross primary productivity over Europe: a systematic study on the effects of using different drivers and terrestrial biosphere models Glob[END_REF][START_REF] Poulter | Impacts of land cover and climate data selection on understanding terrestrial carbon dynamics and the CO2 airborne fraction[END_REF][START_REF] Bonan | Model structure and climate data uncertainty in historical simulations of the terrestrial carbon cycle (1850-2014)[END_REF], Lawrence et al 2019). [START_REF] Hicke | NCEP and GISS solar radiation data sets available for ecosystem modeling: description, differences and impacts on net primary production[END_REF], for example, highlights important biases in net primary productivity (NPP) estimation while using different radiation datasets, and [START_REF] Poulter | Impacts of land cover and climate data selection on understanding terrestrial carbon dynamics and the CO2 airborne fraction[END_REF] conclude that atmospheric forcing results in a large uncertainty compared to land-cover datasets for NPP, heterotrophic respiration (Rh) and net ecosystem exchange. More recently, as part of the assessment and benchmarking of the Community Land Model version 5 (CLM5) [START_REF] Lawrence | The community land model version 5: description of new features, benchmarking and impact of forcing uncertainty[END_REF], climate forcing uncertainty was compared to model structure uncertainty using two forcing data sets (Global Soil Wetness Project 3 and CRUNCEP) and three versions of CLM (CLM4, CLM4.5, and CLM5) that differ markedly in their carbon cycle representation [START_REF] Bonan | Model structure and climate data uncertainty in historical simulations of the terrestrial carbon cycle (1850-2014)[END_REF]. The authors concluded that climate forcing is a large source of uncertainty in the global carbon cycle, especially for GPP, NPP and Rh, and to a lesser extent for Net Biome Productivity. However, this study was performed with different versions of a single model, and only two estimates of atmospheric climate forcing. Arguably therefore, the importance of atmospheric forcing on carbon cycle estimation by TBMs remain insufficiently explored, especially in the context of recent model structural improvements in the 6th Coupled Model Intercomparison Project (CMIP6) generation of models.

In this study, we used the output from the TBMs of three Earth System Models (ESMs) participating in the CMIP6. We performed an analysis of variance similar to [START_REF] Hawkins | The potential to narrow uncertainty in regional climate predictions Bull[END_REF], [START_REF] Lovenduski | Reducing uncertainty in projections of terrestrial carbon uptake[END_REF], [START_REF] Bonan | Model structure and climate data uncertainty in historical simulations of the terrestrial carbon cycle (1850-2014)[END_REF]in order to quantify the relative uncertainty from atmospheric forcing, models and forcing-model interactions in carbon flux estimates made by TBMs. We used a combination of three TBMs from ESM and three atmospheric forcing datasets in our analysis. We present results for the major terrestrial carbon cycle fluxes; gross primary productivity (GPP), autotrophic respiration (Ra) and heterotrophic respiration (Rh), at global and regional scales. We also investigate estimates of CO 2 emissions from natural fires (fFire) and highlight their role in the calculation of the land carbon sink, before focusing on sources of uncertainty in the global net ecosystem productivity (NEP = GPP-Ra-Rh-fFire).

Methods

CMIP6 models and atmospheric forcings

We used the results of three TBMs of CMIP6 ESMs described in table 1, from the Land Surface, Snow and Soil moisture Model Intercomparison [START_REF] Hurtt | Harmonization of global land use change and management for the period 850-2100 (LUH2) for CMIP6[END_REF]. Those three models have been selected because they are the only ones that report output from three different atmospheric forcings and also share data for emissions from fire, on the Earth System Grid Federation archive.

Three climate forcing datasets were used to force the TBMs, according to the LS3MIP protocol. Atmospheric forcing consists of hourly to 6 hourly data on precipitation, short and long wave solar radiation, near surface air temperature, specific humidity, and wind speed. The three forcings described in table 2 have a 0.5 • × 0.5 • spatial resolution.

For reference, results of analysis of this set of model simulations with the International Land Model Benchmarking (ILAMB, Collier et al 2018) package are provided here: www.ilamb.org/land-hist/. ILAMB produces systematic evaluation (plots, summary tables, scoring) of model output against observations for a range of metrics including bias, RMSE, pattern correlation, annual cycle phase, and variableto-variable comparisons. The version of ILAMB presented here assesses 22 land carbon, water, and energy cycle variables model variables against more than 50 observational datasets.

Analysis of variance

To statistically disaggregate the impacts of climate forcing, model choice and their interactions on relevant output variables, we performed an ANOVA. We analysed nine terrestrial carbon cycle simulationsthree TBMs forced by three forcings-for the time period 1960-2012, chosen to coincide with the beginning of atmospheric carbon dioxide measurements in the 1960s, and the end of the Princeton dataset. We quantify relative uncertainty into three sources: models, atmospheric forcing differences and modelforcing interactions which measure the fact that models behave differently to different forcing product. Details can be found in the supplementary material.

Results

Global scale

We first investigated and compared the forcing data in order to highlight the main differences and similarities. The three forcings mostly agree on near surface air temperature, but they present large differences in specific humidity, incoming shortwave and longwave radiation, and windspeed (figure 1). Princeton has much higher specific (and relative-not shown) humidity than CRUJRA or GSWP3 in most parts of the world, except in desert areas. Additionally, specific humidity is generally greater with GSWP3 than with CRUJRA in tropical regions. CRUJRA also has higher shortwave and lower longwave incoming radiation than GSWP3 and Princeton, both having almost identical values. The precipitation forcing differs also with higher values for GSWP3, especially within the Arctic circle related to snowfall uncertainties.

We performed the analysis of variance on the global estimates of GPP, Ra, Rh, GPP-Ra-Rh, fFire and NEP (figure 2). The results illustrate a dominant role of the atmospheric forcing variation on global GPP estimates (figure 2(a)), with large differences between the multi-model mean values according to the forcings (about 15 PgC•yr -1 ), with GSWP3 resulting in lowest GPPs for all three models, Princeton the highest GPPs, and CRUJRA in between. The discrepancies between the models for a given forcing are smaller, especially for CRUJRA and Princeton. Global GPP estimates of the ISBA-CTRIP and CLM5 models are very similar when forced with the CRUJRA forcing, and slightly differ when forced with GSWP3 or Princeton. The global GPP estimates of JSBACH are the largest with all three forcings and the difference with the two other models is the biggest with the GSWP3 forcing. However, the three models obtain very similar results with the Princeton atmospheric forcing. This suggests a general sensitivity of the models to an atmospheric variable that is different between Princeton's and the two other forcings, such as specific humidity (figure 1). These nine simulation results indicate a strong sensitivity of modelled GPP to climate forcing.

Figure 2(b) confirms the dominant role of atmospheric forcing on global GPP uncertainty. It accounts for ∼75% of the variability in the mean estimates, against ∼19% for models and ∼6% for interactions terms.

The dominance of atmospheric forcing contribution over the uncertainty is even more obvious for autotrophic respiration (figure 2(c)), again with similar model mean values for a given forcing but large differences between forcings with no overlap. According to our analysis of variance, ∼90% of the uncertainty can be assigned to forcings, ∼5% to models and ∼5% to the interactions terms (figure 2(d)). The Princeton forcing results in a slightly larger spread between the average autotrophic respiration than the other forcings.

Global NPP (GPP-Ra) behaves differently. Despite differences in multi-model mean values per atmospheric forcings (around 5 PgC•yr -1 ) attesting to the role of forcing in uncertainty, we can see that JSBACH estimates are clearly higher than the two others (more than 10 PgC•yr -1 ) in figure 2(e). This indicates greater model contribution to uncertainty, confirmed by figure 2(f) where models are the dominant source (∼63%) of uncertainty for NPP, but atmospheric forcing remains important with ∼30% of variability explained and ∼7% for modelforcing interactions. Contrary to ISBA-CTRIP and CLM5 that simulate the lowest NPP with GSWP3 and the highest with Princeton (like for GPP and Ra), JSBACH simulates a slightly smaller global NPP with the CRUJRA forcing. Global Rh presents similar results to global NPP (not shown), both for atmospheric contribution to uncertainty and for greater values obtained with GSWP3 than CRUJRA for JSBACH.

We show that atmospheric forcings have an important (NPP, Rh) and even dominant (GPP, Ra) role in uncertainty distribution for the major fluxes of global carbon cycle estimates. It seems that forcing contribution to uncertainty decreases when we consider net flux such as NPP and Rh rather than direct products of photosynthesis (GPP and Ra).

Wildfires are a small contribution to the global land carbon loss fluxes (∼2-3 PgC•yr -1 ) (Van der Werf et al 2010), especially when compared to plant and soil respiration. Nonetheless, fire carbon losses correspond to approximately one-third of the GPP-Ra-Rh net flux. Moreover, this flux has particularly high seasonal and inter-annual variability, in response to both land cover change, as well as climatic events like drought. Despite its importance, anthropogenic biomass burning is not yet represented by all TBMs. Overall, high uncertainty remains in the carbon emissions associated with fire (figures 2(g) and (h)). ISBA-CTRIP tends to overestimate global fire emissions compared to the global fire emission database (2.14 PgC•yr -1 for 1997-2015 van der Werf et al, 2017), while JSBACH is slightly lower (especially with Princeton). Part of these differences may be explained by the intermediate complexity fire module GlobFirm [START_REF] Thonicke | The role of fire disturbance for global vegetation dynamics: coupling fire into a dynamic global vegetation model[END_REF] used in ISBA-CTRIP that is known to overestimate fire emissions [START_REF] Li | A process-based fire parameterization of intermediate complexity in a dynamic global vegetation model[END_REF]. Figure 2(g) shows a large discrepancy between model values for a given forcing, particularly for the Princeton one, and figure 2(h) confirms a fractional contribution to uncertainty of ∼85% for the models, ∼6% for forcings and ∼9% for interactions. A closer look at fire results shows indeed interesting interactions between the forcings and fire modules that explain the bigger dispersion for Princeton. We saw that GPP and NPP (GPP-Ra) is enhanced with the Princeton forcing, most probably because of the high level of air humidity. For the GlobFirm module implemented in ISBA-CTRIP, this leads simply to lower evapotranspiration rates, less moisture stress, higher productivity and litterfall (not shown) resulting in more biomass burned and greater associated carbon emissions. However, with the more mechanistic fire module of CLM5, the higher air humidity results directly in lower level of fuel flammability (see equation ( 8) Li et al 2012). In JSBACH, finally, the three applied forcings cause very different fire-vegetation feedbacks, such that the natural vegetation cover in the three simulations differs strongly, particularly in fire prone regions such as Africa and Asia (not shown). Such strong fire-vegetation feedbacks have already been observed for JSBACH in previous idealised model simulations [START_REF] Lasslop | Multiple stable states of tree cover in a global land surface model due to a fire-vegetation feedback[END_REF]. Emissions from fire appear to be quite sensitive to atmospheric forcing when we compare individually the relative dispersion for each model, but this dependency is overshadowed in this analysis by the even larger uncertainty related to the model structure.

Here, we first present results for GPP-Ra-Rh estimates. This flux takes into account the most important and probably best represented processes of the terrestrial biosphere carbon sink, by excluding fire emissions. As expected, this net flux is smaller and has larger interannual variability (standard deviation) relative to the mean than the gross fluxes GPP, Ra and Rh. Figure 2(i) suggests a more important contribution to uncertainty by forcing than models: model mean values are close to each other for a given forcing and there is a small difference between forcings. The analysis of variance results shows effectively a dominant role of forcings (∼52%) over models (∼8%) but also a substantial contribution of the interaction terms (∼40%) (figure 2(j)). Indeed, the models do not show similar patterns with respect to the forcings. CLM5 results in greater values for CRUJRA (6.53 PgC•yr -1 ) than Princeton (6.14 PgC•yr -1 ), while ISBA-CTRIP simulates slightly bigger fluxes with Princeton (6.33 PgC•yr -1 ) than with CRUJRA (6.23 PgC•yr -1 ), as does JSBACH with a greater difference (6.06 for CRUJRA against 6.68 PgC•yr -1 with Princeton). In addition, we observe only slightly bigger results for CRUJRA than GSWP3 for JSBACH, consistent with the results obtained from Rh and NPP, while the difference is more pronounced for ISBA-CTRIP and CLM5. This highlights differences in model sensitivity to atmospheric forcings, such as air humidity or downward radiation flux, for the land sink calculation.

Larger differences are obtained when we take into account fFire in the NEP calculation (figures 2(k) and (l)). There is only a small difference between the multi-model mean values per forcing, and a larger spread between the values of the different models, especially with the Princeton forcing (figure 2(k)). Moreover, global NEP is almost independent of the forcing with CLM5 and ISBA-CTRIP, which is not the case with JSBACH. The variance analysis results show that models contribute ∼75% to uncertainty, forcings account for ∼9% and interactions for ∼16%. This contrasts with the significant forcing contribution we obtained when not accounting for fFire, and indicates clearly a high model uncertainty related to fire emission.

To complete the global picture we analyze the contribution of atmospheric forcing to uncertainty in carbon stocks (figure 3). Atmospheric forcing plays a role in global carbon stock uncertainties but not a dominant one like for gross carbon fluxes. There are large differences in live (cVeg) and dead biomass (cLitter), and soil organic carbon content (cSoil) between multi-model means per forcing but the differences between models are even larger. The analysis of variance for cVeg shows that ∼26% of the uncertainty can be assigned to forcings, ∼67% to models and ∼7% to the interactions terms (figure 3(b)). For cLitter and cSoil, these numbers are respectively ∼20% and ∼9% for forcing, ∼68% and ∼81% for models and ∼13% and ∼10% for interactions (figures 3(d) and (f)). CLM5 and ISBA-CTRIP both simulate the smallest live and dead stocks with GSWP3 while JSBACH simulates the smallest stocks with CRUJRA(figures 3(a), (c) and (e)). Princeton results in the largest stocks with every model. Live biomass from CLM5 and ISBA-CTRIP are fairly similar while JSBACH's is lower. The litter reservoir is the highest with ISBA-CTRIP and the lowest with CLM5 while the soil carbon is the largest with CLM5 especially for CRUJRA and Princeton. The higher soil carbon content with CLM5 can be explained by the vertically discretized soil carbon module with deeper soils (down to a maximum of 8.5 m). In the case of ISBA-CTRIP and JSBACH the soil carbon reservoir only represents the carbon in the first meter.

Regional scale

We performed the same analysis of variance on the 30 sub-continental regions defined in the Intergovernmental Panel on Climate Change (IPCC) Special Report on Managing the Risks of Extreme Events and Disasters to Advance Climate Change Adaptation (SREX; IPCC SREX 2012). Figure 4 presents the results for GPP, Ra and Rh. The regions that present the largest uncertainty (largest pies) for GPP (figure 4(a) and supplementary tables S2-S4) are North Asia (NAS), the Amazon forest (AMZ), and to a lesser extent West Africa (WAF), South Asia (SAS), South East Asia (SEA), East Africa (EAF), South Africa (SAF), and West North America (WNA) even if WNA contributes only 2-4 PgC•yr -1 to global GPP. Atmospheric forcing is a dominant source of uncertainty (more than half) in some regions that have fairly large contributions to the global total GPP flux, such as EAF and East Asia (EAS), but also in Central North America (CNA), West Coast South America (WSA), and North East Brazil (NEB). Atmospheric forcing is an important source of uncertainty (more than a third) in AMZ and WAF, tropical forest regions, which contribute greatly to global GPP. Generally, it seems that atmospheric forcings play an important role in tropical regions. Less important contribution in SAS and SEA could be explained by monsoon regimes that prevail on atmospheric forcing differences, but also by model resolution that could explain dominant model uncertainty in island regions (SEA but also Caribbean and Pacific Island Regions). Model-forcing interactions play a secondary but significant role in some regions such as NAS, WAF and South East South America (SSA) for the ones that contribute the most to global GPP.

Slightly different results can be observed for Ra (figure 4(b)). Atmospheric forcing contribution to uncertainty is lower in NAS and AMZ, as well as SEA and NEB. Generally, the forcing and interaction terms seem to play a less important role in Ra than in GPP regional estimates. For Rh (figure 4(c)), forcing contribution is even more reduced in tropical regions, where the fractional distribution never exceeds onethird. However, atmospheric forcing plays an important role in heterotrophic respiration uncertainty in the northern hemisphere's mid and high latitudes (about one-third). Interaction terms become also important in numerous mid and northern regions, particularly in NAS, CNA and Central Europe (CEU). An interpretation could be that models do not predict the same evolution for Rh in the areas where atmospheric variables are changing faster.

It is interesting to probe why atmospheric forcing accounts for more uncertainty at global than regional scale. We found that some differences between the forcings have the same impact on most of the world's area, while model differences induce strong regional uncertainty but not in the same manner everywhere. As an example, GPP estimates are slightly larger with Princeton than with CRUJRA or GSWP3 in all regions (see supplementary material figure S1), regardless of the model and add up to a large difference at global scale. This is consistent with the idea that the higher humidity in the Princeton dataset allows for higher photosynthetic rates for all models forced by it. Conversely, GPP estimates by ISBA-CTRIP and CLM5 appear to be similar at global scale, yet the regional analysis highlights strong opposition between tropical regions with higher GPP for ISBA-CTRIP compared to CLM5, and the northern hemisphere which has higher GPP in CLM5 compared to ISBA-CTRIP (see supplementary material figure S2). This leads to a strong model contribution to uncertainty at regional scale, but those differences compensate each other globally.

Conclusion

The goal of this study was to evaluate the uncertainty related to the use of different global atmospheric forcing datasets in terrestrial carbon cycle modelling. We used nine CMIP6 simulations from three TBMs run with three different forcings to perform our analysis. First, we focused on global averages and found that atmospheric forcings are a dominant source of uncertainty compared to model choice for GPP and Ra, and contribute significantly to the overall uncertainty for Rh and NPP. The contribution of atmospheric forcing to the net fluxes uncertainty is reduced since positive effects of atmospheric forcings on photosynthesis are partially offset by enhancement in the respiration fluxes. We still found an important role of atmospheric forcing on GPP-Ra-Rh temporal mean estimates, but with important interaction terms that translate different model responses to climate variables. Important differences in fire modelling lead to a dominant role of models in NEP uncertainty, even on a global scale, despite a visible relationship and interactions between fire emission and climate forcings.

Secondly, we looked at the partitioning of uncertainty at the regional scale, with the purpose of identifying the regions where atmospheric forcings contribute the most to the variability of GPP, Ra and Rh. We showed that generally, the model structure is the dominant source of uncertainty regionally, in contrast to what we found globally. Atmospheric forcing contribution remains significant and even slightly dominant in some regions, notably in tropical forests for GPP and Ra and in mid and northern high latitudes for Rh, but regional discrepancies among the models are stronger. However, and in contrast with the forcings, it is not the same models that result in the biggest and lowest values everywhere. Those differences add up and offset each other, leading to closer results globally.

While the purpose of the study was to investigate the contribution of forcing uncertainty and to demonstrate the influence of choosing a specific forcing over another forcing, this experimental design, using three distinct sets of forcing data fields is not ideal for identifying the relative sensitivity to which drivers are most critical. While the large difference in humidity between Princeton and the other two datasets points to that field being particularly important in governing GPP and other fluxes, a more specific one-at-a-time propagation of the uncertainty in each field may allow for a more specific attribution of uncertainty to meteorologic variables. Likewise, perturbations to the distributions while holding mean values constant may allow attribution of carbon cycle sensitivity to extremes. Nonetheless, these results point to a focus on constraining humidity values as key to reducing forcing uncertainty. To the extent that humidity is a key driver in GPP uncertainty, another question is whether this uncertainty is exacerbated in offline simulations such as used here, GCP, and elsewhere, because of the inability for surface moisture fluxes to attenuate humidity biases in the forcing fields due to the one-way coupling.

Additional work involving more TBMs and climate forcings should be done in order to better quantify the role of forcings in carbon cycle uncertainty, especially for the interactions between fire, land-cover dynamics and climate forcings. Even so, we conclude that atmospheric forcings are a key source of uncertainty in carbon cycle modelling at the global scale and are a significant source of uncertainty in some regions. Therefore, we suggest that where possible, it would be preferable for future MIPs and assessments (e.g. TRENDY, GCP, CMIP) to run simulations with several alternative climate forcing datasets, and to more specifically generate datasets to allow attribution of carbon cycle sensitivity both to uncertainty in meteorologic fields and to uncertainty in mean versus extreme values of meteorologic fields, when estimating the terrestrial carbon sink in order to correctly represent the uncertainty associated.

Figure 1 .

 1 Figure 1. Annual mean of atmospheric variables used to drive terrestrial biosphere models: GSWP3 (gswp, left-most column), Princeton minus GSWP3 (pgf-gswp; center left column), CRUJRA minus GSWP3 (crujra-gswp; center right column), and zonal averages of each product (right-most column).

Figure 2 .

 2 Figure 2. Model means and standard deviations (color dots and bars) and multi-model means (black crosses) of global carbon flux time series for each atmospheric forcing between 1960 and 2012 (pointplot) and fractional distribution of uncertainty (variability) obtained from the ratios of the sum of squared deviation from the mean (equations (3)-(6) in supplementary material) (pieplot) for GPP (a), (b), Ra (c), (d), NPP (e), (f), fFire (g), (h), GPP-Ra-Rh (i), (j) and NEP = GPP-Ra-Rh-fFire (k), (l). Standard deviations are calculated from the yearly time series of the fluxes.

Figure 3 .

 3 Figure 3. Model means and standard deviations (color dots and bars) and multi-model means (black crosses) of global carbon stocks time series for each atmospheric forcing between 1960 and 2012 (pointplot) and fractional distribution of uncertainty (variability) obtained from the ratios of the sum of squared deviation from the mean (equations (3)-(6) in supplementary material) (pieplot) for cVeg (a), (b), cLitter (including coarse woody debris) (c), (d), and cSoil (e), (f).

Figure 4 .

 4 Figure 4. Regional carbon fluxes for 1960-2012. Total flux value for each SREX region in background, distribution of uncertainty between forcings, models and interactions as overlayed pieplots with size of the pies proportional to the total uncertainty value (equation (3) and tables S2-S4 in supplementary material).

Table 1 .

 1 Model description.

	Terrestrial biosphere	Community Land model	ISBA-CTRIP (Decharme
	model	CLM5 (Lawrence et al 2019)	et al 2019, Delire et al 2020)	JSBACH3.2 (Reick et al 2021)
	Parent climate model	Community Earth system	CNRM-ESM2-1 (Séférian	MPI-ESM1.2-LR (Mauritsen
		model CESM2	et al 2019)	et al 2019)
		(Danabasoglu et al 2020)
	Simulation	0.9 • latitude by 1.25 •	T127 (∼1.4 • )	T63 (∼1.9 • )
	resolution	longitude
	Fire	Process-based (Li et al 2012,	Simple, based on GlobFirm	Mechanistic, SPITFIRE
		Li and Lawrence 2017)	(Thonicke et al 2001)	(Lasslop et al 2014)
	Natural vegetation	No	No	Yes (Reick et al 2013)
	dynamics
	Nitrogen cycle	Revised (Fisher et al 2019)	No	New (Goll et al 2017)
	Plant hydrodynamics	(Kennedy et al 2019)	No	No
	Soil carbon	Yes	Yes	New (Goll et al 2015)
	decompostion
	Soil hydrology	Multilayer	Multilayer	Multilayer (Hagemann and
				Stacke 2015)
	Other	Numerous updates to	Similar	Similar
		hydrology, snow, gas
		exchange and crops

Table 2 .

 2 Atmospheric forcing description.

		Global soil wetness project
		(GSWP3) version 1.09 (Kim		Princeton (Sheffield
	Name	2017)	CRUJRA a (Harris 2019)	et al 2006) v2.2
	Time resolution	3 hourly	6 hourly	3 hourly
	Covered period	1901-2014	1901-2018	1901-2012
	Reanalysis used	20th Century Reanalysis version	Japanese Reanalysis (JRA) by	National Centers for
		2 (Compo et al 2011)	the Japanese Meteorological	Environmental Prediction
			Agency (JMA) b	(NCEP) reanalysis
	Corrected by
	air temperature	Climate Research Unit	CRU TS v4.03	CRU TS v3.2
		Timeseries (CRU TS) v3.21
	Precipitation	Global Precipitation	CRU TS v4.03	CRU, Global Precipitation
		Climatology Center (GPCC) v7		Climatology Project and
				TRMM Multi-satellite
				Precipitation Analysis
	Radiation	Surface Radiation Budget (SRB)	CRU TS v4.03	SRB
		datasets (downwelling radiation
		fluxes)
	Use	LS3MIP	GCP	Hydrology

a Erroneously named CRUNCEP in CMIP6 simulations. b http://hydro.iis.u-tokyo.ac.jp/GSWP3/. Project (LS3MIP, Van den Hurk et al 2016) of CMIP6. For the three models, land use transitions are forced by the Land Use Harmonization 2 (LUH2) time series

© 2022 The Author(s). Published by IOP Publishing Ltd

Acknowledgments

L H, C D, B D, R S, T S and R F acknowledge funding by the European Union's Horizon 2020 (H2020) research and innovation program under Grant Agreement No. 641816 (CRESCENDO) for L H, B D, C D and T S, No. 101003536 (ESM2025 -Earth System Models for the Future) for R S and No. 821003 (4C) for R F. C D and B D are supported by the 'Centre National de Recherches Météorologiques' (CNRM) of Météo-France and the 'Centre National de la Recherche Scientifique' (CNRS) of the French research ministry. D M L is supported by the National Center for Atmospheric Research, which is a major facility sponsored by the NSF under Cooperative Agreement 1852977 and by the U.S. Department of Energy, Office of Biological and Environmental Research Grant DE-FC03-97ER62402/A0101. J E M S N and T S would like to thank Veronika Gayler for conducting many of the CMIP6 MPI-ESM1.2-LR simulations. T S work was funded by the EU CRESCENDO Project. C D K acknowledges support by the Director, Office of Science, Office of Biological and Environmental Research of the U.S. Department of Energy under Contract DE-AC02-05CH11231 through the Early Career Research Program and the Regional and Global Model Analysis Program (RUBISCO SFA).

Data availability statement

The data that support the findings of this study are openly available at the following URL/DOI: https://esgf-node.llnl.gov/search/cmip6/.


Frédéric Ansgar Graf
email: ansgar.graf@universite-paris-saclay.fr

Frédéric Piéchon
email: frederic.piechon@universite-paris-saclay.fr

Piéchon Multifold Hopf

Hal

Multifold Hopf semimetals

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

Introduction.

Three-dimensional (3D) topological semimetals [START_REF] Armitage | Weyl and Dirac semimetals in three-dimensional solids[END_REF][START_REF] Lv | Experimental perspective on three-dimensional topological semimetals[END_REF] are materials with energy band crossing points that act as sources or sinks of Berry curvature [START_REF] Berry | Quantal phase factors accompanying adiabatic changes[END_REF][START_REF] Berry | Aspects of Degeneracy[END_REF][START_REF] Volovik | Zeros in the fermion spectrum in superfluid systems as diabolical points[END_REF][START_REF] Fang | The Anomalous Hall Effect and Magnetic Monopoles in Momentum Space[END_REF], so-called Berry monopoles [Fig. 1(a)]. The simplest example is a Weyl semimetal, with linear twoband crossings described by a Hamiltonian H W (q) = γ q • σ, where γ = ± is the chirality, q = (q x , q y , q z ) the momentum measured from the crossing, and σ a vector of Pauli matrices. The corresponding Berry curvature is given by

Ω α (q) = C α q 2|q| 3 , (1)
with α = ± labeling the two bands involved in the crossing and C α = -γα the Chern number (monopole charge).

Berry monopoles always come in pairs [Fig. 1(a)], which is a manifestation of the Nielsen-Ninomiya theorem [START_REF] Nielsen | A no-go theorem for regularizing chiral fermions[END_REF].

Each pair forms a dipole d 0 in the Brillouin zone, which lies at the heart of exotic phenomena such as Fermi arcs, anomalous Hall effect, and chiral anomaly [START_REF] Armitage | Weyl and Dirac semimetals in three-dimensional solids[END_REF][START_REF] Lv | Experimental perspective on three-dimensional topological semimetals[END_REF]. Any linear twoband (N = 2) crossing is of the Weyl form H W (q), and thus a Berry monopole. A common belief is that any linear multiband (N > 2) crossing also represents a Berry monopole. Indeed, Berry monopoles (1), with high Chern numbers C α , arise from multiband crossings governed by a generalized Weyl Hamiltonian, i.e. a pseudospin Hamiltonian H s (q) = γ q • S [START_REF] Bradlyn | Beyond Dirac and Weyl fermions: Unconventional quasiparticles in conventional crystals[END_REF][START_REF] Ezawa | Chiral anomaly enhancement and photoirradiation effects in multiband touching fermion systems[END_REF], where S is the (2s+1)-dimensional matrix representation of a pseudospin and α = -2s, ..., 2s.

The purpose of this paper is to draw attention to the existence of other linear multiband crossings, which are not of the Berry monopole type. In particular, we propose multifold Hopf semimetals (MHSs), with linear multifold crossing points [Fig. 1(c)-(e)] each of which carries a dipolar Berry curvature

Ω α (q) = κ α (d • q) q |q| 4 .
(

) 2
In contrast to the extended dipole d 0 of a Weyl semimetal, Eq. ( 2) describes a point-like/singular Berry dipole d that resides at a single band crossing point [Fig. 1(b)], with d representing an anisotropy axis but not a distance in momentum space. Note that singular Berry dipoles can also be obtained by merging a pair of Weyl nodes (d 0 → 0) or at topological phase transitions between twoband Hopf insulators [START_REF] Alexandradinata | Teleportation of Berry curvature on the surface of a Hopf insulator[END_REF][START_REF] Nelson | Delicate topology protected by rotation symmetry: Crystalline Hopf insulators and beyond[END_REF], but this always involves quadratic band touching. In contrast, for a MHS, the band crossings are linear. Moreover, the dipole charge κ α in Eq. ( 2), an integer for each band α, is symmetric with respect to zero energy, κ α = κ -α , while Ω α = -Ω -α for any twoband system and for any Berry monopole (1).

Our work illustrates three main points. First, linear multiband crossings can carry peculiar quantum geometric structures beyond Berry monopoles. We only discuss the Berry dipole case (2), but even Berry quadrupoles or octupoles are possible [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF]. Second, due to the Berry dipoles, MHSs have physical properties very different from those of known topological semimetals, despite the same low-energy spectrum. We show that Landau levels, anomalous Hall effect, and magnetoconductivity (all of which have been extensively studied in Weyl semimetals) exhibit distinct signatures of the Berry dipole (2), and it is clear that this extends to a host of other observables. Third, MHSs are governed by a delicate topology [START_REF] Nelson | Multicellularity of Delicate Topological Insulators[END_REF]. To demonstrate this, we propose simple lattice models for multiband (or N -band) Hopf insulators (MHIs) [START_REF] Lapierre | N -band Hopf insulator[END_REF] and show that MHSs appear at their topological phase transitions. This justifies dubbing them "Hopf" semimetals, extends notions of Refs. [START_REF] Alexandradinata | Teleportation of Berry curvature on the surface of a Hopf insulator[END_REF][START_REF] Nelson | Delicate topology protected by rotation symmetry: Crystalline Hopf insulators and beyond[END_REF] to the multiband case, and provides a fertile platform to study delicate topology.

Continuum models for MHSs. We start by introducing continuum models H ξ N (q) for Hopf semimetals, with ξ = ± a valley index. They are constructed to have a linear isotropic spectrum E α (q) = c α |q| with band velocity c α [Fig. 1(c)-(e)], identical to that of a pseudospin system H s (q), but a quantum geometry governed by Eq. ( 2). For a three-, four-, and fivefold Hopf semimetal, we choose

H ξ 3 (q) = 0 Q ξ 3 (Q ξ 3 ) † 0 2 , Q ξ 3 = q ξ --iq z , (3a)
H ξ 4 (q) = 0 2 Q ξ 4 (Q ξ 4 ) † 0 2 , Q ξ 4 = aq ξ -iaq z ibq z bq ξ + , (3b)
H ξ 5 (q) = 0 3 Q ξ 5 (Q ξ 5 ) † 0 2 , Q ξ 5 =    0 i √ 2q z iq z q ξ + √ 2q ξ + 0    , (3c)
respectively, where q ξ ± ≡ ξq x ± iq y , and a > b > 0. The band velocities are c α = 0, ±1 for the model (3a), c α = ±a, ±b for (3b), and c α = 0, ±1, ± √ 2 for (3c). The terms ∼ q ξ ± are familiar from graphene and Weyl semimetals, and the terms ∼ q z provide a third direction fixing the dipole axis as d = (0, 0, ξ) (however it is possible to construct models with tunable d [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF]). Computing the Berry curvature [START_REF] Graf | Berry curvature and quantum metric in N -band systems: An eigenprojector approach[END_REF], one finds κ α = -1, 2, -1 for (3a), κ α = -1, 1, 1, -1 for (3b), and κ α = -3, 1, 4, 1, -3 for (3c) from lowest to highest band [Fig. 1(c)-(e)], with large dipole charges carried by the flat bands.

The models (3) have two important symmetries, linked to the dipole charges and the dipole axis, respectively. Namely, first, a chiral symmetry S -1 H ξ N (q)S = -H ξ N (q) with a diagonal matrix S and S 2 = 1. Second, an axial rotation symmetry

[L d + Σ d , H ξ N (q)] = 0, with L d = d•L the projection of angular momentum L = -i(q × ∇ q )1,
and Σ d a diagonal matrix acting as an effective spin projection [START_REF]More precisely, we have Σ d = 1 diag[END_REF]. In contrast, a pseudospin system (Berry monopole) H s (q) has a charge-conjugation parity symmetry C -1 H s (q)C = -H * s (q) with C = exp(iπS y ), and a full rotation symmetry [L + S, H s (q)] = 0. The chiral and axial rotation symmetries determine general properties of the physical responses studied in the following.

Physical properties of MHSs -continuum perspective. We now illustrate the strong impact of the Berry dipole on Landau levels, anomalous Hall conductivity, and magnetoconductivity. We first consider these effects for a sin- gle band crossing H ξ N (q). Below, we will present tightbinding models for MHSs, with one or more Berry dipoles in the Brillouin zone, and discuss how the effects translate to the lattice setup.

(a) Landau levels (LLs). Consider Eq. ( 3) for a strong magnetic field B = B B, where B = (0, sin θ, cos θ) without loss of generality due to axial rotation symmetry. The LLs form a 1D dispersion in terms of q 0 = B • q, and are particle-hole symmetric due to chiral symmetry. The exact LLs for a threeband crossing (3a) are [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF] n,ξ

α = c α 2eB n + 1 -κ α ξ cos θ 2 + q 2 0 , (4)
with c α = 0, ±1 and n ∈ {0, 1, ...} the LL index. The flat band is maintained and the dispersive bands carry a clear signature of the Berry dipole's charge (κ α ) and orientation (ξ cos θ ≡ B • d). As a consequence, the LL spectrum strongly depends on the magnetic field direction: it is gapped for B d, and gapless for B d, see Fig. 2(a). This is a pure quantum geometric effect, since the zero-field spectrum is fully isotropic. Similarly, we find the exact LLs of a fourband crossing (3b):

n,ξ α = α 1 η + + η -+ α 2 (η + -η -) 2 + ν 2 , η ± = c 2 ± 2 2eB n + 1 -κ ± ξ cos θ 2 + q 2 0 . (5)
Here, we use a band index tuple α = (α 1 , α 2 ) with α i = ±; c + = a and c -= b are the band velocities of the two cones, and κ ± = ∓1 the corresponding dipole charges; moreover, ν = c + c -e(B × d) x = abeBξ sin θ. Again, the Berry dipole (2) explicitly appears and the LLs can be tuned by rotating B, see Fig. 2(b), however the spectrum is now always gapless. Finally, the LLs of the fivefold crossing (3c) are similar to the LLs (4), see [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF]. Note the big difference with the LLs of a pseudospin (Berry monopole) H s (q). Those are independent of B due to full rotation symmetry, exhibit a number of chiral LLs determined by the Chern number [START_REF] Bradlyn | Beyond Dirac and Weyl fermions: Unconventional quasiparticles in conventional crystals[END_REF][START_REF] Ezawa | Chiral anomaly enhancement and photoirradiation effects in multiband touching fermion systems[END_REF][START_REF] Delplace | Berry-Chern monopoles and spectral flows[END_REF], and any flat band is destroyed since B breaks parity symmetry [START_REF] Similarly | the LLs of Dirac fermions HD(q) are independent of B. However they are particle-hole symmetric due to chiral symmetry. Indeed, for a = b, where Eq. (3b) becomes a Dirac semimetal HD(q) = q • Γ with anticommuting matrices Γx,y,z, Eq. (5)[END_REF]. One can also compare to the LLs of an extended Berry dipole d 0 formed from two Weyl nodes. Those obviously depend on B [START_REF] Saykin | Landau levels with magnetic tunneling in a Weyl semimetal and magnetoconductance of a ballistic p-n junction[END_REF], but this dependence is quite distinct from the one of Eqs. ( 4)- [START_REF] Volovik | Zeros in the fermion spectrum in superfluid systems as diabolical points[END_REF].

It is also interesting to analyze the LLs semiclassically. Indeed, Eq. ( 4) can be easily obtained using Onsager quantization [START_REF] Onsager | Interpretation of the de Haas-van Alphen effect[END_REF] plus intraband (Berry curvature and orbital magnetic moment) corrections [START_REF] Roth | Semiclassical Theory of Magnetic Energy Levels and Magnetic Susceptibility of Bloch Electrons[END_REF][START_REF] Mikitik | Manifestation of Berry's Phase in Metal Physics[END_REF][START_REF] Fuchs | Topological Berry phase and semiclassical quantization of cyclotron orbits for two dimensional electrons in coupled band models[END_REF][START_REF] Gao | Zero-field magnetic response functions in Landau levels[END_REF]. More importantly, the origin of ν ∼ |B × d| in Eq. ( 5) can be explained semiclassically via interband coupling of degenerate orbits originating from the two cones [Fig. 1(d)]. This can be seen from an original approach to Landau quantization of degenerate orbits developed in [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF].

(b) Anomalous Hall effect and magnetotransport. Consider now a crossing H ξ N (q) in the presence of weak electric and magnetic fields E and B. We adopt a standard approach to describe the linear response of the system, by solving the semiclassical Boltzmann equation in the relaxation time (τ ) approximation to first order in E and B [START_REF] Ziman | Electrons and Phonons: The Theory of Transport Phenomena in Solids[END_REF]. Taking into account corrections due to Berry curvature and orbital magnetic moment [START_REF] Xiao | Berry phase effects on electronic properties[END_REF], and working at zero temperature, we find several nontrivial effects [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF]:

(i) At zeroth order in B, a single crossing point H ξ N (q) causes a non-dissipative anomalous Hall (AH) current

j AH = σ AH E × d. ( 6
)
This current is similar in form to the AH current j W AH = σ W AH E × d 0 created by a pair of coupled Weyl nodes [START_REF] Klinkhamer | Emergent CPT violation from the splitting of Fermi points[END_REF][START_REF] Burkov | Weyl Semimetal in a Topological Insulator Multilayer[END_REF][START_REF] Yang | Quantum Hall effects in a Weyl semimetal: Possible application in pyrochlore iridates[END_REF] (here d 0 may represent the distance between Weyl nodes or a tilt). Importantly, however, it is opposite in parity: σ AH in Eq. ( 6) is odd in the Fermi level E F , while σ W AH is even. This opposite parity is a direct consequence of the symmetry of the dipole charges (κ α = κ -α ), as proven in [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF].

(ii) At first order in B, there is a dissipative quantum geometric current j geo (B). For example, for model (3a),

j geo (B) = - e 3 τ 12π 2 [(E • B)d + (E • d)B + (B • d)E] , (7)
and similarly for the four-and fivefold crossings (3b), (3c). Just like the AH current, the current ( 7) is similar in form to well-known effects from a pair of coupled Weyl nodes: a current j CA ∼ (E • B)d 0 induced by the chiral anomaly [START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF][START_REF] Zyuzin | Magnetotransport of Weyl semimetals due to the chiral anomaly[END_REF][START_REF] Sharma | Chiral anomaly and longitudinal magnetotransport in type-II Weyl semimetals[END_REF]; a current j CME ∼ δ B due to the chiral magnetic effect [START_REF] Zyuzin | Magnetotransport of Weyl semimetals due to the chiral anomaly[END_REF][START_REF] Vilenkin | Equilibrium parity-violating current in a magnetic field[END_REF][START_REF] Fukushima | Chiral magnetic effect[END_REF][START_REF] Zyuzin | Weyl semimetal with broken time reversal and inversion symmetries[END_REF], where δ is an energy difference between the valleys; a current j MCE ∼ (B • d 0 )E caused by the magnetochiral effect [36, 37]. However, j geo (B) again has opposite parity: it is even in E F , whereas the currents j CA , j CME , j MCE are odd. While this cannot be seen from the continuum formula [START_REF] Nielsen | A no-go theorem for regularizing chiral fermions[END_REF], it can be shown to be true by general symmetry considerations [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF]. Tight-binding models for MHSs. We now move from the continuum to a lattice description, introducing two different classes of tight-binding models for MHSs (Fig. 3) both of which have a low-energy theory (3).

(a) Valley-MHSs. Those have an even number of Berry dipole crossings in the Brillouin zone, with opposite dipole orientation in each valley [Fig. 3(a)]. We consider a version with preserved [h N (k)] and broken [ hN (k)] time-reversal symmetry, respectively:

h N (k) = 0 Q N Q † N 0 , hN (k) = 0 QN Q † N 0 , Q 3 = w k -is z , Q3 = s --is z , Q 4 = aw k ias z ibs z bw * k , Q4 = as -ias z ibs z bs + , Q 5 =    0 i √ 2s z is z w * k √ 2w * k 0    , Q5 =    0 i √ 2s z is z s + √ 2s + 0    .
(8) Here, w k ≡ 2 3 j exp(ik • δ j ), where δ 1,2 = 1 2 (± √ 3, 1, 0) and δ 3 = (0, -1, 0), and we use shorthand notations s i ≡ sin k i and s ± ≡ s x ± is y . The Hamiltonians h N (k) [ hN (k)] represent nearest-neighbor tight-binding models on a hexagonal [tetragonal] lattice, consisting of 2D hon-eycomb [square] layers stacked in a particular way along ẑ. As an example, the case N = 3 is shown in Fig. 3(c) and more details can be found in [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF]. Note that the Bravais period is doubled along ẑ such that nodal points can appear only in the k z = 0 plane. Indeed, for h N (k) the Brillouin zone contains one N -fold nodal point at the K (ξ = +) and one at the K' (ξ = -) valley, where

K ξ = -ξ4π/(3 √ 3)
x. Their low-energy theory is given by H + N (q) (Berry dipole up) and H - N (q) (Berry dipole down), respectively. For hN (k) there are also two N -fold nodal points in the Brillouin zone, located at the Γ and M point and described by H + N (q) and H - N (q), respectively. (b) True MHSs. Those have an odd number of Berry dipoles in the Brillouin zone [Fig. 3(b)]. We consider

h N (k) ≡ hN (k) + g(k)D N , g(k) = ∆ 0 + cos k x + cos k y + cos k z , (9)
where D 3 = antidiag(1, 0, 1), D 4 = antidiag(-a, b, b, -a), and

D 5 = antidiag( √ 2, 1, 0, 1, √ 2
). Such models can be built on a simple cubic lattice with N orbitals per site, see Fig. 3(d) for the case N = 3. To have a semimetal, we allow the parameter ∆ 0 to take only discrete values ∆ 0 = ±1, ±3. Indeed, for ∆ 0 = -3 (∆ 0 = 3) there is one nodal point at the Γ (R) point of the Brillouin zone, and for ∆ 0 = -1 (∆ 0 = 1) there are three nodal points at the inequivalent X (M) points. The low-energy theory of all these points is of the form (3) [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF].

Physical properties of MHSs -lattice perspective. For a Fermi level close to the nodal points, the physical properties of the models (8) and ( 9) are simply obtained by summing the continuum results over all valleys. It is then clear from Fig. 3(a)-(b) that the currents ( 6) and ( 7) cancel for a valley-MHS, but not for a true MHS. To confirm this, we have numerically calculated the anomalous Hall conductivity of h N (k) for parameters such that a single Berry dipole crossing exists in the Brillouin zone [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF]. The N = 3 case is shown in Fig. 3(e). This calculation confirms the previous claim that j AH is odd in E F , while it is even for a pair of Berry monopoles. A similar lattice calculation (not shown) confirms the existence of weak-field magnetocurrents j geo (B) for a Hopf semimetal [START_REF] Ezawa | Chiral anomaly enhancement and photoirradiation effects in multiband touching fermion systems[END_REF], which are even in E F . In contrast to the currents, the LLs are nontrivial also for a valley-MHS. For example, consider h 3 (k) [Fig. 3(c)] under strong B. For θ close to 0 or π, it is clear from Fig. 2(a) that one of the K and K' valleys will be gapped and the other gapless, such that only one valley polarization is accessible in the ultraquantum limit (E F ∼ 1/l B ). This could be interesting, for example, for magnetooptics.

Relation to multiband Hopf insulators. In the MHS Hamiltonian [START_REF] Ezawa | Chiral anomaly enhancement and photoirradiation effects in multiband touching fermion systems[END_REF], we now replace the discrete values ∆ 0 by a continuous parameter ∆. The Hamiltonians

h Hopf N (k) ≡ h N (k, ∆ = ∆ 0 ) (10)
are nearest-neighbor tight-binding models for multiband Hopf insulators (MHIs) [START_REF] Lapierre | N -band Hopf insulator[END_REF], a generalization of wellknown twoband Hopf insulators [START_REF] Moore | Topological Surface States in Three-Dimensional Magnetic Insulators[END_REF][START_REF] Deng | Hopf insulators and their topologically protected surface states[END_REF]. To confirm this, we compute the Hopf number

N Hopf = W 3 [U (k)] = BZ d 3 k 24π 2 χ(k) (11)
as the third winding number W 3 of the unitary matrix U (k) that diagonalizes h Hopf N (k), which can be written as an integral of the Hopf density χ(k) ≡ ijl Tr[u i (k)u j (k)u l (k)] over the Brillouin zone, with u i (k) ≡ U † (k)∂ i U (k) and ∂ i ≡ ∂/∂k i . From Eq. ( 11), we find a quantized Hopf number for ∆ = ∆ 0 [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF], as shown in Fig. 4. At the topological phase transitions (TPTs) ∆ = ∆ 0 , all N -1 band gaps close simultaneously at one or three points of the Brillouin zone. Clearly, the MHSs ( 9) are critical points of the MHIs [START_REF] Alexandradinata | Teleportation of Berry curvature on the surface of a Hopf insulator[END_REF], in analogy to how quadratic band touchings with Berry dipole mediate TPTs between twoband Hopf insulators [START_REF] Alexandradinata | Teleportation of Berry curvature on the surface of a Hopf insulator[END_REF][START_REF] Nelson | Delicate topology protected by rotation symmetry: Crystalline Hopf insulators and beyond[END_REF].

Three remarks are in order. First, h Hopf 3 (k) is actually well known as a chiral topological insulator (TI) in the literature [START_REF] Wang | Probe of Three-Dimensional Chiral Topological Insulators in an Optical Lattice[END_REF][START_REF] Lian | Machine Learning Topological Phases with a Solid-State Quantum Simulator[END_REF]; it was introduced for constructing fractional topological phases [START_REF] Neupert | Noncommutative geometry for threedimensional topological insulators[END_REF], and we identify it as a threeband HI [START_REF]Note in particular that, when chiral symmetry is broken without a[END_REF]. Second, the jump δN Hopf at a TPT can be understood by using a continuum Hopf number and summing over all band crossing points [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF]. Third, we can define a valley-MHI as h v-Hopf N (k) = hN (k) + ∆D N , which has χ(k) = 0 but N Hopf = 0.

Experimental realization. Given the simplicity of the MHS models, their existence in crystals appears possible. Berry dipole crossings could be probed via the peculiar, strongly valley-dependent depence of the Landau levels on the magnetic field orientation [Fig. 2(a)]. AH currents and linear magnetocurrents induced by the Berry dipoles have symmetries opposite to those of a pair of Weyl nodes [Fig. 3(e)], which could be probed by varying the Fermi level close to half filling. A different route involves artificial systems such as ultracold atoms, photonic crystals, or superconducting circuits. Those have been used to realize semimetallic phases with two-and multifold crossings [START_REF] Zhang | Topological quantum matter with cold atoms[END_REF][START_REF] Zhu | Emergent pseudospin-1 Maxwell fermions with a threefold degeneracy in optical lattices[END_REF][START_REF] Fulga | Geometrically protected triple-point crossings in an optical lattice[END_REF][START_REF] Hu | Topological Triply Degenerate Points Induced by Spin-Tensor-Momentum Couplings[END_REF][START_REF] Lu | Experimental observation of Weyl points[END_REF][START_REF] Chen | Photonic crystals possessing multiple Weyl points and the experimental observation of robust surface states[END_REF][START_REF] Riwar | Multi-terminal Josephson junctions as topological matter[END_REF][START_REF] Wang | Optical Interface States Protected by Synthetic Weyl Points[END_REF][START_REF] Tan | Topological Maxwell Metal Bands in a Superconducting Qutrit[END_REF], and a threeband crossing [START_REF] Palumbo | Revealing Tensor Monopoles through Quantum-Metric Measurements[END_REF] similar to Eq. (3a) was recently observed using a transmon in a cavity [START_REF] Tan | Experimental Observation of Tensor Monopoles with a Superconducting Qudit[END_REF]. Moreover, various schemes for observing Hopf numbers in insulators have been proposed [START_REF] Deng | Probe Knots and Hopf Insulators with Ultracold Atoms[END_REF][START_REF] Ünal | Hopf characterization of two-dimensional Floquet topological insulators[END_REF][START_REF] Schuster | Realizing Hopf Insulators in Dipolar Spin Systems[END_REF][START_REF] Schuster | Floquet engineering ultracold polar molecules to simulate topological insulators[END_REF] and could be extended to the MHIs [START_REF] Alexandradinata | Teleportation of Berry curvature on the surface of a Hopf insulator[END_REF]

HI h Hopf

3

(k) (viewed as a chiral TI) there already exists not only a proposal based on ultracold atoms [START_REF] Wang | Probe of Three-Dimensional Chiral Topological Insulators in an Optical Lattice[END_REF], but also an experimental realization based on machine learning analysis of a nitrogen-vacancy center in diamond [START_REF] Lian | Machine Learning Topological Phases with a Solid-State Quantum Simulator[END_REF]. The MHS ( 9) is reached at the topological phase transitions of such experiments. While it appears challenging to implement the magnetic responses discussed above in artificial systems, we expect the Berry dipoles to also affect more accessible phenomena such as interferometry of atomic wave packets or quench dynamics [START_REF] Zhang | Topological quantum matter with cold atoms[END_REF].

Conclusions. We proposed multifold Hopf semimetals (MHSs), featuring N -fold linear band crossings with Berry dipole, distinct from the Berry monopole crossings of known topological semimetals [START_REF] Armitage | Weyl and Dirac semimetals in three-dimensional solids[END_REF][START_REF] Lv | Experimental perspective on three-dimensional topological semimetals[END_REF]. The Berry dipole strongly affects observables, e.g. Landau levels, anomalous Hall effect, and magnetoconductivities. Gapping out MHSs, we found the first simple lattice models for multiband Hopf insulators [START_REF] Lapierre | N -band Hopf insulator[END_REF], avoiding the second-neighbor hoppings of twoband Hopf insulators [START_REF] Moore | Topological Surface States in Three-Dimensional Magnetic Insulators[END_REF][START_REF] Deng | Hopf insulators and their topologically protected surface states[END_REF] and providing an ideal platform to study delicate topology [START_REF] Nelson | Multicellularity of Delicate Topological Insulators[END_REF].

To close, we mention interesting perspectives. First, the systems introduced here are, in a sense, 3D analogs of 2D Dirac semimetals and Chern insulators, see Table I and [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF]. To make such an analogy more precise, one could systematically analyze all possible ways to perturb MHSs and establish the corresponding phase diagrams, similar to recent work on the conversion between Weyl points, nodal lines, quadratic Berry dipole touchings and twoband Hopf insulators [START_REF] Nelson | Delicate topology protected by rotation symmetry: Crystalline Hopf insulators and beyond[END_REF][START_REF] Liu | Symmetry-protected topological Hopf insulator and its generalizations[END_REF][START_REF] Sun | Conversion Rules for Weyl Points and Nodal Lines in Topological Media[END_REF][START_REF] Bouhon | Non-Abelian reciprocal braiding of Weyl points and its manifestation in ZrTe[END_REF]. Second, the Hamiltonians (3a)-(3b) seem to be related to 4D semimetals with tensor monopoles [START_REF] Palumbo | Revealing Tensor Monopoles through Quantum-Metric Measurements[END_REF][START_REF] Zhu | Fourdimensional semimetals with tensor monopoles: From surface states to topological responses[END_REF] by dimensional reduction. Finally, the models with Berry quadrupole and octupole provided in [START_REF]See Supplemental Material for: (i) Continuum models with Berry quadrupole and octupole structure. (ii) A MHS continuum model with a tunable Berry dipole vector. (iii) A derivation of the exact LL spectrum for MHS continuum models. (iv) An original semiclassical approach to LLs in MHSs, going beyond Onsager quantization and taking into account interband coupling between degenerate orbits. (v) A review of semiclassical magnetotransport in multiband systems. This is used to compute anomalous Hall and magnetoconductivities of MHS continuum and lattice models, and to explain their symmetries with respect to EF . (vi) Details on the lattice models for MHSs defined by Eqs. (8)[END_REF] indicate the possibility to establish a full hierarchy of Berry multipole crossings.

FIG. 1 .

 1 FIG. 1. (a) Known topological semimetals are based on linear crossings each of which is a Berry monopole (left). Those appear in monopole-antimonopole pairs (right). (b) Multifold Hopf semimetals are based on linear crossings each of which is a Berry dipole. (c)-(e) Energy spectrum and dipole charges of the models (3a)-(3c).

FIG. 2 .

 2 FIG. 2. (a) LL spectrum (4) for ξ = + for three different θ, in units of the inverse magnetic length 1/lB ≡ √ eB. It can be tuned from gapped to gapless by rotating B. (b) LL spectrum (5) for ξ = +, a/2 = b = 1. It is gapless for any θ.

FIG. 3 .

 3 FIG. 3. (a)-(b) Valley-MHSs (true MHSs) have an even (odd) number of Berry dipoles. (c)-(d) Lattice models for threeband HSs, with lattice constant a0 = 1. (e) Red: AH conductivity for h3(k) with ∆0 = -3, caused by the Berry dipole at the Γ point. Black: Typical AH conductivity for a Weyl semimetal with a single pair of Weyl points [12].

FIG. 4 .

 4 FIG. 4. (a) Hopf invariant[START_REF] Nelson | Delicate topology protected by rotation symmetry: Crystalline Hopf insulators and beyond[END_REF] for the models[START_REF] Alexandradinata | Teleportation of Berry curvature on the surface of a Hopf insulator[END_REF]. Topological phase transitions occur for ∆ = ∆0, where gaps close at the Γ, X, M or R points of the Brillouin zone, see (b). Red (blue) colors denote a positive (negative) sign of δN Hopf .

  . Most notably, for the threeband

	valley-DS ↔ valley-MHS	valley-CI ↔ valley-MHI
	n even, top. dens. singular	top. dens. trivial, no topology
	graphene ↔ h N	Semenoff insulator [59] ↔ h v-Hopf N
	top. DS ↔ MHS	CI ↔ MHI
	n odd, top. dens. trivial	top. dens. nontrivial, topology
	Haldane boundary [60] ↔ h N	Haldane insulator [60] ↔ h Hopf N
	TABLE I. Analogy between 2D systems ↔ 3D systems,
	based on the number n of linear crossing points and the
	topological density (Berry curvature ↔ Hopf density). The
	third line of each box gives an example.	DS = Dirac
	semimetal, MHS = Multifold Hopf semimetal, CI = Chern in-
	sulator, MHI = Multiband Hopf insulator.

Acknowledgements. We are indebted to A. Mesaros for very enlightening discussions on Hopf insulators.