Vicente García-Navas
email: vicente.garcianavas@gmail.com.

Wilfred Thuiller

Farmland bird assemblages exhibit higher functional and phylogenetic diversity than forest assemblages in France

Keywords: avifauna, agriculture, biodiversity, community dynamics, functional traits Running head: Diversity dynamics in land bird communities

Aim: Under a global change scenario, research focused on changes in assembly patterns over spatial and temporal axes is more than timely as it will improve our understanding about how biological communities respond to anthropogenic disturbance. Despite an increasing need to assess whether associations among diversity and community metrics change in relation to environmental heterogeneity, the way in which community assembly rules vary across habitats has been hardly explored. Here, we tested for differences in patterns of functional diversity (FD) and phylogenetic diversity (PD) between farmland and forest bird communities.

Taxon: 107 species of breeding common birds.

Location: continental France.

Methods:

We used an extensive dataset (13 years; 7,115 bird communities) from the French Breeding Bird survey in conjunction with a matrix of 142 functional traits (including information on habitat, diet, life-stories, behaviour, and morphology) to compute different metrics of FD and PD.

Results: We found that farmland assemblages showed higher FD and PD than forest assemblages, which were phylogenetically clustered. Both FD and PD of forest assemblages increased with increasing species richness, whereas in farmland assemblages the relationship turned out to be asymptotic in both cases. It may be due to the accumulation of generalists, which can end up displacing specialist species when the environment becomes oversaturated triggering a decline in diversity. Contrary to expectations, FD and PD of farmland assemblages increased over the study period, whereas forest assemblages showed a non-linear pattern.

Farmland and forest assemblages also showed divergent trajectories over time in relation to FD metrics.

Main conclusions:

We conclude that, although farmland intensification has led to a sharp decline in population of farmland birds, agriculture landscapes in Southern Europe still harbor diversity-rich communities probably due to the legacy effects of past land-use (traditional practices). Our study highlights the need to take into account the influence of historic landscape configurations when assessing the effect that contemporary land uses have on biotic communities.

Introduction

Community phylogenetics seeks to characterize the structure of communities and unravel the relative importance of chance (stochastic processes) and ecological rules (deterministic processes) in shaping community composition (reviewed in [START_REF] Pausas | The jungle of methods for evaluating phenotypic and phylogenetic structure of communities[END_REF].

Specifically, at a local scale, the quantification of functional and phylogenetic structures of assemblages allows recognition of factors leading to competitive exclusion or to the coexistence of ecologically similar taxa (environmental filtering vs. limiting similarity: e.g. [START_REF] Stevens | Phylogenetic structure illuminates the mechanistic role of environmental heterogeneity in community organization[END_REF]). Thus, examining both biodiversity dimensions -functional and phylogenetic-turns out to be necessary for understanding the complete composition, structure and dynamics of communities (Devictor et al. 2010a).

Although several studies have assessed the phylogenetic or functional structure of bird assemblages (most of them focused on specific bird families or clades: e.g. [START_REF] Graham | Phylogenetic structure in tropical hummingbird communities[END_REF][START_REF] Gómez | A phylogenetic approach to disentangling the role of competition and habitat filtering in community assembly of Neotropical forest birds[END_REF], there is a paucity of studies examining the relationship between the phylogenetic and functional structure of assemblages and the underlying assembly processes in contrasting habitats [START_REF] Thuiller | Predicting global change impacts on plant species distributions: future challenges[END_REF]; but see [START_REF] Sobral | Functional and phylogenetic structure of forest and savanna bird assemblages across spatial scales[END_REF][START_REF] Morelli | Evidence of evolutionary homogenization of bird communities in urban environments across Europe[END_REF]. This is noteworthy because it's expected that species pools and the structure of assemblages in, for example, forests and grasslands (the two extremes along a close-open continuum) are shaped by distinct ecological, evolutionary and historical mechanisms. In addition to spatial heterogeneity, the existence of temporal variability in habitat quality constitutes an important factor to consider when addressing the response of biological communities to disturbance (e.g. [START_REF] Méndez | Functional diversity across space and time: trends in wader communities on British estuaries[END_REF]. While some habitats have undergone significant modifications during the last decades due to drastic changes in land use, some others remain unaltered and thus, they may act as true reservoirs of biodiversity. Specifically, agriculture landscapes have changed at a great pace during the last few decades due to farming intensification [START_REF] Tscharntke | Landscape perspectives on agricultural intensification and biodiversity ecosystem service management[END_REF]. New agricultural practices have led to habitat simplification through the expansion of monocrops, the disappearance of structural elements (hedgerows, tree lines, etc.) and the abandonment of the once common practice of letting fields lie fallow. Each one of these factors seems to have contributed to the observed decline in farmland bird populations [START_REF] Benton | Farmland biodiversity: is habitat heterogeneity the key?[END_REF][START_REF] Donald | Further evidence of continentwide impacts of agricultural intensification on European farmland birds, 1990-2000[END_REF]. Besides species richness and abundance, agricultural intensification may have an impact on other functional and structural properties of biological communities. For instance, the progressive homogenization of habitat structure may decrease the proportion of specialist species [START_REF] Julliard | Spatial segregation of specialists and generalists in bird communities[END_REF]. High Nature Value (HNV) farming systems, characterized by low intensity farming and diverse, small-scale mosaics of land-use types have been shown to be effective in reverting these negative effects on biodiversity (e.g. [START_REF] Morelli | Relative importance of marginal vegetation (shrubs, hedgerows, isolated trees) surrogate of HNV farmland for bird species distribution in Central Italy[END_REF][START_REF] Aue | How High Nature Value (HNV) farmland is related to bird diversity in agro-ecosystems -Towards a versatile tool for biodiversity monitoring and conservation planning[END_REF].

Here, we take advantage of long-term breeding bird dataset on the French avifauna to investigate the processes that govern community assembly in different forest types and farmland habitats managed with a different degree of intensification. Specifically, we analyze spatial and temporal variation in the structure of local communities in order to address several hypotheses.

Firstly, in terms of phylogenetic structure, it has been argued that forests contain a greater proportion of evolutionarily distinctive species compared to most human-modified habitats due to its more ancient origin [START_REF] Frishkoff | Loss of avian phylogenetic diversity in neotropical agricultural systems[END_REF]. Furthermore, forests are more stable habitats, favoring the colonization and successful establishment of many clades [START_REF] Dreiss | Taxonomic, functional, and phylogenetic dimensions of rodent biodiversity along an extensive tropical elevational gradient[END_REF]. On the other hand, farmland habitats represent more unpredictable and stressful conditions, which may result in physiological constraints that restrict the presence of some clades in these environments [START_REF] Clavel | Worldwide decline of specialist species: toward a global functional homogenization?[END_REF]. Accordingly, we hypothesize that: (1) farmland bird assemblages posse a lower phylogenetic diversity and are phylogenetically clustered, and (2) forest bird assemblages show a pattern of phylogenetic overdispersion.

Secondly, it is expected that the availability and distribution of resources largely determines the functional structure of communities [START_REF] Weiher | Assembly rules, null models, and trait dis persion: new questions from old patterns[END_REF]. Under harsh or stressful conditions, environmental filtering operating on convergent traits may generate functional clustering (Cavender-Bares et al. 2009), whereas in extremely productive habitats in which species primarily face adversity from competition, limiting similarity may favor the coexistence of species with different functional traits [START_REF] Weiher | Community assembly rules, morphological dispersion, and the coexistence of plant species[END_REF][START_REF] Gotelli | Macroecological signals of species interactions in the Danish avifauna[END_REF].

Given that forests should provide a greater availability of niches and food resources than human-modified landscapes like those resulting from agricultural practices, we predict: (3) a stronger effect of competition in forest sites resulting in functionally overdispersed assemblages; and (4) a prevalence of functional clustering in farmland sites due to the higher environmental harshness.

Thirdly, because the probability of adding new evolutionary lineages/functional groups to an assemblage decreases as assemblages become more taxonomically diverse, phylogenetic/functional diversity is expected to be positively associated with species richness in a decelerating fashion [START_REF] Cumming | Contrasting spatial patterns of taxonomic and functional richness offer insights into potential loss of ecosystem services[END_REF][START_REF] Kluge | Phylogenetic diversity, trait diversity and niches: species assembly of ferns along a tropical elevational gradient[END_REF]. However, the strength of this association is expected to vary depending on the evolutionary heritage of communities and the evolutionary distinctiveness or 'originality' of species conforming a given assemblage [START_REF] Pavoine | Is the originality of a species measurable?[END_REF]. As forests represent the potential natural vegetation in the temperate zone and thus, they are considered older environments than more anthropic habitats, we expect: (5) a steeper relationship (i.e. diversity level-off faster) between taxonomic and phylogenetic/functional diversity in forest assemblages in comparison with farmland assemblages.

Finally, temporal trends in diversity at the community level are expected to mirror environmental changes and the anthropogenic influence on biodiversity [START_REF] Magurran | Long-term datasets in biodiversity research and monitoring: assessing change in ecological communities through time[END_REF].

In France, farmland bird populations have fallen by one third in the last two decades [START_REF] Jiguet | Suivi Temporel des Oiseaux Communs. Bilan 2007 pour la France[END_REF]. Reduction in effective numbers may be accompanied by a depauperation of farming bird assemblages. Meanwhile, large-scale forest maturation and afforestation have augmented in most of the Southern Europe due to increased abandonment of rural land during the last century and the restriction of agricultural practices to highly productive areas. Consequently, species associated to forested areas seem to show a positive trend (e.g. [START_REF] Gil-Tena | Effects of forest landscape change and management on the range expansion of forest bird species in the Mediterranean region[END_REF]. Based on this, we predict (6) the existence of differences in temporal trends between forest and farmland assemblages.

Material and Methods

Bird data

We used abundance data on breeding birds from the French Breeding Bird Survey (Suivi Temporel des Oiseaux Communs, STOC), a standardized monitoring program launched in 2001, in which skilled volunteer ornithologists identify breeding birds by song or visual contact each spring. Briefly, a 2×2 km plot consisting of 10-point counts -separated by at least 300 m-is randomly assigned to each volunteer. Volunteers record year after year all birds seen or heard in each of the ten sites during a 5-min period. Point counts are visited at dawn twice within three weeks around mid-May to ensure the detection of both early-and late-breeders. The maximum value of the two censuses is retained for each species as yearly abundance index. More details about the scheme and census protocols are given in [START_REF] Monnet | Asynchrony of taxonomic, functional and phylogenetic diversity in birds[END_REF].

We analyzed data spanning a 13-yr period, from 2002 to 2014. Since there exist large differences among the study plots in terms of sampled years, we only included those plots located in farmland or forest habitats that were monitored for more than 5 years in order to avoid potential biases associated to variation in sampling effort. Forest habitats comprised broad-leaved, coniferous and mixed forests, whereas farmland plots were located in areas in which predominates one of the following uses according to the CORINE Land Cover inventory: non-irrigated arable land, complex cultivation patterns, agricultural land with significant areas of natural vegetation and pastures or permanent grasslands under agricultural use. In total, we gathered records from a total of 812 sites, of which 279 were located in forest habitats and 533 in agricultural land. Species with a large home range (raptors) and species that only breed in wetlands were excluded from the analysis because they are not properly monitored by the program. We also excluded species that were rarely observed and whose presence cannot be considered regular, but anecdotal (i.e. rare species). Our final data set included 107 species, which account for 99% of the total number of counted individuals (excluding raptors and water birds).

Functional traits

A matrix of 142 traits was compiled for the 107 bird species using an extended version of the dataset described in [START_REF] Thuiller | The European functional tree of bird life in face of global change[END_REF]. This extensive dataset includes information on habitat, diet, life-stories, behaviour, morphology and reproduction providing a comprehensive characterization of the ecology of each species and how they interact with the environment and with other organisms (Table S1). Traits were divided into seven categories, which included both multi-trait and single-trait variables. Although we used the whole set of variables to quantify functional diversity (hereafter 'generic functional diversity') we also computed functional diversity indices separately for habitat and feeding traits (see more below). We paid special attention to these two subsets because they clearly asses different aspects of the niche, which defines the level of specialization of a given species (Devictor et al. 2010b). Habitat traits reflect the environmental requirements of each species and thus, the extent to which species overlap in their use of space (Grinnellian or scenopoetic niche). Feeding traits reflect the trophic position and the guild to which each species belongs to. That is, its place in the biotic environment and the way in which each species interacts with food and other species (Eltonian or bionomic niche). These groups of traits correspond to different dimensions of the niche and thus, do not necessarily should exhibit congruent patterns over space or time [START_REF] Cisneros | Season-specific and guild-specific effects of anthropogenic landscape modification on metacommunity structure of tropical bats[END_REF].

Phylogenetic and functional diversity indices (a) Phylogenetic diversity

Phylogenetic diversity (PD) reflects the accumulated evolutionary history of an assemblage, and is based on the evolutionary (patristic) distance between tips in a phylogenetic tree [START_REF] Faith | Conservation evaluation and phylogenetic diversity[END_REF]).

To determine overall PD, we computed the Faith's PD index [START_REF] Faith | Conservation evaluation and phylogenetic diversity[END_REF]) and associated standardized effect sizes (SES) for each of the 7,115 analyzed assemblages. This index assesses species relatedness using the sum of branch lengths of the tree connecting all taxa within an assemblage. Ultrametric phylogenetic distances between the 107 species were retrieved from a Maximum Clade Credibility tree obtained from a sample of 100 dated, calibrated molecular phylogenetic trees assembled by [START_REF] Thuiller | Consequences of climate change on the tree of life in Europe[END_REF]. To assess differences that may exist in phylogenetic structure between habitat types, we also computed the SES values of the mean pairwise phylogenetic distance (MPD) and the mean nearest taxon distance (MNTD), respectively using the independent swap algorithm [START_REF] Gotelli | Null model analysis of species co-occurrence patterns[END_REF]. MNTD and Faith's PD are considered "terminal" indices as they are a more sensitive to patterns occurring at the tips of the tree, whereas MPD is more strongly influenced by the basal structure of the phylogenetic tree (i.e. deeper branches) [START_REF] Mazel | Influence of tree shape and evolutionary time-scale on phylogenetic diversity metrics[END_REF]. Negative SESMPD/MNTD values indicate phylogenetic clustering where coexisting taxa are more related to each other than expected by chance.

Negative values indicate an overdispersed phylogeny where coexisting taxa are less related to each other than would be expected at random. SESPD, SESMPD and SESMNTD values were calculated using the R package picante [START_REF] Kembel | picante: R tools for integrating phylogenies and ecology[END_REF].

(b) Functional diversity

Functional diversity (FD) represents variability in ecological traits among species, reflecting the 'ecological robustness' or resilience of assemblages, which is relevant in order to assess the capacity of an ecosystem to respond effectively to global change [START_REF] Villéger | New multidimensional functional diversity indices for a multifaceted framework in functional ecology[END_REF]. As an overall measure of FD, we used the Rao quadratic entropy index (Botta-Dukát 2005), a generalization of the Simpson's index of biodiversity. Rao's index represents the average dissimilarity between all co-occurring species in the same assemblage and is greater as the number of functionally unique species increases [START_REF] Laliberté | A distance-based framework for measuring functional diversity from multiple traits[END_REF]. We calculated the standardized effect size of Rao (SESRao) using the 'melodic' function [START_REF] De Bello | Functional diversity through the mean trait dissimilarity: Resolving shortcomings with existing paradigms and algorithms[END_REF].

SESRao quantifies the number of standard deviations (SD) that observed Rao values are above or below the mean Rao of random assemblages and it has proven an accurate metric of FD [START_REF] Mason | A guide for using functional diversity indices to reveal changes in assembly processes along ecological gradients[END_REF][START_REF] Mouchet | Functional diversity measures: an overview of their redundancy and their ability to discriminate community assembly rules[END_REF]. Under the null hypothesis of random trait distribution in species assemblages, mean SESRao equals zero. As a result, SESRao negative values indicate trait convergence (i.e. assemblages formed by taxa sharing more similar traits than expected by chance), while positive values indicate trait divergence (i.e. assemblages formed by taxa with more distinct traits than expected at random). We tested whether annual values of FD were on average significantly different from values expected at random by means of Student's t-tests, comparing observed SESRao values with the null hypothesis (μ = 0) [START_REF] Kembel | picante: R tools for integrating phylogenies and ecology[END_REF]). In addition, we obtained a functional dendrogram to compute tree-like FD indices (mean functional distance, MFD, and mean nearest taxon distance, MNTD) and associated SES values (SESMFD and SESMNTD) in an analogous way to that previously reported.

We also computed three complementary multidimensional FD metrics: functional richness (FRic), functional divergence (FDiv) and functional evenness (FEve) [START_REF] Mason | Functional richness, functional evenness and functional divergence: the primary components of functional diversity[END_REF][START_REF] Villéger | New multidimensional functional diversity indices for a multifaceted framework in functional ecology[END_REF]). An interpretation of these metrics is provided in Table 1. As the trait matrix included binary and continuous trait variables, we previously performed the following steps in order to obtain a new set of 'traits' suitable to be used as input when computing the FD indices. First, a generalization of Gower's distance was employed to calculate multivariate distances between species based on the raw trait data [START_REF] Pavoine | On the challenge of treating various types of variables: application for improving the measurement of functional diversity[END_REF]). We made sure the seven trait categories had a similar weight using the 'dis.ktab' function in the R package ade4 [START_REF] Dray | The ade4 package: Implementing the duality diagram for ecologists[END_REF]. The resultant species × species dissimilarity matrix was then subject to a Principal Coordinates Analysis (PCoA) and the six first axes of this ordination were used to build a six-dimensional convex hull from which we computed FRic, FDiv and FEve using the function 'multidimFD' developed by Sebastien Villéger. We computed these three metrics from the complete matrix of 142 traits as well as separately for Grinnellian and Eltonian traits using a matrix of habitat and feeding dissimilarity respectively. We used generalized linear mixed models [START_REF] Pinheiro | Mixed-effects models in S and S-PLUS, 2nd edn[END_REF], in which FRic/FDiv/FEve was considered the response variable and habitat type the explanatory variable to test for differences between forest and farmland assemblages. Models included study year, taxonomic richness, and study plot (to account for variability in the level of indices among sites) as covariates and random effect, respectively.

Relationship between species richness and structure of assemblages

In order to explore the relationships between species richness (SR) and the phylogenetic and functional structure of assemblages, we correlated mean values (averaged over years) of species richness, SESPD and SESRao for each of the 812 study plots by means of Pearson's correlations.

We obtained similar results when considering each community as an independent sampling point (i.e. one point per community; n = 7,115 communities) (analyses not shown).

Temporal trends in phylogenetic and functional diversity

We modeled bird diversity trends over the 2002-2014 timeframe to test whether forest and farmland assemblages exhibit similar temporal patterns. We considered five response variables: SESPD, SESRao, FRic, FDiv and FEve. All these response variables were scaled to null mean and unit standard deviation (z-scores) prior to modeling in order to facilitate direct comparisons among them. We used generalized additive mixed models (GAMM) [START_REF] Wood | Generalized additive models. An introduction with R[END_REF]) in which trends were quantified by a continuous penalized spline with a degree of smoothness estimated by generalized cross validation (constrained to a maximum of 4 to avoid unjustified complexity). We included study plot as random effect to account for variability in the level of indices among sites and controlled residual spatial autocorrelation with an unconstrained twodimensional (latitude, longitude) spline. We further added taxonomic richness as a spline covariate in all models with FRic, FDiv or FEve as response variable to correct statistically the sampling-mediated positive relationship between the number of species and the trait range in an assemblage [START_REF] Barnagaud | Temporal changes in bird functional diversity across the United States[END_REF].

Results

Differences in phylogenetic and functional structure between farmland and forest communities

Although farmland and forest assemblages harbor a quite similar number of breeding common species (mean ± SE; farmland: 34.54 ± 7.00, range: 11-56; forest: 32.97 ± 7.25, range: 10-56), farmland communities show a considerably higher phylogenetic diversity in comparison with forest communities (t = -11.07, p<0.001) (Fig. 1a). Annual average values of phylogenetic diversity (PD) were significantly greater than zero in farmland assemblages throughout the study period, which suggests a prevalence of phylogenetic overdispersion in these communities (Fig. 2). Instead, all annual averages for forest assemblages were negative indicating a trend towards phylogenetic clustering (Fig. 2). Accordingly, standardized effect sizes of mean pairwise phylogenetic distance (MPD) and mean nearest taxon distance (MNTD) values were lower in forest assemblages in comparison with farmland assemblages (SESMPD, forest: -0.49 ± 0.71, farmland: 0.10 ± 0.38; SESMNTD, forest: -0.53 ± 0.75, farmland: -0.22 ± 0.71; Fig. 1a) (both p-values <0.001). In farmland communities, SESPD values increased over the study period, while these showed a 'hump-shaped' pattern in forest assemblages (Fig. 3a).

In terms of functional diversity (FD), although mean SESRao values were negative in almost all cases (Fig. 2), forest assemblages showed significantly lower values in comparison with farmland assemblages (forest: -0.48 ± 0.40, farmland: -0.04 ± 0.33; t = -19.26, p <0.001).

In both habitats, temporal trends for SESRao values followed a similar pattern to that observed for PD (Fig. 3a), indicating a tight correlation between both facets (farmland: ρ = 0.27, p <0.001; forest: ρ = 0.51, p<0.001). The obtained values for both metrics of functional structure (SESMFD and SESMNTD) indicated that forest assemblages exhibit a stronger signature of functional clustering (SESMFD, forest: -0.39 ± 0.35, farmland: -0.11 ± 0.72; SESMNTD, forest: -0.33 ± 0.39, farmland: -0.20 ± 0.70; Fig. 1a) (both p-values <0.001).

Farmland communities also differed significantly from those in forests with respect to two of the three complementary functional diversity indices; functional richness (FRic) and divergence (FDiv) (Table 1). Overall, farmland communities occupy a greater amount of functional space than forest assemblages irrespective of the category to which they belong to (Fig. 4). Temporal trends in FRic and FDiv for farmland and forest assemblages showed divergent patterns (Fig. 3b). Trait range (FRic) and spread (FDiv) within the functional space decreased in the forest assemblages, whereas both indices tended to increase in farmland communities. Bird assemblages appeared to become less functionally even over the study period (Fig. 3b).

When grouping functional traits into different categories (subsets), we observed that, similar to generic FD, the FRic and FDiv of both feeding and habitat traits were greater in farmland assemblages (Table 1). In contrast, the regularity of species' abundance distribution in the functional space (functional evenness) did not differ significantly between farmland and forest assemblages in either case (Table 1). It means that, although the amount of niche space filled tend to vary between farmland and forest assemblages, the way in which this space is filled does not differ between habitats.

Species richness versus structure of assemblages

The sign of the relationship between species richness (SR) and both functional and phylogenetic structure varied between habitat categories. Both PD and FD of forest assemblages increased with increasing SR, indicating that poor assemblages were phylogenetic and functionally clustered and that random patterns have a greater influence in richer assemblages (PD: ρ = 0.37, p <0.001; FD: ρ = 0.38, p<0.001; Fig. 5, Fig. S1). In contrast, in farmland assemblages, the correlation between SR and phylogenetic and functional structure was non-significant in either case (both p-values > 0.1). Rather, the relationship turned out to be asymptotic in both cases (Fig. 5, Fig. S1). In fact, albeit non-significant, the model that included the quadratic term of SR provided a better fit to the data than did the model including the linear term of SR (PD, SR: estimate = 0.37 ± 0.27, t = 1.35, p = 0.17; SR 2 : estimate = -0.79 ± 0.27, t = -2.88, p <0.01; FD, SR: estimate = 0.06 ± 0.26, t = 0.23, p = 0.81; SR 2 : estimate = -0.43 ± 0.26, t = -1.66, p = 0.09).

Discussion

Phylogenetic structure

Although it is expected that severe habitats such as savannas [START_REF] Almeida | The influence of biogeographic history on the functional and phylogenetic diversity of passerine birds in savannas and forests of the Brazilian Amazon[END_REF]) and agricultural landscapes [START_REF] Deikumah | Land use types influenced avian assemblage structure in a forest-agriculture landscape in Ghana[END_REF] tend to harbor assemblages with lower phylogenetic diversity (PD) because the persistence of some lineages is jeopardized, we found the opposite pattern; forest assemblages exhibited lower PD than those located in farmlands. This result indicates that, although farmland intensification has led to a sharp decline in populations of farmland birds, agriculture landscapes still harbor phylogenetically rich communities probably due to the legacy effects of past land-use [START_REF] Scherson | Landscape structure and breeding bird distribution in a sub-mediterranean agroecosystem[END_REF]. The Mediterranean basin has experienced thousands of years of human landscape modification for livestock and agricultural purposes, with forest destruction the most obvious consequence of such an action [START_REF] Blondel | The 'Design' of Mediterranean Landscapes: A Millennial Story of Humans and Ecological Systems during the Historic Period[END_REF]. As a result of traditional landscape design and management by humans, and in discordance with general assumptions on the relationship between biodiversity and habitat degradation, Mediterranean bird communities do not reach their maximum biodiversity in pristine oak woodland, but in agro-sylvo-pastoral landscape mosaics [START_REF] Blondel | Biology and wildlife of the Mediterranean region[END_REF]. In addition, it has been suggested that most of the few bird groups that differentiated within the Mediterranean region are fundamentally birds of open habitats and scrublands (e.g. Sylvia species) [START_REF] Covas | Biogeography and history of the Mediterranean bird fauna[END_REF]. All these factors may have contributed to the existence of a higher diversity in farmland assemblages. In this respect, it should be noted that the Mediterranean region only comprises a relatively small portion of continental France. However, traditional agroforestry systems in other French regions share key features with those typical from Mediterranean environments sensu stricto. For example, original bocages, prevalent in Western France (Atlantic zone), consist of permanent meadows used for livestock farming divided into small hedged fields interspersed with groves of trees, which provide shelter and food for a host of animals. This patchwork-like environment dating back from the Iron Age is evolving rapidly in a trend towards simplification [START_REF] Pointereau | L'arbre des champs : haies, alignements et prés-vergers ou l'art du bocage[END_REF]. Thus, intermediate levels of human-induced disturbance throughout civilizations seem to have contributed to keeping rural landscapes heterogeneous and promoting biological diversity. At the present time, agricultural communities may still bear the imprint of historical landscape configurations, which is thought to persist for decades or centuries due to the lagged response of species to landscape change ("ghosts of landscape past"; With 2007).

On the other hand, most forests in the temperate zone -including Southern Europe-are secondary forests that developed after logging of primeval forests or abandonment of agricultural lands [START_REF] Blondel | Origins and dynamics of forest birds of the Northern Hemisphere[END_REF]. Species inhabiting these forests constitute an impoverished subset of the pool of forest specialists from the Nearctic-Palearctic region, many of them being habitat generalists able to colonize a wide variety of forest environments [START_REF] Covas | Biogeography and history of the Mediterranean bird fauna[END_REF]. Hence, Southern Europe forests would have lost most of the forest specialists a long time ago. In fact, we found that forest communities are made up of species with a lower degree of specialization than those conforming farmland assemblages (Gini specialization index; forest: 0.20 ± 0.03, farmland: 0.22 ± 0.03; see Fig. S2). Such a circumstance, coupled with habitat fragmentation, could explain the existence of strong phylogenetic underdispersion in forest communities. On the contrary, the level of clustering was much lower in farmland assemblages, mostly in plots with significant areas of natural vegetation. It suggests that the presence of structural elements such as field margins, hedgerows and non-cultivated areas (small woodland patches) increases substantially the diversity of lineages that colonizes and persists in these environments (e.g. [START_REF] Doxa | Preventing biotic homogenization of farmland bird communities: the role of High Nature Value farmland[END_REF]. The existence of a subtler signature of environmental filtering may have to do with the fact that domestications of animals and plants began about 4,000 years ago in this region and agriculture environments had been occupied repeatedly by phylogenetically-divergent lineages in multiple colonization events along its long and dynamic history [START_REF] Blondel | The 'Design' of Mediterranean Landscapes: A Millennial Story of Humans and Ecological Systems during the Historic Period[END_REF].

Functional structure

In terms of functional diversity (FD), and in accordance with the phylogenetic results, we found a higher FD in farmland assemblages. A similar result has been previously observed in northern Spain [START_REF] Clavero | Functional homogenization of bird communities along habitat gradients: accounting for niche multidimensionality[END_REF], central Italy [START_REF] Morelli | Associations among taxonomic diversity, functional diversity and evolutionary distinctiveness vary among environments[END_REF]) and tropical montane Ethiopia [START_REF] Gove | Structurally complex farms support high avian functional diversity in tropical montane Ethiopia[END_REF], which indicates that bird communities in certain agricultural landscapes sustain distinct species from an ecological perspective. Whilst, in concert with that reported in terms of phylogenetic structure, forest assemblages showed a higher degree of functional clustering, which means that they harbor a greater number of functionally similar species.

The highest values of functional richness (FRic) were found in communities from lowintensity farmlands (agricultural land with significant proportions of natural vegetation), which represent traditional land uses. At the opposite end, coniferous forests constitute the habitat type supporting the lowest values of FRic probably due to its low spatial heterogeneity in vertical and horizontal structure (Fig. 4). Differences in FRic between farmland and forest assemblages were less pronounced when considering only Eltonian traits because the low richness in terms of species' feeding traits observed in communities located in lands dominated by complex cultivation patterns or non-irrigated land. It suggests a reduced trophic-niche space in the most anthropogenized and simple agricultural landscapes.

Obtained values of functional divergence (FDiv) indicate a higher degree of niche differentiation among species within farmland communities in comparison with forest communities. That is, the most abundant species in farmland assemblages are very dissimilar which may be due to -as above mentioned-limiting similarity precludes the co-existence of functionally redundant taxa, or the availability of a great variety of feeding resources (open fields for seed-eating species, bushes for frugivorous species and groves of trees for those that are canopy-gleaners) and nesting substrates in some agricultural systems. Differences in FDiv between habitats remained when analyzing Eltonian and Grinellian traits separately, which implies the existence of a higher level of differentiation in farmland assemblages at different dimensions of the niche.

Regarding the way in which the biomass of the species assemblage is distributed in niche space (functional evenness, FEve), our results suggests that ecological resources are exploited in a pretty similar manner in both habitats. It means that the availability of niches potentially available but unused by the bird species of the community does not differ between farmlands and forests, a similar result to that reported by [START_REF] Morelli | Associations among taxonomic diversity, functional diversity and evolutionary distinctiveness vary among environments[END_REF]. According to the obtained FEve values, both habitat types exhibit a moderate uniformity of species in functional space (intermediate FEve values; Fig. 4). Thus, by using functional evenness as proxy of community resilience, we hypothesize that farmland and forest assemblages are equally vulnerable to modifications (e.g., climate change or biological invasions).

In terms of functional structure, forest bird assemblages showed a pattern of clustering, contradicting our initial hypothesis of functional repulsion as result of the prevalence of limiting similarity in more benign habitats. This result agrees with that reported in forest bird assemblages from the Cerrado [START_REF] Sobral | Functional and phylogenetic structure of forest and savanna bird assemblages across spatial scales[END_REF] and the Brazilian Atlantic Forest [START_REF] Rocha | The conservation value of cacao agroforestry for bird functional diversity in tropical agricultural landscapes[END_REF]. Although temperate forests are not as structurally complex and dense as rainforests from the tropical region, one possibility is that closed habitats (i.e. forests) act as local environmental filter restricting the occurrence of large birds or those that require open country (e.g. ground-feeding insectivores like pipits and larks).

Relationship between taxonomic richness and structure of assemblages

The shape of the relationships between metrics of biodiversity and taxonomic richness differed between habitats. While in forest habitats, PD and FD increased linearly with increasing species richness with no evidence of a saturating effect, within farmlands such a relationship seems to be asymptotical indicating that beyond a certain threshold (~35 spp.) both dimensions of biodiversity tend to decay (Fig. 5) (see also [START_REF] Hanspach | Functional diversity and trait composition of butterfly and bird communities in farmlands of Central Romania[END_REF]. It suggests that forests are less-redundant systems than farmlands, in which an increase in the number of species not necessarily increase the functional diversity of the community. Such a pattern may be due to the accumulation of generalist species, which can end up displacing specialist species when the environment becomes oversaturated triggering a decline in diversity [START_REF] Clavero | Functional homogenization of bird communities along habitat gradients: accounting for niche multidimensionality[END_REF].

Temporal trends

Phylogenetic and functional diversity of farmland bird assemblages increased over the period indicating that decreasing bird abundances in agricultural landscapes are not associated with a phylogenetic impoverishment or functional homogenization of local bird assemblages. This tendency was not observed in forest assemblages, where diversity trends showed a non-linear pattern. In a similar way, multidimensional FD indices presented distinct temporal trends in both habitat types. Whereas the volume of functional space and the distribution of species within the convex hull have decreased over time in forest assemblages, these metrics tend to show a positive trajectory in farmland communities. It indicates that the drastic declines documented in bird populations across the French countryside have not led to a depaupuration of local communities in terms of PD and FD. This finding could be due to the existence of a time-lag or a decoupling between the mechanisms regulating species populations (e.g. niche processes) and community-level properties that emerge from these species, which heavily depend on biotic interactions. Hence, compensatory dynamics could explain the existence of a decrease in functional richness of forest communities despite most forest species show stable trends if, for example, increases in the populations of forest species are offset by decreases in the populations of transient species [START_REF] Santini | Assessing the suitability of diversity metrics to detect biodiversity change[END_REF]. In this sense, species populations and diversity indices could exhibit different responses to disturbances, which highlights the need for a better understanding of how the impacts of shifts at the species level relate to community-level descriptors (Supp and Ernest 2014).

Lastly, a caveat of this study is that we did not account for imperfect detection of species in our analyses, a factor that may depend on the habitat type in which the survey is conducted. For instance, in a recent study, [START_REF] Kułaga | Bird species detection by an observer and an autonomous sound recorder in two different environments: Forest and farmland[END_REF] compared the number of bird species detected through human observers (in the field) and recordings by autonomous sound recorders (manually analyzed by observers in the lab) within two different habitats -forest and farmland-in Poland. They reported that observers detected more species than recorders in farmland, but not in the forest [START_REF] Kułaga | Bird species detection by an observer and an autonomous sound recorder in two different environments: Forest and farmland[END_REF]. Regarding this, we have to point out that although our results should be interpreted with caution, this potential bias does not affect the main findings of this study, that is, the different trajectory that farmland and forest assemblages exhibit over time.

Conclusions

Several intriguing ideas for both community ecologists and conservationists arise from this study considering that agro-ecosystems constitute one of the most common landscapes throughout the world and represent the main habitat for many bird species [START_REF] Pain | Birds and farming in Europe: The Common Agricultural Policy and its implications for bird conservation[END_REF]. First, our study supports the notion that heterogeneous landscapes in which environmental conditions change significantly within short distances resulting in small-scale mosaics typical of the Mediterranean Basin and Western France, may benefit generalist species and impose limits to the potential ranges of species requiring large continuous areas of forest.

Second, phylogenetic diversity can be a surrogate for high-dimensional trait diversity and thus, by protecting more phylogenetic diversity, we should also protect a greater amount of total functional diversity [START_REF] Faith | Phylogenetic diversity and conservation evaluation: perspectives on multiple values, indices, and scales of application[END_REF]. Third, some of the observed patterns suggest that past land use may generate strong legacy effects on biotic communities and thus, it's necessary to consider the influence of historic landscape configurations when assessing the effect that contemporary configurations have on them. Finally, and regarding this latter, community-level properties may be particularly resilient to land use alterations and thus, their accuracy as indicators of change may be limited. Consequently, gaining a better understanding of how changes at the species population level scale up to shape community properties constitutes a major challenge for biodiversity research.

32

Diversity indices were scaled to mean = 0 and SD = 1 to ease comparisons.

Figure captions 9 10 Figure 1 .

 101 Figure captions 9 10 Figure 1. Differences in phylogenetic diversity (PD) and phylogenetic structure (left), and

Figure 2 .Figure 3 .

 23 Figure 2. Annual mean values (± 95% CI) of phylogenetic (SESPD) and functional diversity

Figure 4 .Figure 5 .

 45 Figure 4. Functional diversity quantified in terms of functional richness, functional divergence,

and functional evenness of farmland and forest bird assemblages in France over the period

2002-2014. Values were standardized as z-cores to facilitate comparison among the different

Acknowledgements

We are indebted to all volunteers who took part in the French breeding bird survey. Julien Renaud helped us to compile information on habitat types and land uses. VGN was funded by a grant from the 'Make Our Planet Great Again' initiative (Campus France) during his stay at the

Data accessibility:

The dataset in this study will be deposited in a public repository upon acceptance.

Conflict of interest:

The authors declare that they have no conflict of interest.

Biosketch

Vicente García-Navas is a FCT researcher at the University of Lisbon. His work addresses how historical, phylogenetic and spatial patterns of eco-morphological and lineage diversity are influenced by biotic and abiotic factors and the interaction between them.

Wilfred Thuiller is a senior scientist at CNRS based in Grenoble. His research is focused in investigating the impacts of global change on biodiversity.

1 the amount of the functional space filled by the community; functional divergence (FDiv) measures the degree to which the abundance of a community is 2 distributed toward the outer margins of occupied trait space and can thought of as functional 'specialization'; functional evenness (FEve) quantifies the 3 regularity of species' abundance distribution in the functional space to allow effective utilization of the entire range of resources available.


E Baudin

C Voisin

B Plaçais

Hyperbolic Phonon Polariton Electroluminescence as an Electronic Cooling Pathway

Engineering of cooling mechanism is of primary importance for the development of nanoelectronics. Whereas radiation cooling is rather inefficient in nowadays electronic devices, the strong anisotropy of 2D materials allows for enhanced efficiency because their hyperbolic electromagnetic dispersion near phonon resonances allows them to sustain much larger (∼ 10 5 ) number of radiating channels. In this review, we address radiation cooling in 2D materials, specifically graphene hexagonal boron nitride (hBN) heterostructures. We present the hyperbolic dispersion of electromagnetic waves due to anisotropy, and describe how the spontaneous fluctuations of current in a 2D electronic channel can radiate thermal energy in its hyperbolic surrounding medium. We show that both the regime of (i) thermal current fluctuations and (ii) out-of-equilibrium current fluctuations can be described within the framework of transmission line theory leading to (i) superPlanckian thermal emission and (ii) electroluminescent cooling. We discuss a recent experimental investigation on graphene-on-hBN transistors using electronic noise thermometry. In a high mobility semimetal like graphene at large bias, a steady-state out-of-equilibrium situation is caused by the constant Zener tunneling of electrons opening a route for electroluminescence of hyperbolic electromagnetic modes. Experiments reveal that, compared to superPlanckian thermal emission, electroluminescent cooling is particularly prominent once the Zener tunneling regime is reached: observed cooling powers are nine orders of magnitude larger than in conventional LEDs.

I. INTRODUCTION

Heat transfer is a key issue in nano-electronics. This is true for the silicon MOSFET technology where cooling limits device performance above modest power density ∼ 100W cm -2 [1] corresponding to a thermal conductance G th < 1 kWm -2 K -1 . With the rise of 2D (bi-dimensional) electronics [2], this problem becomes even more acute because of the minute heat capacity of the channel and the reduced thermal links to the environment. In this respect, the case of graphene-based devices is emblematic : the absence of a bandgap requires driving graphene transistors in the velocity saturation regime where Joule dissipation is prominent. In high-mobility devices, the very limited electron-phonon coupling (nonpolar material, low defect density) associated to a drop of the Wiedemann-Franz heat conduction mechanism (due to a vanishing differential conduction in the saturation regime), yields considerable electronic temperatures, up to several thousands of kelvins [2][3][4]. In this peculiar situation where phonons remain relatively cold and thermal conduction is blocked, heat release through radiative processes becomes a major issue.

It is usually recognized that black-body radiation cooling is rather inefficient in electronic devices. This essentially stems from the large impedance mismatch at optical (infrared) frequencies between vacuum and the electron gas. From a microscopic point of view, this mismatch is related to the small density of electromagnetic modes within the light cone as compared to that of the electron gas at high temperature. This issue has been overcame in recent studies [3,5], by using a so-called hyperbolic substrate. Within the Reststrahlen (RS) bands (which roughly corresponds to the optical phonon modes) hyperbolic materials sustain long-distance propagating hyperbolic phonon-polariton (HPhP) modes that can carry a large momentum [6][7][8][9][10][11][12][51]. When the 2D channel becomes coupled to the HPhP modes (through their evanescent field), the whole density of states can be used to radiate power : this is the so-called super-Planckian regime [7,14]. The enhancement of the radiated power as compared to the regular blackbody can reach up to 5 decades and is not limited by hot phonon effects thanks to the relatively long propagation length of HPhP modes.

Nevertheless, the absolute cooling power remains rather limited (100 kW cm -2 ) in the case of graphene supported on hexagonal Boron Nitride (hBN) due to the narrow acceptance spectral band of HPhPs (∼ 30meV in hBN) compared to the broadband spectrum of the thermal emission of hot electrons in the channel. This is very similar to the well-known poor efficiency of a bulb lamp (defined as the fraction of power emitted in the visible window) because of the broadband emission of thermal fluctuating charges in the infrared. This limitation is overcame in light-emitting diodes (LEDs) that achieve up to 90% electroluminescence conversion efficiency using a non-equilibrium distribution of electrons and holes that recombine across the band gap, thereby funneling efficiently the electrical power into the relevant visible window [15][16][17][18]. Similarly, emission in the narrow spectral RS bands where the substrate sustains hyperbolic modes can be targeted by driving the transistor in the Zener-Klein (ZK) tunneling regime. The Zener current becomes equivalent to an electrical pumping of non-equilibrium electron-hole pairs that efficiently recombine in the hyperbolic modes of the RS bands. This striking effect is readily seen in Figure 1, where the steady-state electronic temperature of a graphene transistor is monitored (by means of noise thermometry [3,[19][20][21][22][23][24]) as a function of injected Joule power [52]. The outstanding observation is a drop of the temperature for increasing injected power above a certain threshold that corresponds to the ignition of HPhP electroluminescence in the hyperbolic hBN substrate.

In this review, we show -based on both experimental results and theoretical proposals -that radiative heat transfer can be considerably enhanced to become a major cooling pathway by either engineering of optical impedance matching between the transistor channel and a hyperbolic substrate or by operating the device in a Zener-like regime where electrical pumping of electron-hole pairs funnels efficiently the injected power into the propagating hyperbolic bands of the substrate. The discussion is based on the emblematic graphene on hBN transistors but can easily be extended to other prospective active material/hyperbolic substrate couples where superPlanckian radiative cooling could yield a major improvement. The paper is organized as follows : in a first section we briefly introduce hyperbolic materials and the concepts of thermal coupling between a 2D electronic channel and a bulk hyperbolic substrate; then we investigate the case of non-equilibrium situations and the electroluminescence of HPhPs. Finally, we present experimental data for a bilayer graphene on hBN transistor where the electronic temperature is monitored through noise thermometry.

II. COUPLING BETWEEN GRAPHENE AND HYPERBOLIC MATERIALS

Computing the thermal radiation of an electronic channel in the electromagnetic modes requires the knowledge of the electromagnetic response of the dielectric encapsulating medium which is provided by the dielectric permittivity. For an isotropic polar dielectric, the optical phonons couple to light and provide an effective dielectric response of the form

(ω) = ∞ ω 2 LO -ω 2 -jωγ ω 2 T O -ω 2 -jωγ , (1)
where ω LO is the longitudinal optical frequency, ω T O is the transverse optical frequency and γ is the decay rate of the optical phonon resonance. For high-quality materials, Re( (ω)) is negative in the range [ ω T O , ω LO ], as a consequence, the material is opaque in this frequency range called Reststrahlen band (RS) usually in the mid-infrared ([170, 200]meV for hBN). A RS band is characterized by its central frequency Ω = (ω LO + ω T O ) /2 and its width ∆ω = ω LO -ω T O . An anisotropic dielectric material posesses a tensorial dielectric response. For uniaxial crystals, a class including most dielectric 2D materials, in-plane (transverse t) and out-of-plane (z-axis) optical phonon energies are sufficiently separated to give rise to a splitting of the Reststrahlen band into two bands. As an example, hBN possesses two well separated bands corresponding to out-of-plane OPs (type-I band, [97,103]meV), and in-plane OPs (type-II band, [170,200]meV). The energy and width of hBN RS band II makes it a much more efficient heat conduction channel, therefore, by sake of simplicity, we shall focus specifically on this band in this review (with Ω ≡ Ω II ). [START_REF]A discussion of hyperbolicity type can[END_REF] As a consequence, for layered 2D heterostructures, the dispersion of electromagnetic extraordinary waves takes the form :

k 2 z t + k 2 t z = ω 2 c 2 (2)
and, within a type II RS band, t < 0 < z which means that the dispersion relationship is hyperbolic. Note that this occurs only for materials with low-loss optical phonon modes ∆ω/γ 1; a category named hyperbolic phonon materials. The important consequence is that k t > ω/c becomes possible within the dispersion relationship, which opens up propagative modes with high momenta that can help radiate energy.

The efficiency of hyperbolic radiative cooling depends critically on two parameters: (i) the ability of the dielectric material to efficiently couple to the current fluctuations of the metal, and (ii) the propagation depth of this radiation within the hyperbolic dielectric. These two aspects are best understood within the framework of transmission line theory [54]: Considering a 2D conductor with random [55] current density fluctuations j(ω, q), where q is the inplane wavector, this source interacts with the TM electromagnetic modes of the substrate of optical impedance Z(ω, q) and the metal conductivity σ(ω, q) in parallel. As a consequence, the power cast by this source in its electromagnetic environment is dP = 1 2 Re(Z -1 ) |Z -1 +σ| 2 |j| 2 . Note that the optical impedance of the dielectric environment Z(ω, q) = Z xy = |E x |/|H y | is defined as the ratio between the longitudinal electric field collinear to the currents in the material, and the transverse magnetic field of the extraordinary wave. It fully determines the electromagnetic response of the environment as seen from the channel's perspective. In terms of optics, the power dP can be interpreted as the average flux of the Poynting vector through the channel/dielectric interface.

According to the fluctuation-dissipation theorem, for a single band, the surface density of the current fluctuations in a material of conductivity σ at temperature T reads

|j| 2 = 1 (2π) 3 4Re (σ) ω n ph (T, ω) d 2 qdω, (3)
where

n ph (T, ω) = [exp( ω/k B T ) -1] -1
is the Bose occupation number. In our experiment, the hBN layer remains cold due to the presence of the gold backgate electrode. Therefore, we can neglect the back thermal flux from the hBN to the channel. The total radiative power is then [56]

P rad = 1 (2π) 2 ˆ∞ 0 dω ˆ∞ 0 qdq M(ω, q) ω n ph (T, ω), (4)
where 0 ≤ M ≤ 1 is the impedance matching factor :

M(ω, q) = 4Re Z -1 Re(σ) |Z -1 + σ| 2 . ( 5
)
The matching factor (5) restricts the phase space available for radiative heat exchange. Indeed the optical impedance reads Z = Z0 t kz k0 , where k 0 = ω/c, and Z 0 = 377 Ohms is the vacuum impedance. When the photon energy is outside the RS band, k z becomes imaginary when q > k 0 √ ε t , and thus M = 0. As a consequence, radiation is limited to the very narrow lightcone as represented on Figure 2 and the emitted power is very small. For monolayer graphene at 2000 K in the vacuum, it represents a radiated power Mσ SB T 4 1.5 Wcm -2 , where σ SB = 5.67 10 -8 Wm -2 K -4 is the Stefan-Boltzmann constant, and M ∼ 2.3% is the average far-field emissivity of graphene [25]. This corresponds to a characteristic thermal conductance P T 10 Wm -2 K -1 . In contrast, within the RS band, k z remains real, and hence M remains sizeable for much larger q vectors, enhancing considerably the emission efficiency.

At this point we distinguish two possible causes for current fluctuations: (a) thermal fluctuations of charges (blackbody radiation) and (b) non-equilibrium fluctuations of charges (electroluminescence). The next two sections will address each of them.

A. Thermal hyperbolic cooling

Matching Factor

Let us first consider thermal emission by a single band. Since RS bands are narrow ∆ω Ω, the photon population is nearly constant over a RS band and given by the occupancy at mid-band n ph (T, ω) n ph (T, Ω). The radiated power (4) in the RS band energy then reads

P rad = A F k 2 F Ω n ph (T, Ω)∆ω , (6)
with

A F = 1 (2π) 2 ∆ω ˆωLO ω T O dω ˆ∞ 0 xdxM(ω, k F x), (7)
where k F is the Fermi wavevector. A F can be interpreted as an average impedance matching factor over the RS band. For a metal with parabolic dispersion, A F π Ωn ph (T, Ω)∆ω also corresponds to the average power radiated per electron in the conduction band. Let us note that in the degenerate case, integral over wavevector ( 6) is bounded in a window due to a cut-off defined by Pauli blocking. Since most optical modes lie out of the lightcone k F √ ε t Ω/c, the optical impedance is well approximated by

Z Z0 √ ε ∞ ⊥ ε ∞ z ω 2 LO -ω 2 ω 2 -ω 2 TO q k0
in the infinite dielectric-thickness limit [57]. Regarding the 2D electronic channel response, due to the significant extension of the integration window in the first Brillouin zone, the nonlocal conductivity is needed. We rely on the linear response theory for the case of non-interacting electrons in a single parabolic band:

σ(ω, q) = iG 0 g s g v /2 (2π) 2 ω q 2 ˆdk f k -f k+q ω + E k -E k+q + i η | χ k |χ k+q | 2 , (8)
where G 0 = 2e 2 /h = 77 µS is the quantum of conductance, g s = 2 and g v = 2 are the spin and valley degeneracy,

f k = 1 exp(E k -µ)/k B T +1
is the Fermi-Dirac occupation factor for state of wavector k, E k is the band energy dispersion, µ is the chemical potential, η is the effective relaxation rate of electrons, and |χ k is the periodic part of the Bloch wavefunction.

By sake of simplicity, we make the following assumptions:

1. Isotropic parabolic dispersion: Except for monolayer graphene, isotropic parabolic dispersion covers most naturally occurring cases. The specific case of graphene has been dealt in detail in the literature [14].

2. Unitary Bloch wave projection: Due to the small extent of the Fermi sea compared to the Brillouin zone, we can approximate

| χ k |χ k+q | 2 ∼ 1.
The computation of the nonlocal conductivity is given in appendix. In the degenerate case, one has σ(ω, q) = G 0 gsgv 4π S(µ/ ω, q/k F ), where S is a simple geometric function. Consequently, the prefactor A F is a function of the reduced doping µ r = µ/ Ω II and of an effective light-matter coupling factor that reads

α HP hP = Re(Z(Ω, q c ))Re(σ(Ω, q c )) = α F S g s g v 2π ε ∞ ⊥ ε ∞ z k F k 0 , (9)
where α F S = G 0 Z 0 /4 = 1/137 is the fine-structure constant, and

q c = k F / √ µ r = 2m * Ω/
, where m * is the effective mass of the parabolic electronic band. The doping dependence of α HP hP reflects the ability of the channel's electrons to couple to the HPhP modes of the substrate. Figure 3 represents the matching factor taking into account only the conduction band of bilayer graphene for increasing doping. The allowed transitions in the degenerate regime restrict the HPhP radiation window in the wavevector range

[( 1 + µ -1 r -1)k F , 1 + µ -1 r + 1)k F ], which is centered on q c = k F / √ µ r = √ 2m * Ω II / 4.43 10 8 m -1 .
For semimetals (like few layer graphene), the wavevector window for impedance matching is larger and extends over [0, √ 2q c ] in the degenerate case, due to the allowed interband transitions in the thermal regime. For moderately doped graphene, µ < Ω II /2, the weak electron-photon coupling α F S is balanced by the characteristic electronic transition wavevector with respect to optical wavevector q c /k 0 ∼ 460, so that α HP hP ∼ 1. We conclude that graphene is naturally matching the naturally hyperbolic material hBN. The number of layers of graphene itself impacts obviously the radiated power due to different electronic dispersion. As a rule of thumb, we estimate a radiative cooling power of 1.6 GWm -2 at a temperature of 2000 K which translates into a thermal conductance 0.8 MWm -2 K -1 intermediate between the experimental values in Figure 1. Finally, note that experimental observations reveal that HPhP cooling is suppressed when graphene is submitted to a quantizing magnetic field [24] because electrons localize in small cyclotron orbits forbidding wavevector matching with HPhPs.

Propagation depth of the HPhPs

The above evaluation of the radiated power doesn't provide information on where heat is dissipated. To obtain an efficient heat release, the propagation depth of HPhPs must be large enough to release heat remotely in a large volume and avoid hot-phonon effects. The characteristic propagation depth of the HPhPs is given by l = (Im(K z )) -1 . For large in-plane k-vector compared to k 0 , the characteristic wavevector in a RS band is K z -ε ⊥ εz q. At center frequency Ω, field propagation is maximal and reaches

l ∆ω γ ε ∞ z ε ∞ ⊥ 1 q . ( 10
)
In the case of bilayer graphene on hBN, the characteristic wavector of near-field radiative heat exchange is q = k F ∼ 200k 0 , for which the propagation depth reaches 240 nm. Note that the presence of a nearby metallic back gate, with a low optical impedance, will reflect HPhPs and eventually lead to the reabsorption of the emitted power after one round trip. This situation occurs in thin hBN samples such as in Reference [3]. This is avoided in the back-gated device analyzed below where the hBN thickness (200 nm) exceeds δ/2 120 nm, so that a full HPhP cooling can be observed.

B. Electroluminescence in a Hyperbolic Material

Electroluminscence is associated with interband optical conductivity occurring when electron-hole pairs are electrically pumped from the valence to the conduction band, a situation which occurs spontaneously in high mobility graphene when a large bias is applied [3]. The electroluminescence results from the interband recombination of electron-hole pairs. Whereas classical electroluminescence is limited to vertical transitions, this emission process can also be considerably boosted when high-k propagating modes are available in the dielectric. In this section we focus on the excess radiation power related to electroluminescence.

The theoretical treatment of electroluminescence is closely related to near-field thermal emission, except that out-ofequilibrium current fluctuations are given by the nonlocal van Roosbroeck-Shockley relation for interband transitions [26][START_REF] Rosencher | Optoelectronics[END_REF][START_REF] Greffet | [END_REF]:

|j inter | 2 = 1 (2π) 3 4Re (σ inter ) ω n ph (T, ω) d 2 qdω, (11)
where

n ph (T, ω) = [exp( ω -µ ph )/k B T -1]
-1 with µ ph = µ λ -µ λ is the photon chemical potential associated to the difference in chemical potential between the two bands of interest (denoted as λ and λ ). Note that this approximation assumes an ultra-fast electron-electron thermalization, and the existence of two pseudo Fermi distributions in the valence and conduction bands. The interband nonlocal conductivity, σ inter (ω, q), is evaluated on the interband transitions between the bands of interest. The total radiated power is the sum of power radiated by all intraband and interband transitions. Focusing on the case of a pair of valence and conduction bands, the fraction of radiated power cast by electroluminescence is deduced from Eqs.( 4) and (5) using j inter as the source. Note that, in this case, the impedance matching coefficients read

M(ω, q) = 4Re Z -1 Re (σ inter ) |Z -1 + σ inter + σ intra | 2 , ( 12
)
where σ intra is the intraband conductivity. With gaped semiconductors, electroluminescence becomes important when the difference in conduction and valence chemical potential reaches the bandgap energy µ c -µ v E g . In contrast to this well-known configuration, graphene has no bandgap, however, the threshold energy is now the RS band energy µ c -µ v Ω II . For unbiased graphene, out-of-equilibrium electronic distribution can be transiently created by optical pumping and have been observed to relax on sub-ps timescales by optical pump-probe techniques [2,5,29,30]. In the next section, we show that electrical pumping of charges occurs similarly in high mobility graphene under large bias and leads to a favorable situation for HPhP electroluminescence cooling. Noise temperature TN = SI / (4G ds kB) deduced from the microwave excess current noise SI and the differential conductance G ds . Noise reveals a striking difference between the intraband and interband regimes with a superlinear bias dependence characteristic of "hot" electrons (dotted black line) in the former, dropping toward a linear dependence characteristic of "cold" electron regime in the latter. The theoretical value for the threshold between electroluminescent cooling and heating is sketched as the dashed line. As sketched in the inset, the tunneling process creates an out-of-equilibrium electron-hole distribution which relaxes under HPhP emission.

To illustrate HPhP cooling we analyze here one of the graphene ZKT transistors investigated in Reference [3] (Figure S4 of the Supplementary Material). It is a high-mobility bilayer graphene flake, of dimensions L × W = 3.6 × 3 µm, deposited on a 200 nm-thick hBN crystal (inset of Figure 4-a). A global metallic back gate is used to tune electron density in the range -n = [0, 2.1 × 10 12 ] cm -2 . Figure 4-a depicts the transport properties of the ZKT transistor which is characterized by an initial increase of the intraband current up to a doping-dependent saturation value I sat followed by an interband contribution with a doping-independent differential conductance G zk 0.45 mS. Figure 4-b shows the noise temperature characteristics.

A. Intraband current saturation

At moderate bias in high mobility devices, we observe a current partial saturation [31] which is well described by the semi-empirical expression of the conductivity

g ds = G ds × L/W = |n|eµ/[1 + (E/E sat ) 2 ] [3]
. Fits of the differential conductivity allows obtaining mobility µ 1 m 2 V -1 s -1 and a doping-dependent saturation electric field E sat 0.3 V µm -1 . Mobility is large enough for intraband current saturation to be achieved at low to moderate voltage V ds 1 V, leaving a broad experimental window for inter-band transport. The current saturation characteristic energy E sat is deduced from the saturation electric field E sat with E sat = hneµ e(4k F /π) E sat . [3] At the low to moderate doping of our experiment, we remark that current saturation is limited by the Fermi energy E sat = E F [3].

B. Zener-Klein tunneling

We recall here the properties of the high-bias Zener-Klein regime introduced in Reference [3]. Assuming a uniform electric field, one can express the doping-independent Zener-Klein resistance as the series addition of L/l zk incoherent tunneling processes taking place over a length l zk with an elementary conductance G 1 = 1 2 β zk M G 0 , where M = k F W/π is the number of transverse electronic channels involved in electronic conduction, and β zk 1 is a smooth junction factor [58] (bilayer graphene [3]). This leads to an apparent Zener-Klein conductivity σ zk = β zk G 0 × M l zk W ∝ β zk k F l zk , because each electron-hole pair spontaneously created by Zener tunneling diffuses over a length l zk before recombining separately with conjugated charge carriers. Experimentally, one can extract the product β zk l zk ∝ 1/k F from the constant differential conductance G zk 0.45 mS. Zener-Klein tunneling being controlled by Pauli blocking, its ignition is subject to a threshold field E zk 2|E F |/el zk . If mobility is large enough, the intraband current is saturated prior to the onset of Zener-Klein tunneling so that J(E zk )

J sat . This defines the intraband/interband transport boundary as the line I sat = 2 β zk G zk V zk , highlighted as a black dotted line in Figure 4-a. Note in passing that saturation currents do reach the intrinsic limit (symbols in the dotted line of Figure 4-a) in this sample. From the slope of the line we extract β zk 0.5, and finally deduce l zk (|n|) (lower inset of Figure 4-a). The l zk ∝ 1/ |n| dependence, indicates that the electron-hole diffusion length is proportional to the average distance between charges. Remarkably enough, transport does not show any indication of thermal degradation at high bias in spite of the large Joule power P = 200 kWcm -2 , which is three orders of magnitude larger than the limit for thermal degradation in MOSFETS [1]. To get a deeper insight on the physical origin of this hyper-cooling we rely on the noise analysis below.

C. Hyperbolic phonon polariton cooling

As explained in the experimental section below, a Johnson-Nyquist noise temperature T N = S I (ω)/4k B G ds can be deduced from the GHz current noise S I and differential conductance G ds . Its bias voltage dependence T N (V ds ) in Figure 4-b reveals very contrasted behaviors for intraband and interband transport. The former corresponds to a a superlinear bias dependence T N ∝ V 1.5 ds at low bias [3] and a quasi-linear dependence T N ∝ V ds at high bias. Remarkably the transition occurs at the Zener-Klein tunneling onset defined above according to Pauli blocking (black dotted line), and gives rise to a dramatic drop of the noise temperature (or a noise plateau for graphene transistors with thin hBN dielectric [3]). The neutrality regime provides a valuable insight on the cooling mechanism, by featuring a threshold at V ds = 0.2 V, which is invariantly observed in graphene allotropes (single-layer, bilayer and trilayer graphene), pointing to a phonon activation mechanism with an energy 0.2 eV matching the type-II HPhP band of hBN [3]. This activation energy corresponds to the "optical gap" Ω II expected for HPhP electroluminescence.

If the low-bias regime can be interpretaed in the framework of superplanck HPhP black-body emission, the high-bias regime cannot (for details see Reference [3]). In particular, superplankian emission cannot account for the drop of temperature in increasing Joule power. An out-of-equilibrium emission of HPhPs by hot electrons in the inter-band regime is needed.

D. Electroluminescence hyperbolic phonon polariton cooling

A parallel with LED is helpful to understand why electroluminescence is associated to a cooling of the electron gas. In LEDs, electron-hole pairs are injected from the electrodes and recombine in the active region. The energy per injected electron-hole pair E inj is slightly larger than the effective bandgap at the junction, and consequently, radiative recombination energy E g < E inj is insufficient to sink the injected energy (the situation is made even worse considering the competing nonradiative recombination). As a consequence, electroluminescence is usually associated with a small heating of the device. However, it has been reported that, for high quality devices at subthreshold bias for which E g > E inj , the LED acts as a refrigerator as the missing energy E g -E inj is provided by the thermal energy of the electron gas. To date, LED refrigerators are very inefficient with a record cooling power of only 8 pW [16]. The graphene transistor presented in this paper operates in refrigerator mode at moderate bias where E inj < Ω II , therefore, the sudden drop in temperature at the Zener-Klein tunneling threshold is a consequence of a competition between intraband heating, and interband HPhP electroluminescent cooling. The cooling power reaches ∼ 10mW, i.e., nine orders of magnitude larger than in conventional LEDs. However, since the electron gas temperature is much larger than the phonon bath one, as a whole this graphene transistor is not a refrigerator and cooling power is not bounded by a Carnot inequality, in contrast with LEDs where the device gets colder than its environment.

The injected energy per electron-hole pair is the electrical work done by the charge during diffusion and reads E inj = P zk / ṅzk , where P zk = σ zk (E -E zk )E is the Joule power associated with Zener tunneling processes, and ṅzk = G 0 M (E -E zk )/e is the rate of electron-hole pair creation by Zener-Klein tunneling. We deduce the injection energy per electron-hole pair to be E inj = β zk l zk eE. This allows obtaining the characteristic electric field for which injection energy balances HPhP emission energy Ω II by electron-hole pair recombination. This electric field E th = Ω II β zk l zk e defines the boundary between electroluminescent cooling and heating, it is indicated in Figure 4-b (dashed line) and it matches well the voltage for which the minimal temperature is reached.

The chemical potential imbalance between the valence and conduction band can be estimated from the data in Figure 1. For -n = 1.25 × 10 12 cm -2 doping, the maximum temperature T N 3200 K is reached at ZK threshold for an injected Joule power 0.25 mW µm -2 , and drop to a minimal temperature T N 1400 K for 1 mW µm -2 . Assuming the cooling power variation with the electronic bath state is mainly governed by the photon occupation factor P rad → n ph (a legit assumption because of its singular behavior at the electroluminescence-superluminescence threshold Ω II = µ c -µ v ), we deduce that the chemical potential imbalance µ c -µ v 0.16 eV ω T O reaches the lower boundary of the hBN RS band II. Below Zener-Klein electric field threshold, the electronic temperature results from a power balance between intraband processes only. At very high electric field, with Zener-tunneling and electroluminescence recombination increasing, the conduction and valence band reach a doping state which is mainly independent of gate-doping:

µ c -µ v ∼ Ω II (µ c + µ v )/2 = µ 0
, where µ 0 is the Fermi energy at vanishing bias. As a consequence, the power balance tends to be dominated by Joule heating originating from ZK tunneling processes, which is independent of doping, and HPhP emission by both intraband thermal emission and interband electroluminescence defined by a similar electronic distribution. Therefore, noise temperatures obtained at finite doping tend to line up on the temperature obtained at null doping at large bias as observed on Figure 1.

IV. CONCLUSIONS AND PERSPECTIVES

The above noise thermometry measurements in the Zener-Klein regime of a hBN-supported bilayer graphene transistor, supported by the theoretical analysis, demonstrate the superior cooling power of radiative HPhP electroluminescence. As shown in Reference [3] the mechanism is quite generic and consistently observed in single layer, bilayer and trilayer graphene. The role of hBN thickness is emphasized here with a stronger cooling effect, and a negative thermal conductance regime reported in the semi-infinite hBN substrate regime. Theoretical modeling provides useful keys to assess the difference between supported and encapsulated graphene, and constitutes a base to tackle more involved situations such as the rich moiré phases of hBN encapsulated graphene [32][33][34]. Hyperbolicity of the phonon modes adds to the list of hBN qualities which include dielectric strength, compatibility with high graphene mobility [35,36], and optical properties [37,38].

HPhP radiative cooling is not limited to the hBN-graphene couple. It should also be observed with graphene on calcite [39] or mica [40] which have similar hyperbolic Restrahlen bands. Graphene on mica transistors have a lower mobility possibly precluding the electroluminescent regime but still constitute are valuable route for high performance flexible electronics [41]. Conductive materials such as cuprates or tetradymites (e.g. Bi 2 Se 3-x Te x ) also sustain hyperbolic modes [42]. The case of tetradymites deserves special attention as HPhP energy lies in the near infrared to optical range [43], which is larger than the bandgap. Their finite conductivity comes from topological surface states which can be used to realize field-effect transistors in the bulk depleted regime [44]. Remarkably they combine, in one and the same material, the two properties valued in graphene-hBN heterostructures: a gapless Dirac Fermion transport at the surface and hyperbolic HPhP modes in the bulk. Finally, and possibly most impactful for applications, superplanckian HPhP cooling can be envisioned in silicon on insulator (SOI) technology [START_REF] Taur | Proc. IEEE[END_REF], provided that the 10 nm-thin active channel is nested in the near-field of a hyperbolic dielectric such as α-sapphire or αquartz [START_REF] Winta | [END_REF]. The inclusion of an HPhP cooling stage in MOSFET technology would lift the thermal bottleneck at ∼ 100 W cm -2 reported in Reference [1].

V. EXPERIMENTAL SECTION A. Sample

Sample fabrication is described in Reference [3]. The sample analyzed here is made of an as-exfoliated bilayer graphene flake, of dimensions L × W = 3.6 × 3 µm, deposited on a 200 nm-thick hBN crystal (inset of Figure 4-a).

The use of a thick hBN dielectric (i.e. thicker than the HPhP propagation depth) avoids HPhP reflections at the metallic back gate and simplifies analysis. The sample is equipped with high-transparency contacts, a global metallic backgate, and is embedded in a coplanar wave guide (CPW) for high-frequency current noise measurement. DC current and noise are measured at liquid helium temperature (4.2K), in the hole doped regime -n = 0-2.1 × 10 12 cm -2 , where we obtain a low-field mobility µ 1 m 2 V -1 s -1 .

B. Noise thermometry principles

Noise thermometry principles, described in Reference [3] (Supplementary), are introduced here. Electrons in a metal sustain thermal fluctuations leading to the well-known Johnson-Nyquist white current noise S I = 4Gk B T e [47], proportional to the sample conductance G and electronic temperature T e . This Nyquist formula serves as a primary thermometry and is used to estimate the out-of-equilibrium electronic temperature at finite bias [3,4,[19][20][21][22]24], when the heated electron gas decouples from the lattice and gets thermalized thanks to efficient electron-electron relaxation (rate ∼ 10 fs). To account for self-heating in the Nyquist formula, one substitutes T e by a noise temperature T N T e L,W [59], and the linear conductance by the differential conductance G ds = ∂I ds /∂V ds , taken at the noise measuring frequency.

Experimentally the finite currents entail an additional 1/f -like conductance noise, S I = α H I 2 /N f with α H ∼ 2 10 -3 is the Hooge constant and N = |n|LW the total number of carriers participating to transport [48], obscuring thermal noise of interest at low frequency. Writing thermal noise as S I = 2eFI, where F ∼ 0.1 is a typical thermal Fano factor, one estimates a 1/f -corner frequency, f c ∼ α H 2F × v d /L. The large drift velocities v d 3 × 10 5 m s -1 of our high-mobility graphene yields f c 1 GHz, whence the need to measure thermal noise in the microwave band f = 1-10 GHz.

VI. SUPPORTING INFORMATION A. Impedance matching in the degenerate case

In this appendix, we compute a generic expression for the nonlocal conductivity of a single parabolic electronic band in the degenerate regime k B T µ. We then compute the coupling factor A F in the case of HPhP thermal radiation within a RS band of energy Ω and width ∆ω.

In the degenerate case, σ(ω, q) = iG 0 g s g v /2 (2π)

2 ω q 2 ˆdk n k -n k+q ω + E k -E k+q + i0 + , (14)
σ(ω, q) = iG 0 ω 1 q 2 g s g v /2 (2π) 2 (I 1 (q, ω) + I * 1 (q, -ω)) .

Focusing now on I 1 (q, ω), for a single isotropic massive band

E k = 2 |k| 2
2m , and at null temperature and Fermi energy µ, n k = Θ(µ -E k ), and we have

I 1 (q, ω) = ˆ2π 0 dθ ˆ∞ 0 kdk Θ(µ -E k ) ω + E k -E k+q + i0 + . ( 16
)
Defining x = k/k F , and y = q/k F , we have

I 1 (q, ω) = k 2 F µ ˆ1 0 xdx ˆ2π 0 dθ (ω+i0 + ) µ -y 2 -2xy cos θ , ( 17
)
The integral is easily done and yield

I 1 (q, ω) = 2π k 2 F µ 1 2y    ω µ -y 2 2y - ω µ -y 2 2y 2 -1    (18)
for ω µ > (q/k F ) 2 , whereas the case ω µ < (q/k F ) 2 is obtained by symmetry of the integral

I 1 (q, ω) = 2π k 2 F µ 1 2y    ω µ -y 2 2y + ω µ -y 2 2y 2 -1    . ( 19
)
Figure 5 represents the computed A F factor at null temperature function of the channel/insulator properties. Corresponding illustration of the nonlocal impedance matching factor is represented on Figure 3.

The nonlocal conductivity in the degenerate case has been computed in detail for monolayer [49] and bilayer [50] graphene and allows computing the radiated power using the above formalism (see [14] for the case of monolayer graphene).

VII. ACKNOWLEDGMENTS

The research leading to these results have received partial funding from the European Union "Horizon 2020" research and innovation program under grant agreement No. 785219 "Graphene Core", and from the ANR-14-CE08-018-05 "GoBN".

FIG. 1 :

 1 FIG.1: Hyperbolic phonon polariton cooling of a bilayer graphene transistor. The noise temperature TN shows a nonmonotonic dependence in the applied Joule power with a characteristic drop at a doping-dependent threshold between a lowbias hot electron regime and a high-bias cold-electron regime signaling the onset of a new and very efficient cooling mechanism: the electroluminescence of mid-infrared hyperbolic phonon polaritons. The cooling is characterized by the heat conductance G th = P/Te P/TN which switches from 50 kWm -2 K -1 at low bias to 1.3 MWm -2 K -1 at high bias. Colors represent doping with respective carrier densities -n = 0 (red), 0.18 (tangerine), 0.35 (orange), 0.53 (amber), 0.70 (yellow), 0.88 (chartreuse), 1.05 (green), 1.23 (turquoise), 1.40 (cyan), 1.58 (cerulean), 1.75 (blue), 1.93 × 10 12 cm -2 (dark blue). Electroluminescent cooling starts when Zener tunneling begins. Zener tunneling threshold is represented as a dashed line and increases with doping because of Pauli blocking. Dotted line represents the Joule power for which electroluminescent cooling balances Joule heating due to Zener current, see section III.

FIG. 2 :

 2 FIG.2: Real part of the optical admittance (Z -1 ) of hBN (omitting RS band I) on two in-plane wavevector scales to illustrate the far field and HPhP extension of radiation phase space. As discussed in section II A, the impedance matching window in the RS band is upper bounded by at least 2kF .

FIG. 3 :

 3 FIG. 3: Impedance matching factors M (color scale) in RS band II of hBN (single band approximation for bilayer graphene with m * = 0.03m0) for increasing doping (from bottom to top, µr = µ/ ΩII = {0.1, 0.5, 1.0, 5.0}). Wavevector is given in units of qc = 2m * ΩII / 4.4 10 8 m -1 . The emission window increases with doping up to ∼ 2kF at large doping, while the average matching is reduced over the RS band.

2 β

 2 FIG. 4: Transport and noise of a hBN-supported graphene ZKT transistor. Panel a) current saturation and Zener-Klein regimes in a high-mobility bilayer graphene flake of dimensions L×W = 3.6×3 µm. Sample (upper inset) is deposited on a 200 nm thick hBN flake acting both as a bottom gate dielectric (capacitance Cg = 0.14 mF m -2 ) and hyperbolic phonon polariton radiator. Carrier density is tuned in the hole doping range -n = 0 (red), 0.35 (orange), 0.70 (yellow), 1.05 (green), 1.40 (cyan), 1.75 (blue), 2.10 × 10 12 cm -2 (purple). At low bias intraband current dominates up to a threshold field Esat = vsat/µ set by the saturation velocity vsat 3 × 10 5 m s -1 and mobility µ 1 m 2 V -1 s -1 . At high bias, above a Pauli unblocking voltage V zk (see text), intraband current is fully saturated and an additional current arising from interband Zener-Klein tunneling with a constant conductance G zk = 0.45±0.05 mS sets in. The boundary between intraband and interband transport regimes follows theoretical expectation Isat 2 β zk G zk V zk (black dotted line with the squares corresponding to theoretical values of current). Panel b)Noise temperature TN = SI / (4G ds kB) deduced from the microwave excess current noise SI and the differential conductance G ds . Noise reveals a striking difference between the intraband and interband regimes with a superlinear bias dependence characteristic of "hot" electrons (dotted black line) in the former, dropping toward a linear dependence characteristic of "cold" electron regime in the latter. The theoretical value for the threshold between electroluminescent cooling and heating is sketched as the dashed line. As sketched in the inset, the tunneling process creates an out-of-equilibrium electron-hole distribution which relaxes under HPhP emission.

2 ω

 2 mid -∆ω/2 dω ˆ∞ 0 xdxM(ω, k F x),

FIG. 5 :

 5 FIG. 5: SuperPlanckian material dependent factor at null temperature AF function of the HPhP fine structure constant α HP hP and the reduced electronic chemical potential µr. The prefactor is AF = 0.023 at α HP hP = 1, µr = 1.