domain: sde sdv

text:
François Duchenne
email: francois.duchenne@wsl.ch

Emmanuelle Porcher

Jean- Baptiste Mihoub

Grégoire Loïs

Colin Fontaine

Tim Coulson

Gabor L Lovei

Controversy over the decline of arthropods: a matter of temporal baseline?

Keywords: insects, population trends, shifting baseline syndrome, temporal trends, time series, abundance

Recently, a number of studies have reported somewhat contradictory patterns of temporal trends in arthropod abundance, from decline to increase. Arthropods often exhibit non-monotonous variation in abundance over time, making it important to account for temporal coverage in interpretation of abundance trends, which is often overlooked in statistical analysis. Combining four recently analysed datasets that led to contrasting outcomes, we first show that temporal abundance variations of arthropods are non-monotonous. Using simulations, we show non-monotony is likely to bias estimated linear abundance trends. Finally, analysing empirical data, we show that heterogeneity in estimated abundance trends is significantly related to the variation in temporal baseline of analysed time series. Once differences in baseline years, habitats and continents are accounted for, we do not find any statistical difference in estimated linear abundance trends among the four datasets. We also show that short time series produce more stochastic abundance trends than long series, making the dearth of old and long-term time series a strong limitation in the assessment of temporal trends in arthropod abundance. The lack of time series with a baseline year before global change acceleration is likely to lead to an underestimation of global change effects on biodiversity.

Introduction

Over the last decades, many studies have documented a biodiversity crisis on the basis of high extinction rates [START_REF] Dirzo | Global State of Biodiversity and Loss[END_REF][START_REF] Ceballos | Accelerated modern human-induced species losses: Entering the sixth mass extinction[END_REF] and population losses or declines [START_REF] Butchart | Global Biodiversity: Indicators of Recent Declines[END_REF][START_REF] Ceballos | Biological annihilation via the ongoing sixth mass extinction signaled by vertebrate population losses and declines[END_REF], mainly in vertebrates. Worldwide erosion of biodiversity is caused by anthropogenic global change [START_REF] Sala | Global Biodiversity Scenarios for the Year 2100[END_REF], i.e. a set of pressures including among others land use change, climate warming, overexploitation, or pollution [START_REF] Steffen | Global Change and the Earth System: A Planet Under Pressure[END_REF]. In recent years, invertebrates, and especially arthropods, have been at the centre of a debate [START_REF] Dornelas | Nuanced changes in insect abundance[END_REF][START_REF] Mcdermott | News Feature: To understand the plight of insects, entomologists look to the past[END_REF] regarding the magnitude and even the directionality of the temporal trends in their abundance. Some studies showed a strong decline on the basis of standardized inventories [START_REF] Hallmann | More than 75 percent decline over 27 years in total flying insect biomass in protected areas[END_REF]Seibold et al. 2019), while analyses of meta-datasets assembling heterogeneous time series evidenced a decline of terrestrial, but an increase of aquatic arthropods [START_REF] Pilotto | Meta-analysis of multidecadal biodiversity trends in Europe[END_REF][START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF] or no overall decline [START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF]. Finally, an analysis of opportunistic occurrence data revealed nonmonotonous dynamics and no overall decline [START_REF] Outhwaite | Complex long-term biodiversity change among invertebrates, bryophytes and lichens[END_REF].

Heterogeneity in population trends among studies is not surprising but the underlying causes need to be differentiated, particularly to tell apart abiotic/biotic factors from methodological factors. The dynamics of biodiversity changes are well-documented and their heterogeneity relatively well understood for vertebrates [START_REF] Antão | Temperaturerelated biodiversity change across temperate marine and terrestrial systems[END_REF]Daskalova et al. 2020b, a;[START_REF] Leung | Clustered versus catastrophic global vertebrate declines[END_REF] or for some specific functional groups of invertebrates, such as pollinators [START_REF] Grab | Agriculturally dominated landscapes reduce bee phylogenetic diversity and pollination services[END_REF][START_REF] Duchenne | Long-term effects of global change on occupancy and flight period of wild bees in Belgium[END_REF][START_REF] Soroye | Climate change contributes to widespread declines among bumble bees across continents[END_REF][START_REF] Millard | Global effects of land-use intensity on local pollinator biodiversity[END_REF]. However, the heterogeneity in population trends remains poorly explained for most invertebrates. Several reasons may explain the contrasting patterns revealed by the studies involved in the arthropod decline debate. Global change pressures can vary in space among locations or ecological habitats, and species abilities to respond to environmental changes may depend on their traits or their evolutionary history [START_REF] Helmus | Communities contain closely related species during ecosystem disturbance[END_REF][START_REF] Grab | Agriculturally dominated landscapes reduce bee phylogenetic diversity and pollination services[END_REF]; hence spatial, ecological and taxonomic coverages are obvious sources of heterogeneity among studies and are widely discussed in the recent literature [START_REF] Blowes | The geography of biodiversity change in marine and terrestrial assemblages[END_REF][START_REF] Pilotto | Meta-analysis of multidecadal biodiversity trends in Europe[END_REF]. For example, some studies used global datasets with a bias towards the northern hemisphere [START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF] while others considered national or even more local spatial extents [START_REF] Hallmann | More than 75 percent decline over 27 years in total flying insect biomass in protected areas[END_REF]. Some focused on terrestrial arthropods [START_REF] Hallmann | More than 75 percent decline over 27 years in total flying insect biomass in protected areas[END_REF]Seibold et al. 2019) while others included aquatic groups [START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF][START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF]. Finally, some studies analysed average trends by pooling multiple taxonomic groups [START_REF] Hallmann | More than 75 percent decline over 27 years in total flying insect biomass in protected areas[END_REF][START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF], while others reported trends for each taxonomic group [START_REF] Outhwaite | Complex long-term biodiversity change among invertebrates, bryophytes and lichens[END_REF]. Differences in temporal coverage (i.e. the time period during which a taxon or community was monitored) across studies are at the core of the debate over biodiversity crisis [START_REF] Pauly | Anecdotes and the shifting baseline syndrome of fisheries[END_REF][START_REF] Cardinale | Is local biodiversity declining or not? A summary of the debate over analysis of species richness time trends[END_REF][START_REF] Loreau | Do not downplay biodiversity loss[END_REF], but efforts to account for temporal coverage in statistical analysis remain limited [START_REF] Cardinale | Is local biodiversity declining or not? A summary of the debate over analysis of species richness time trends[END_REF][START_REF] Didham | Interpreting insect declines: seven challenges and a way forward[END_REF][START_REF] Loreau | Do not downplay biodiversity loss[END_REF]. Heterogeneity in temporal coverage is likely to influence the conclusions of studies assessing temporal trends in arthropod abundance for two reasons. First, there is variation in baseline years across the studies fuelling the debate: some analysed time series of identical temporal coverage and relatively old baseline year [START_REF] Hallmann | More than 75 percent decline over 27 years in total flying insect biomass in protected areas[END_REF][START_REF] Outhwaite | Complex long-term biodiversity change among invertebrates, bryophytes and lichens[END_REF], while others used heterogeneous datasets consisting mostly of time series with recent baseline years [START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF][START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF]. Second, non-monotonous dynamics repeatedly reported in arthropods [START_REF] Macgregor | Moth biomass has fluctuated over 50 years in Britain but lacks a clear trend[END_REF][START_REF] Baranov | Complex and nonlinear climate-driven changes in freshwater insect communities over 42 years[END_REF][START_REF] Høye | Nonlinear trends in abundance and diversity and complex responses to climate change in Arctic arthropods[END_REF][START_REF] Schowalter | Arthropods are not declining but are responsive to disturbance in the Luquillo Experimental Forest, Puerto Rico[END_REF], due to temporally variable environmental pressures [START_REF] Baranov | Complex and nonlinear climate-driven changes in freshwater insect communities over 42 years[END_REF][START_REF] Duchenne | Long-term effects of global change on occupancy and flight period of wild bees in Belgium[END_REF][START_REF] Schowalter | Arthropods are not declining but are responsive to disturbance in the Luquillo Experimental Forest, Puerto Rico[END_REF], can affect linear estimates of abundance trend beyond the original description of the shifting baseline syndrome (Fig. 1). With non-monotonous dynamics, estimated linear trends can vary from positive to negative when considering different baseline years for the same time series (e.g. Fig. 1h). As most arthropod trends were estimated assuming linear trends and available data often do not make it possible to account for all factors influencing temporal variation in abundance, it appears critical to account for differences in baseline years when making comparison within and among studies. Here, we evaluated the effect of the baseline years of time series on estimated linear trends in arthropod abundance, using the four largest datasets of arthropod time series recently analysed in studies fuelling the current debate [START_REF] Dornelas | BioTIME: A database of biodiversity time series for the Anthropocene[END_REF][START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF][START_REF] Outhwaite | Complex long-term biodiversity change among invertebrates, bryophytes and lichens[END_REF][START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF]. We first characterized non-monotony in temporal variations of arthropod abundance. Then, using simulations we assessed how non-linearity can bias linear abundance trends. Finally, using a sliding baseline method on the empirical dataset aggregated here (cf. Methods, Fig. 2), we measured how nonlinearity in arthropod abundance trends can produce a statistical dependency of estimated linear abundance trends on the baseline year.

Methods

The four source datasets (step 1)

We aggregated four source datasets from recent publications evaluating abundance trends in arthropods (Table S1): (i) annual occupancy estimates (the proportion of 1km 2 grid cells in a region occupied by a species, a proxy for abundance) at species level for a wide diversity of arthropods from Great-Britain, produced by [START_REF] Outhwaite | Annual estimates of occupancy for bryophytes, lichens and invertebrates in the UK, 1970-2015[END_REF], (ii) annual estimates of arthropod abundances mostly at species level from American Long-Term Ecological Research (LTER) sites from [START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF], (iii) annual estimates of arthropod abundances from the meta-analysis of van [START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF], aggregated at the resolution of taxonomic order, and (iv) abundance time series from the BIOTIME database [START_REF] Dornelas | BioTIME: A database of biodiversity time series for the Anthropocene[END_REF]), mostly at the species level (Fig. S1).

Fig. 2: Schematic description of data aggregation, and workflow of analyses to assess the effect of the temporal baseline on abundance trends. In Step 1, we aggregated empirical datasets from four different source datasets. Then we characterized non-monotony of abundance variations over time (Step 2), and we assessed how non-monotony can affect estimated linear abundance trends using simulations (Step3, Goal 1). Finally, we created an array of subset datasets using a sliding baseline (Step 4) to estimate how baseline years influenced estimated linear abundance trends in empirical data (Step 5, Goal 2). Abundance trends are expressed as growth rates for all time series, allowing comparisons among common and rare species, as well as among datasets (cf. Supplementary Methods, Fig. S2 &S3). Statistical analyses were then performed accounting for habitat (aquatic vs terrestrial), continent, data source and taxonomy.

We focused on well covered continents and habitats, using only arthropod time series from North America and Europe with information on habitat (aquatic vs. terrestrial). We homogenized taxonomy using the Global Biodiversity Information Facility (GBIF) taxonomy backbone. For data from [START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF] and from van [START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF], some time series describe the temporal variations of a wide diversity of species pooled together by summing their respective abundances. For these datasets, we removed all time series with taxonomic resolution coarser than taxonomic order, except for non-insect arthropods that are often grouped at taxonomic class level in available datasets (Chilipoda, Diplododa, Collembola, Branchiopoda etc.). We retained these groups and consider them to be the same rank as taxonomic orders in the following for simplicity. Details about aggregation and filtering step are available in Supplementary Methods.

To study baseline year effects on linear abundance trends estimated from time series with consistent ending dates, we removed time series ending before 2005 (n = 14,717). We further removed the few abundance values before 1970 (n = 1,039, 0.4%) to focus on the time period when most of the data were collected. Finally we removed time series shorter than 3 years (n = 47). This led to 14,130 original time series (Table S2, Fig. S1).

Assessing the monotony of temporal variation in abundance (step 2, goal 1)

Because estimated linear abundance trends depend on the baseline year only if abundance varies nonlinearly, especially non-monotonically, over time (Fig. 1), we estimated non-monotony using a Generalized Additive Model (GAM) for each original time series. We estimated the strength of non-monotony as the number of local extrema, hereafter turning points, observed in the non-linear trend predicted by the GAM, as detailed in supplementary Methods.

Simulations of temporal declines for different population dynamics (step 3, goal 1)

We used simulations to assess how non-linearity in temporal variation of abundance (i.e. population dynamics) can bias estimated linear abundance trends, depending on their baseline year. We used four different shapes of population dynamics (Fig. 4a), using four functions describing temporal variation in average population size over 42 years (cf. R script available in supplementary materials). For each shape of population dynamics, we simulated time series with different growth rates r from stable (r = 1) to declining (r = [0.95, 0.9, 0.85]). The three latter values of r correspond to declines of 5%, 10% and 15% per year, respectively. We simulated 100 abundance time series for each shape of population dynamics and growth rate. To do so, for each year, we sampled observed abundance values from a Poisson distribution with a mean equal to the average population size of the corresponding year. Then, we estimated a linear abundance trend over the entire time period (baseline year at t = 1), as well as over truncated time series with different baseline years (t =10, 20, 30), using a Generalized linear model (GLM) with a Poisson error structure and a log link function, and accounting for temporal autocorrelation (equation (1) below).

Generating an array of subset datasets using a sliding baseline (step 4)

To study the effect of the baseline year of time series on abundance trends estimated from the 14 130 empirical original time series, we created an array of 41 datasets, hereafter called subset datasets, corresponding to 41 different baseline years, from 1970 to 2010 by steps of one year (Fig. 2). For each of the subset datasets, time series were either truncated to start at the given baseline year, or removed if they did not include this specific baseline year. By construction, each of the original time series therefore appears several times in the array of 41 datasets, corresponding to n (1 ≤ n ≤ 41) truncated time series. Since time series with old baseline years are rare, the number of time series included in the subset datasets decreases with earlier baseline years.

Estimating abundance trends (step 5)

We estimated arthropod abundance trends using one GLM per truncated time series. We considered only truncated time series with at least three annual estimates of abundance, including the abundance estimate in the year used to truncate the time series.

To obtain comparable abundance trends among the various sources, expressed in the same unit, we used a model structure that allows the estimation of growth rates. To do so we used GLMs with Poisson error structure with a log link function for count data, from van Klink et al. (2020), [START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF] and BIOTIME, GLMs with a binomial error structure with a logit link function for occupancy estimates from [START_REF] Outhwaite | Annual estimates of occupancy for bryophytes, lichens and invertebrates in the UK, 1970-2015[END_REF], and GLMs with a gaussian error structure with a log link function for density estimates from BIOTIME. Trends estimated with a log or logit link functions are expressed as the logarithm of a growth rate (see Supplementary Methods, Fig. S2 &S3), allowing comparisons between common and rare species, but also among datasets. Therefore, this approach with appropriate link functions gives the same importance to rare and to common species, in contrast to classic standardization of abundance by mean and standard deviations, which biases average trends by giving more weight to species with lower inter-annual variability in abundance (Fig. S4).

The GLM used for each truncated time series explains the abundance estimates of each year 𝑖 by a Poisson or Binomial distribution, of parameter 𝜆 𝑖 and 𝑝 𝑖 respectively, which depend on a linear year effect (𝛽):

𝑙𝑜𝑔(𝜆 𝑖 ) 𝑜𝑟 𝑙𝑜𝑔𝑖𝑡(𝑝 𝑖 ) = 𝛼 + 𝛽 × 𝑦𝑒𝑎𝑟 𝑖 + 𝜑 𝑖 + 𝑒 𝑖 (1)
Evaluating the importance of the effect of baseline year in arthropod abundance trends (goal 2) Pooling together the slopes of the year effect from the linear models presented above, across the 41 subset datasets, we obtained 192,244 abundance trends. We removed abundance trends estimated from time series with a single non-zero yearly estimate of abundance. We did so because growth rates estimated from such time series are likely to be extreme (i.e. strongly negative or positive) if this positive abundance estimate is by chance at the end or at the beginning of the truncated time series (Fig. S5). We thus kept the 175,796 abundance trends derived from truncated time series with ≥2 non-zero abundance estimate.

We examined the effect of baseline year used to truncate the original time series on abundance trends. To this end, we used a Bayesian linear mixed-effects model explaining abundance trends with spatial and temporal variables. Since we do not expect a linear relationship between abundance trends and baseline year (Fig. 1), the effect of baseline year was modelled as a polynomial of order three. We added a continent effect, a habitat effect, and all two-and three-way interactions between these effects and the polynomial baseline year effect (equation ( 2)), because recent results suggest that terrestrial and aquatic arthropods exhibit differences in their abundance trends and because population trends can vary over space [START_REF] Outhwaite | Complex long-term biodiversity change among invertebrates, bryophytes and lichens[END_REF][START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF]. To assess whether differences among source datasets [START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF][START_REF] Outhwaite | Complex long-term biodiversity change among invertebrates, bryophytes and lichens[END_REF][START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF]) persist after spatio-temporal variables have been taken into account, we included a dataset effect, with four levels.

We also accounted for pseudo-replication in species belonging to the same taxonomic orders, as well as for the fact that several truncated time series are obtained from the same original time series, by adding a random effect of taxonomic order and time series ID on the intercept. Data from van [START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF] and [START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF] originated from multiple sites and different data sources. Thus, we accounted for this structure by adding a random site effect nested in a random data source effect. For data from [START_REF] Outhwaite | Annual estimates of occupancy for bryophytes, lichens and invertebrates in the UK, 1970-2015[END_REF], the data source corresponds to the groups used to calculate occupancy estimates in the original dataset and the site corresponds to the country coverage of the data for each group (UK or GB), both extracted from the Online-only Table 1 [START_REF] Outhwaite | Annual estimates of occupancy for bryophytes, lichens and invertebrates in the UK, 1970-2015[END_REF].

The linear mixed-effect model is thus the following:

𝐴𝑇 𝑗𝑏𝑐ℎ𝑠𝑜𝑑𝑙𝑖 = 𝛼 𝑐ℎ + 𝛽 𝑠 + 𝜑 𝑏𝑐ℎ + 𝜃 1 𝑜 + 𝜃 2 𝑑 + 𝜃 3 𝑙 + 𝜃 4 𝑖 + 𝑒 𝑗𝑏𝑐ℎ𝑠𝑜𝑑𝑙𝑖 (2)
where 𝐴𝑇 𝑗𝑏𝑐ℎ𝑠𝑜𝑑𝑙𝑖 is the abundance trend of truncated time series 𝑗 with baseline b, from continent 𝑐, habitat ℎ, dataset 𝑠, order 𝑜, data source 𝑑, site 𝑙, from original time series 𝑖. 𝛼 𝑐ℎ is the intercept for all combination of continent 𝑐 and habitat ℎ, while 𝛽 𝑠 is the effect of the source dataset s. 𝜑 𝑏𝑐ℎ is baseline effect that depends on baseline b, continent 𝑐 and habitat ℎ. It is modelled as a temporal random walk of order one (𝜑 𝑏𝑐ℎ ~𝒩(𝜑 (𝑏-1)𝑐ℎ , 𝜎 𝜑 2 ), with 𝜑 1𝑐ℎ = 0). The 𝜃s denote random effects on the intercept: 𝜃 1 𝑜 for taxonomic order, 𝜃 2 𝑑 for data source, 𝜃 3 𝑙 for site and 𝜃 4 𝑖 for time series ID. 𝑒 𝑗𝑏𝑐ℎ𝑠𝑜𝑑𝑙𝑖 is an error term, independent and identically distributed, following 𝒩(0, 𝜎 2 ).

To account for the fact that our response variable is estimated, and thus each value has an associated standard error (𝑠𝑑𝑒 𝑗 ), we modelled the residual variance as a function of this error. In addition, since the four source datasets have different taxonomic scopes or different spatial scales, we expect that residual variance (𝜎 2 ) will be strongly structured by the source dataset. Finally, the baseline year also affects the number of abundance estimates in time series, which is expected to affect the stochasticity of abundance trends [START_REF] Bahlai | The broken window: An algorithm for quantifying and characterizing misleading trajectories in ecological processes[END_REF]. We modelled the dependence of variance of the residuals (σ²) on standard error associated to abundance trends of truncated time serie j (𝐴𝑇 𝑗 ), on source dataset s and on the number of years with data in each truncated time series (ny):

log(𝜎 2 𝑗𝑠 ) = 𝛼 𝑠 + 𝛽 1𝑠 × 𝑛𝑦 𝑗 + 𝛽 2𝑠 × 𝑛𝑦 𝑗 2 + 𝛽 3 × log (𝑠𝑑𝑒 𝑗 ) (3)
where 𝛼 𝑠 is the intercept, which depends on source dataset s, and 𝛽 1𝑠 and 𝛽 2𝑠 the polynomial effects of the number of years with data in truncated time series j, for each dataset, on the residual variance 𝜎 2 𝑠𝑗 . 𝛽 3 is the effect of standard error associated to 𝐴𝑇 𝑗 on the residual variance 𝜎 2 𝑠𝑗 . Parameters of the two models, one for the mean (equation ( 2)) and one for the residual variance (equation ( 3)), are estimated simultaneously. On this model (equation ( 2)), we estimated the variance of abundance trends explained by each random or fixed effect, as the ratio of its variance on the sum of all these variances plus the variance of the residuals. Priors used are detailed in the R script available in supplementary material. The model was fitted using the R2jags R package [START_REF] Su | Butterfly monitoring in Europe: methods, applications and perspectives[END_REF], with 3 chains using 60,000 iterations with a burnin of 50,000 and a thin rate of 3, which was enough to reach convergence (all parameters with Rhat<1.1).

Results

The number of turning points per empirical time series increases with the number of years with data in the time series, indicating that the vast majority of population trends are non-monotonous when time series are long enough (Fig. 3a). As expected from Figure 1, the strength of the non-monotony, measured as the number of turning points per time series, affects trend estimation. Non-monotonous abundance trends are characterized by multiple changes in trend direction over time, such that estimated linear trends may have opposite signs depending on the baseline year considered as the start of the time series (Fig. 3b). This pattern is consistent across the four source datasets (Fig. S6).

Using simulated population dynamics (Fig. 4a), we show that similar rates of declines lead to different estimated growth rates, depending on the shape of the population dynamics (Fig. 4a-c). Assuming linearity to estimate population trends from non-linear dynamics can induce a strong bias when the growth rate departs from stability (r=1), by either overestimating or underestimating the genuine simulated trends (Fig. 4c). As expected, truncating the time series towards more recent baseline years increases the bias in the estimated growth rate, except in simulations with linear population dynamics (Fig. 4d). In addition to the magnitude of the bias, uncertainty (i.e. dispersion of values around the median) also increases when truncating the time series towards more recent baseline years (Fig. 4d). Short and recent time series tend to produce extreme estimated growth rates, which can even be opposite to the simulated long-term decline (Fig. 4d). These results suggest that not accounting for non-linearity in population trends can induce strong biases, and that short-term and recent time series cannot be used to infer long-term population change, even if population dynamics are linear, due to the strong uncertainties on estimated values.

Fig. 4: Non-linear population dynamics can bias estimated abundance trends. (a)

Raw patterns of population dynamics when the population is stable, i.e. with a growth rate r = 1, or when the population declines by 5% (r=0.95), 10% (r=0.9) or 15% (r=0.85) each year. (b) Scaled patterns of population dynamics shown in (a) but using a log10 y-axis scale illustrating that the red curve corresponds to linear dynamics (cf. Supplementary methods, Fig. S2). (c) Estimated values of log(r) as a function of the shape of the population dynamics over 100 stochastic simulations per type of population dynamics. (d) Estimated values of log(r) as a function of the shape of the population dynamics and of the baseline year used to truncate the time series. In (c) and (d) the dashed horizontal line shows the value of the logarithm of the true (simulated) growth rate. Boxplots have the same meaning as in Figure 3. Right panel of (d) without outliers is represented in Figure S7.

When looking for general patterns in all empirical truncated time series together, statistical analyses show that average linear abundance trends estimated from the truncated time series strongly depend on the baseline year used for truncation, in interaction with habitats and continents (Fig. 5). This suggests that the non-linearity in temporal variations of abundance (Fig. 3b) leads to a strong dependence of the Arthropod population trends to the considered period (Fig. 5a). Results show a higher uncertainty in North America likely because time series with old baseline years are scarce for this region (Fig. 5b). This highlights that data are missing to estimate long-term abundance changes of Arthropods, and that recent data cannot help to fill this blank because of non-linearity. Importantly, for the two source datasets including short and long time series (i.e. those from Outhwaite et al. and van Klink et al.,Fig. 6a), the residual variance of the model strongly decreases with the number of years with data in time series. This is consistent with our previous results based on simulations of truncated time series (Fig. 4c), and suggests that the stochasticity in abundance trends estimated from short time series is much greater than that from long time series. For data from Crossley et al. and BIOTIME, the relationship between the residual variance in abundance trends and the number of years with data is strongly parabolic. However, those datasets are strongly biased towards recent and short time series (Fig. 6a), prompting caution in interpreting this signal.

Although statistically significant, the effects of baseline year explain only a small fraction of the total variation in abundance trends (Fig. 6b). Overall, random effects controlling for methodological issues, spatial and taxonomic heterogeneity explain 56% of the variation in abundance trends (Fig. 6b). A high proportion of this variation in abundance trends is explained by the random effect of time series ID, controlling for the artificial dependencies among abundance trends generated by the truncation procedure (cf. Methods). Most of the remaining variation in abundance trends is explained by taxonomic order and local site random effects, highlighting that abundance trends are strongly heterogeneous among clades and among sites. This questions the relevance of global estimates and stressing the need to carefully identify the drivers of such heterogeneity. Some groups, such as Trombidiformes (mites), Orthoptera, Collembola and Isopoda, exhibit more positive average abundance trends while other groups, such as Archaeognatha (jumping bristletails), Hymenoptera, Coleoptera, Dermaptera and Blattodea are associated with more negative trends (Fig. S8). Strikingly, once the various sources of heterogeneity are controlled for, we do not find any significant difference among the four source datasets from which the abundance estimates were extracted (Fig. 6c).

Discussion

Accounting for differences in baseline years contributes to settling the debate on arthropod decline Our analysis shows that the non-monotony of empirical abundance time series induces a strong dependency of estimated abundance trends to the baseline year, as well as a strong uncertainty in abundance trends from short time series. Our findings therefore bring statistical support to the fact that most of the available data regarding arthropods, which are biased towards recent and short time series, are not appropriate to extrapolate long-term population trends, as suggested recently in discussions over vertebrates and arthropods population trends [START_REF] Daskalova | Accounting for year effects and sampling error in temporal analyses of invertebrate population and biodiversity change: a comment on Seibold et al[END_REF][START_REF] Loreau | Do not downplay biodiversity loss[END_REF][START_REF] Mehrabi | Shifting baselines and biodiversity success stories[END_REF].

We show that baseline year is statistically linked to abundance trends, and we highlight that average linear abundance trends, if they make any sense, should be interpreted in the light of the temporal window covered by the analyzed time series. By statistically assessing the effect of shifting baseline, we provide a general picture of how abundance trends change as function of baseline year, which may help to reposition the findings of past and future studies in a broader context, and hopefully to make between-studies comparison easier. Importantly, once we controlled average temporal trends of studies by the various sources of heterogeneity we did not find any statistical difference among source datasets, suggesting that the contradiction among recent results regarding arthropod decline comes from methodological issues related to temporal coverage as well as taxonomical and geographical bias.

The importance to define temporal baseline and use common spatial yardsticks when evaluating temporal change was previously emphasized in the wider context of biodiversity decline [START_REF] Lotze | Historical baselines for large marine animals[END_REF][START_REF] Mihoub | Setting temporal baselines for biodiversity: the limits of available monitoring data for capturing the full impact of anthropogenic pressures[END_REF][START_REF] Cardinale | Is local biodiversity declining or not? A summary of the debate over analysis of species richness time trends[END_REF]Stouffer et al. 2021), but here we bring further statistical support to this caveat. Indeed, although the importance of comparing results with common baselines is well known, the baseline effect is rarely explicitly accounted for in quantitative analysis (but see [START_REF] Macgregor | Moth biomass has fluctuated over 50 years in Britain but lacks a clear trend[END_REF] for an example), in contrast to other sources of heterogeneity in abundance trends (e.g. space and taxonomy). Although [START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF] tested the effect of the starting year by truncating their time series (see Fig. 3 of van Klink et al. 2020), they did not formally test for a baseline effect: discarding data earlier than a given baseline threshold, from 1960 to 2005, did not result in constraining the baseline years to be equal across all time series. Since their dataset was mainly composed of time series with a baseline year after 1990, discarding data before a giving baseline threshold did not affect the overall distribution of baseline years much, except for recent thresholds (post-1990).

We also show that shorter time series exhibit much more stochasticity (i.e. higher residual variance) than long term series, which increases the uncertainty of results from these short-term time series that are commonly used in assessing arthropod abundance trends (Seibold et al. 2019;[START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF][START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF]. In other words, our results show that short-term series need to be much more replicated than long-term time series to reach the same level of confidence in the results. This is consistent with previous trends assessment regarding moths in Great-Britain [START_REF] Macgregor | Moth biomass has fluctuated over 50 years in Britain but lacks a clear trend[END_REF]) and more generally with the fact that arthropod abundance trends estimated from short time series are strongly sensitive to year to year variations [START_REF] Daskalova | Accounting for year effects and sampling error in temporal analyses of invertebrate population and biodiversity change: a comment on Seibold et al[END_REF]). Here we show that uncertainty decreases exponentially with the length of time series, highlighting the importance to maintain existing biodiversity monitoring schemes. Obtaining long enough time series is critical for assessing reliable abundance trends, which echoes hot and recent debates about abundance trends in vertebrates [START_REF] Leung | Reply to: Do not downplay biodiversity loss[END_REF][START_REF] Loreau | Do not downplay biodiversity loss[END_REF].

Limitations and future challenges

In addition to the dramatic lack of data for some regions of the world, our study suggests that available arthropod abundances suffer from temporal limitations that should be carefully kept in mind when assessing their trends. First, the significant effects of baseline year and study area (continent and local site) suggest that the estimation of the arthropod abundance trends suffers from large uncertainties, mainly due to the lack of historical data. Here, 1970 is the oldest baseline we can tackle with the data at hand, but it cannot be considered as a reference before the rise in global change pressures [START_REF] Mihoub | Setting temporal baselines for biodiversity: the limits of available monitoring data for capturing the full impact of anthropogenic pressures[END_REF], leading to a likely underestimation of global change effects on arthropod abundance. Since we cannot go back in time to sample biodiversity, this issue will remain difficult to solve. Museum collections and other sources of historical data could help fill this gap, although extracting reliable information from such data is still challenging [START_REF] Isaac | Statistics for citizen science: extracting signals of change from noisy ecological data[END_REF][START_REF] Bartomeus | Historical collections as a tool for assessing the global pollination crisis[END_REF][START_REF] Duchenne | Long-term effects of global change on occupancy and flight period of wild bees in Belgium[END_REF][START_REF] Outhwaite | Complex long-term biodiversity change among invertebrates, bryophytes and lichens[END_REF].

Whether or not scientists can manage to obtain the necessary data and apply relevant methods to effectively turn back the clock , our analysis stresses the critical need to maintain long-term monitoring and secure appropriate archiving of related data [START_REF] Millar | The presence of citizen science in sustainability reporting[END_REF].

Second, we used the same datasets as previous studies, so our analysis does not introduce new elements to assess the reliability of these abundance trends with respect to potential biases related to spatial and taxonomic coverage or data quality [START_REF] Desquilbet | Comment on "Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF]. For example, some of the time series incl ded in van lin et al 's and Crossle et al 's datasets ere rod ced b ex eri ental st dies manipulating environmental conditions likely to influence abundance trends [START_REF] Desquilbet | Comment on "Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF][START_REF] Desquilbet | Adequate statistical modelling and data selection are essential when analysing abundance and diversity trends[END_REF], leading to criticisms about the use of these time series to assess temporal trends [START_REF] Cardinale | Is local biodiversity declining or not? A summary of the debate over analysis of species richness time trends[END_REF]. Similarly, using only or mostly data from research stations, such as Long-Term Ecological Research sites, could bias estimated abundance trends upward as these locations are often partially protected from disturbances. Our results show that local site strongly explains heterogeneity in abundance trends, which is consistent with the strong influence of local changes in environmental conditions on arthropod abundance trends (Seibold et al. 2019). This has consequences for the interpretation of global trends obtained from a non-representative sample of sites. This potential bias stresses the need for standardized protocols to monitor arthropod abundance in numerous sites, representative of the areas covered by different habitats and land-uses, to handle the diversity of anthropic pressures on biodiversity, some remaining restricted to particular areas while others apply widely over space. Monitoring schemes based on citizen sciences are one way to tackle this challenge [START_REF] Su | Butterfly monitoring in Europe: methods, applications and perspectives[END_REF][START_REF] Jeliazkov | Large-scale semi-automated acoustic monitoring allows to detect temporal decline of bush-crickets[END_REF], since they can produce protocoled or semi-protocoled datasets over a large set of habitats and landscapes, over seasons and years, at national or even continental scales. Citizen science monitoring schemes are often recent, but would be of considerable help to ensure long-term monitoring of species abundances, should they be maintained in the future.

Moreover, by expressing abundance trends as growth rates, we gave the same importance to rare species as to common species, which could be debated. Decline of extremely rare species could be poorly informative and less likely to affect ecosystem functioning than decline of common species, while rare species can exhibit extreme abundance trends (Fig. S9), thus affecting average abundance trends in a nonnegligible way. On the other hand, rare species can contribute greatly to some biodiversity metrics, such as species richness, phylogenetic diversity or functional diversity, making it relevant to weigh them similarly as common species for some purposes. Consistent with previous comments ( 'Hara & ot e 2010; [START_REF] Desquilbet | Adequate statistical modelling and data selection are essential when analysing abundance and diversity trends[END_REF], we show that transforming abundance counts with log(x+1) before statistical analysis of the data, as done by [START_REF] Crossley | No net insect abundance and diversity declines across US Long Term Ecological Research sites[END_REF][START_REF] Van Klink | Meta-analysis reveals declines in terrestrial but increases in freshwater insect abundances[END_REF], instead of using model structures adapted to the data (i.e. GLM instead of linear models), can introduce an asymmetrical bias, flattening the abundance trends of rare species (Fig. S2). In a similar way, standardizing data by dividing by standard deviation also strongly biases the relative values of abundance trends among species (Fig. S4), stressing the need to explicitly test the effect of any transformations performed on the data.

Finally, the large variation in abundance trends across sites, taxonomic groups, habitats and continents brings into question the relevance of producing global multitaxon linear trends. Global multitaxon trends are likely to be disconnected from the ecological causes and consequences of biodiversity changes. Losses in one place or one taxonomic group cannot be balanced by gains in another place or taxonomic group. Losses and gains can have contrasting ecological and evolutionary consequences, that need to be assessed at relevant ecological scales, e.g. at community level or for a given functional group. Moreover, our results show that variation in arthropod abundance over time is non-linear and sometimes non-monotonous. This suggests that the use of linear analyses is inadequate, despite being the most straightforward and used analysis, and should at least always be associated to the temporal coverage of the data. Disconnecting arthropod decline assessments from temporal yardsticks can affect the understanding of published results making them apparently contradictory. This is particularly important for topics of interest for the general public such as arthropod decline, as it could lead to undermining trust in science [START_REF] Dornelas | Nuanced changes in insect abundance[END_REF].

Fig. 1 :

 1 Fig. 1: Schematic examples of how baseline year may affect estimated linear abundance trends. (a-d) Hypothetical time-series showing abundance variations over years (black) and three possible temporal coverages with contrasting baseline years: old (blue), intermediate (red) and recent (yellow). (e-h) Linear abundance trends estimated with time series starting from these three baseline years. (i-l) Pattern of estimated abundance trends against baseline year, reconstructed by estimating abundance trends over all possible baseline years. The dashed black line shows the zero-value delimiting positive (above) vs. negative (below) abundance trends.

Fig. 3 :

 3 Fig. 3: Empirical time series of arthropod abundance are non-monotonous, leading to unstable abundance trends over a gradient of baseline years. (a) Distribution of the number of original time series as a function of the number of turning points (a proxy for non-monotony, cf. Methods) and of the number of years with data in the corresponding time series. (b) Variation in direction of abundance trends (positive vs negative) of truncated time series, as a function of the number of years with data in the corresponding time series. Truncated time series are from the same original time series which is truncated for every possible baseline year (cf. Methods, Fig. 2). Boxplots represent minimum and maximum values (bottom and top of vertical lines), first and third quartiles (Q1 and Q3, bottom and top of boxes) and median (thick horizontal lines); colours indicate sample size (number of original time series). Points with values outside of the range [Q1-1.5(Q3-Q1), Q3+1.5(Q3-Q1)] are considered as outliers and represented as full circles.

Fig. 5 :

 5 Fig. 5: Effect of the baseline year on average abundance trends. (a) Relationships between the baseline year of times series and estimated abundance trends (log of the growth rate) relative to the abundance trends estimated with the oldest baseline of the dataset (1970), for each continent and habitat. Error bars are 95% confidence intervals. (b) Number of truncated time series used for each baseline value.

Fig. 6 :

 6 Fig. 6: Variation in abundance trends is strongly structured by methodological effects. (a) Predicted (lines) and observed (points) residual variance of the model for each source dataset. Predictions are from equation (3). The histogram at the top shows the distribution of truncated time series along the x-axis, per dataset. (b) Variation in abundance trends explained by the model effects, and residual variation. (c) Black dots show the effect (± CI95%) of the source dataset on abundance trends, relatively to the reference level (BIOTIME).

where 𝛼 is the intercept, 𝜑 𝑖 is a temporal random walk of order one (𝜑 𝑖 ~𝒩(𝜑 𝑖-1 , 𝜎 𝜑 2 ), with 𝜑 1 = 0) to account for temporal autocorrelation and 𝑒 𝑖 an error term. We fitted these GLMs using the INLA R package[START_REF] Rue | Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations[END_REF]).

Acknowledgements

We thank all authors of previous studies on the topic, especially van Klink et al. (2020), Crossley et al. (2020), Outhwaite et al. (2019), Daskalova et al. (2021) and Dornelas et al. (2018) for making the data underlying their work available, and for building step by step our knowledge about the way biodiversity change over time. We also thank an anonymous reviewer from a previous submission for his/her insightful comments and opinions on that manuscript. The simulations were performed at the HPCaVe centre at Sorbonne Université. François Duchenne was funded by the European Research Council (ERC) under the E ro ean Union's Hori on 0 0 researc and innovation rogra (grant agree ent N° 7876 8). Version 3 of this article (https://doi.org/10.1101/2022.02.09.479422) has been peer-reviewed and recommended by Peer Community In Ecology (https://doi.org/10.24072/pci.ecology.100098).

Funding

E ro ean Researc Co ncil (ERC) nder t e E ro ean Union's Hori on 0 0 researc and innovation program (grant agreement N° 787638).

Data, scripts and codes availability

All data used here were publicly available (cf. Methods). Scripts used for analyses are available in Supplementary materials (https://doi.org/10.1101/2022.02.09.479422) and here: https://github.com/fduchenne/Controversy-over-the-decline-of-arthropods.

Supplementary material

Supplementary material are available online: https://doi.org/10.1101/2022.02.09.479422

Conflict of interest disclosure

The authors declare no competing interests.


Bassem Fahs
email: fahsb@rpi.edu

Kefei Wu

Walid Aouimeur

Muhammad Waleed Mansha

Christophe Gaquière

Patrice Gamand

Wojciech Knap

Mona M Hella

About 250/285 GHz push-push oscillator using differential gate equalisation in digital 65-nm CMOS

This study presents a push-push oscillator architecture based on differential gate equalisation to enhance the oscillation frequency while providing relatively high output power with ultra-compact layout form factor. The frequency enhancement is derived as a function of the equivalent RLC model of the oscillator's main constituents. The proposed principle is applied to a terahertz oscillator in the 200-300 GHz range to mitigate the excessive substrate and skin effect losses in standard digital 65-nm complementary metal-oxide-semiconductor technology at such high frequencies. The design concept is validated using two single-stage push-push oscillators. The first oscillator shows -8.1 dBm output power at 250 GHz oscillation frequency and -106.8 dBc/Hz phase noise at 10 MHz offset while consuming 76 mW power from 1.5 V DC supply voltage. The chip area is 200 × 250 μm 2 . The second oscillator provides -14.8 dBm output power at 285 GHz and -106 dBc/Hz phase noise at 10 MHz offset with 80 mW power consumption from 1.5 V DC supply. The chip area is 200 × 200 μm 2 .

Introduction

The unique properties for the electromagnetic (EM) waves in the terahertz and sub-mmWave frequency ranges (penetration, absorption, scattering etc.) promise disruptive technological advances in short-distance high data-rate communication, spectroscopy, bio-medical diagnosis, imaging and others. Existing THz systems are mostly bulky, expensive and based on photonic realisations [START_REF] Nagatsuma | THz communication systems'. Optical Fiber Communication Conf. (OFC)[END_REF]. Integrated systems using low-cost standard integrated circuit (IC) technologies have become a requirement to provide compact, portable, cost-effective and adaptive commercial solutions [START_REF] Park | Semiconductor-based terahertz photonics for industrial applications[END_REF]. It remains that most of the used standard IC processes [complementary metal-oxide-semiconductor (CMOS) or bipolar] have their own frequency limitations and increased parasitic losses at such high frequencies, which seriously reduce the output transmission power and impede the THz solutions from moving towards low-cost commercial products.

One approach of generating power at THz frequencies is based on amplifier-multiplier chains with the input taken from a stable low-frequency source to generate the output signal by mixing effects at higher harmonic frequencies as in [START_REF] Muralidharan | A 165-230 GHz SiGe amplifierdoubler chain with 5 dBm peak output power[END_REF][START_REF] Wu | A wideband SiGe BiCMOS frequency doubler with 6.5-dBm peak output power for millimeter-wave signal sources[END_REF][START_REF] Ergintav | An integrated 240 GHz differential frequency sixtupler in SiGe BiCMOS technology[END_REF][START_REF] Lin | A 200-245 GHz balanced frequency doubler with peak output power of +2 dBm[END_REF]. This approach suffers from high DC power and large chip area. In addition, the phase noise of the output is a function of the input signal's phase noise and thus a stable input frequency source is required. Another approach for THz signal generation is direct generation using highfrequency coupled oscillators to extract the higher order harmonics of the fundamental oscillation frequency. A small chip size could be achieved with low DC power and the signal source in this case is self-contained. Push-push [START_REF] Koo | Design and analysis of 239 GHz CMOS push-push transformer-based VCO with high efficiency and wide tuning range[END_REF][START_REF] Kananizadeh | High-power and high-efficiency millimeterwave harmonic oscillator design exploiting harmonic positive feedback in CMOS[END_REF][9] and self-feeding harmonic [START_REF] Han | A CMOS high-power broadband 260-GHz radiator array for spectroscopy[END_REF] oscillators use balanced structures and symmetric layouts to extract the second harmonic. Triple-push [START_REF] Grzyb | A 288-GHz lens-integrated balanced triple-push source in a 65-nm CMOS technology[END_REF][START_REF] Momeni | High power terahertz and millimeter wave oscillator design: a systematic approach[END_REF][START_REF] Tousi | A 283-to-296 GHz VCO with 0[END_REF] and quad-push [START_REF] Shirazi | A quadcore-coupled triple-push 295-to-301 GHz source with 1.25 mW peak output power in 65 nm CMOS using slow-wave effect[END_REF] oscillator topologies have also been utilised to achieve higher oscillation frequencies at the third and fourth harmonics, respectively. However, triple-push oscillators suffer from complexity in generating a symmetric layout while the output power of the fourth harmonic in the case of quad-push oscillators is relatively small. To increase the output power, either multiple core oscillators are coupled [START_REF] Tousi | A 283-to-296 GHz VCO with 0[END_REF] or arranged in an array format [START_REF] Jameson | A 300 GHz wirelessly locked 2 × 3 array radiating 5.4 dBm with 5.1% DC-to-RF efficiency in 65 nm CMOS[END_REF].

In THz and sub-mmWave frequency ranges, substrate losses and skin effects decrease the quality factor of passive devices and varactors making it more challenging to extract the highest power from the transistors within the oscillator while targeting a small layout form factor. This is particularly true for digital CMOS processes which tend to have a thin back end of line. Design examples proposed in the literature include using a slow-wave coplanar waveguide (CPW) which has lower loss than a standard or grounded CPW to implement the tank [START_REF] Tousi | A 283-to-296 GHz VCO with 0[END_REF]. Varactor-based techniques are applied using hyper-abrupt junction varactor to increase the linearity of C-V characteristics or the parasitic capacitance of metal-oxide-semiconductor (MOS) transistor in [START_REF] Muralidharan | A 110-132 GHz VCO with 1.5 dBm peak output power and 18.2% tuning range in 130 nm SiGe BiCMOS for Dband transmitters[END_REF][START_REF] Nikpaik | A 219-to-231 GHz frequency-multiplier-based VCO with ∼3% peak DC-to-RF efficiency in 65nm CMOS[END_REF], respectively, with both oscillators operated in the 100 GHz band.

In this paper, a push-push oscillator architecture based on a differential structure is demonstrated. The proposed approach is capable of achieving a compact form factor through a simplified routing with less parasitics in addition to reduced DC power consumption. We show that the use of differential equalisation between the cross-coupled metal-oxide-semiconductor field-effect transistor gates reduces the equivalent capacitance at the gate level which results in increasing the frequency of oscillation. A frequency enhancement factor (FEF) is derived as a function of the equivalent RLC model of the oscillator design parts. The design is validated using two single-stage push-push oscillators implemented in 65 nm standard digital CMOS technology. The first oscillator shows -8.1 dBm output power at 250 GHz oscillation frequency while consuming 76 mW power from 1.5 V DC supply with 200 × 250 μm 2 chip area. The second oscillator gives -14.8 dBm output power at 285 GHz with 80 mW power from 1.5 V DC supply and 200 × 200 μm 2 chip area.

The paper is organised as follows. In Section 2, the general push-push design principle with differential gate equalisation is presented. In Section 3, the implementation and the different oscillator parts modelling are detailed. Section 4 gives the experimental setup description and the measurement results. Finally, conclusions are drawn in Section 5.

Design principle

A typical oscillator architecture is based on the use of a resonating part and a transconductance gm-core providing the negative resistance to compensate for the tank losses. Fig. 1a shows a differential oscillation mode at a frequency f 0 , where the signal goes through an internal loop without the need to be referenced to the ground (GND). The push-push behaviour can be described as the circulation of the common-mode (CM) signal at a frequency 2f 0 (or generally at even harmonics of f 0 ) through a loop going from GND to the output load and back to GND. The output signal could be either electrical or radiated (output terminated with an antenna). In the case of an electrical signal, a bias-T or a signal blocker at 2f 0 could be used to provide the DC supply voltage (VCC). Transmission lines (TL) use at the intended frequency range is preferred as they present much lower losses than lumped elements, inductor and capacitor-based tanks [START_REF] Wu | A 103-GHz voltage controlled oscillator with 28% tuning range and 4.2 dBm peak output power using SiGe BiCMOS technology[END_REF].

To increase the output power from the oscillator, the conversion from differential mode to CM signal should be optimised. Since part of the CM power at even harmonics (mainly at 2f 0 ) travelling through the gm-core from GND to the top load will be absorbed by the gm-core transistors' gates, this reduces the power delivered to the output. Adding a 'signal blocking' impedance or transmission line before the gates can prevent such losses if the realised impedance transformation goes from capacitive to inductive. This can also be expressed as a function of gate-blocking TL (TLG) parameters, i.e. the characteristic impedance Z 0,TLG and the electrical length φ TLG to maximise the CM power transfer to the output at 2f 0 , as was also reported in [START_REF] Han | A CMOS high-power broadband 260-GHz radiator array for spectroscopy[END_REF]. Different (Z 0,TLG , φ TLG ) can be determined through simulations to satisfy the maximum output power condition. In this paper, the goal is to build on this understanding and also realise a compact layout with short TLG to minimise the associated losses with the series resistance. On the top side, the TLtop transmission lines depicted in Fig. 1b also contribute to the differential-mode circulation. However, they do not influence the frequency enhancement effect which will be introduced later. Further considerations for TLtop sizing are presented in the implementation section.

Based on the previous discussion, a TLG is added to the crosscoupled oscillator structure, as shown in Fig. 1b. The used MOS transistor (MOST) size is selected to get maximum transition frequency (f T ) performance. With f T = g m /2πC in , where g m is the MOS transconductance and C in is its input gate capacitance, f T allows selecting optimum MOST size and bias point settings regardless of the input/output loading conditions. The MOS size was accordingly set to a width of 40 μm, a gate length of 60 nm and folded into 60 fingers. The gate impedance characterisation, represented by Z MOS in Figs. 2a-c The oscillation frequency (f 0 ) is dictated by two parallel paths. The first path is through TLtop and the second one is through the

f 0 = f 1 a 1 / / f 2 a 2 (1)
where a 1 and a 2 are coefficients that can be related to the relative impedances of TLtop and the (L G , C M ) delay cell, respectively, and f 1 and f 2 are the relative frequencies set by the delays through these TLs, such that

f 1 = 1 2 1 2π L G C M (2)
and

f 2 = 1 4 v 2πl top (3)
where l top is the physical length of TLtop and v is the light velocity in the considered material. In this paper, the introduced frequency enhancement effect is relative to the oscillation frequency increase through f 1 . Thus, the following developments will be only in regard to this frequency. Note that the calculated frequency (2) does not take into account the equivalent capacitance at the MOS drain, which would even result in a lower fundamental frequency. As mentioned earlier, opting for longer TLG to present higher impedance or ideally an open circuit at the gate would also increase L G and lower f 0 . This behaviour defines a trade-off between maximising the output power and achieving a high oscillation frequency.

A similar design approach has been demonstrated in [9,[START_REF] Han | A CMOS high-power broadband 260-GHz radiator array for spectroscopy[END_REF][START_REF] Shirazi | A quadcore-coupled triple-push 295-to-301 GHz source with 1.25 mW peak output power in 65 nm CMOS using slow-wave effect[END_REF]. In [9], the gate blocking is applied on eight single-ended stages to extract the second harmonic through direct power combining. In [START_REF] Han | A CMOS high-power broadband 260-GHz radiator array for spectroscopy[END_REF], the gate blocking is employed in eight differential oscillators each generating the second harmonic along with radiated power combining. In [START_REF] Shirazi | A quadcore-coupled triple-push 295-to-301 GHz source with 1.25 mW peak output power in 65 nm CMOS using slow-wave effect[END_REF], the gate-blocking impedance is used in a differential structure and connected in an oscillator loop of four cells. In this case, the fourth harmonic (4f 0 ) of the oscillation frequency is exploited to generate the output signal (quadruple push oscillator). In all these cases, having a large number of stages gives higher routing complexity, more parasitics, higher DC power consumption and an increasing chip area.

The proposed oscillator design is based on a differential structure to achieve compact form factor, simplified routing with less parasitics and reduced DC power consumption. A differential transmission line (TLd) is added between the cross-coupled transistor gates. Its equivalent series inductance (L D ) cancels part of the gate capacitance (C M ) in addition to the TLd's own distributed parasitic capacitance (C D ), as shown in Figs. 3a andb. This effect could be seen as a differential equalisation where by reducing the equivalent capacitance at the gate level, a higher oscillation frequency is obtained. Note that TLd has a negligible effect on the common mode behaviour as long as its own parasitic capacitance (C D ) can be made negligible compared to the equivalent gate capacitance (C M ).

The transmission line (TLd), assumed to be lossless, can be represented by a double section C D L D C D π-model, as depicted in Fig. 3b. The modified delay cell given in Fig. 3c includes L D and C D elements in parallel to C M , and refers to a half-length of TLd with respect to the middle AC-ground. The parallel combination of L D , C D , and C M is represented by an equivalent capacitance C eq in the delay cell of Fig. 3d, such that

C eq = C M + C D - 1 4π 2 f 2 L D (4)
From ( 2), the new frequency (f′ 1 ) could be written as

f ′ 1 = 1 2 1 2π L G C eq = 1 2 1 2π L G C M + C D -(1/4π 2 f ′ 1 2 L D ) (5)
Rearranging (5), we get

f 1 ′ = 1 + 4L G /L D 1 + C D /C M f 1 = FEF ⋅ f 1 (6)
where f 1 is the initial frequency calculated in [START_REF] Park | Semiconductor-based terahertz photonics for industrial applications[END_REF]. From ( 6), it can be inferred that a FEF between f′ 1 and f 1 has resulted from the proposed equalisation method. For a given TLG design, the MOST size, for a higher FEF would require lowering both L D and C D which leads to a relatively short TLd. Sizing TLd width is a compromise between several factors. For example, finding the minimum series resistance (giving wide TLd), maximising FEF, and setting the right characteristic impedance. Thus, an acceptable value for TLd could be determined through simulations.

To illustrate the potential of the proposed technique, Fig. 4 gives the evolution contours of FEF as a function of L G /L D and C D /C M ratios. As C D can be considered negligible or lower than C M , this lets us consider further the left side in Fig. 4. A large range of L G /L D settings can thus be chosen to easily obtain a FEF larger than one, and this is mainly due to the factor 4 in the numerator in [START_REF] Lin | A 200-245 GHz balanced frequency doubler with peak output power of +2 dBm[END_REF]. If we assume, for instance, that TLd has similar characteristics to the previously considered TLG (L G ≃15 pH, C G ≃6 fF), the resulting FEF gets roughly equal to two.

Implementation and modelling

The oscillator is implemented in a standard eight-metals 65 nm digital CMOS technology from ST Microelectronics. The top metal height is about only 5.7 μm from the substrate, which means that the parasitic capacitances (estimated to ∼20 aF/μm 2 between M8 and M1) is quite high at the target frequency range. The oscillator design uses TL exclusively to avoid the power loss associated with degraded quality factors of passive elements. Fig. 5a gives a descriptive schematic of the oscillator design and Fig. 5b shows a layer-level schematic illustrating TLG and TLd layout and connections to the MOSTs (T 1 , T 2 ).

The modelling of the different oscillator circuit components is critical to accurately predict the oscillation frequency and the output power. This might be challenging as the routing of the MOSTs and TLs has many superposed layers and shared zones to isolate distributed and close mutual capacitances or inductances [START_REF] Beattie | On-chip induction modeling: basics and advanced methods[END_REF].

The MOSTs are sized to 40 μm width, 60 nm length and folded into 60 fingers. The substrate ring in M1 (part of the default library layout) is removed on one side and continued only with active (P+) layer to allow direct connection (with M1 and M2) from the gate to both TLG and TLd. This is done to achieve a minimum series resistance at the cost of slightly increased additional capacitance to the bulk. The source connection is made through parallel M2 and M3 layers to the circuit ground (GND) as shown in Fig. 5b. The capacitive parasitics due to the MOSTs' routing is extracted using Calibre software. The overlapping vias and metal-layer series resistance connected to either the drain or the source is simulated with Momentum software to take into account the resistance increase by the skin effect. The skin depth is estimated to be approximately 0.13 μm at 300 GHz. The extracted RC elements for each routed MOST are represented in Fig. 6a. To provide a minimum ohmic ground to the gm-core despite the very tight skin depth, two wide planes in M1 and M2 are connected over the complete circuit width (∼ 200 μm) from the GND pads to the gmcore cell.

The TLG design is based on a compact and direct connection between the respective MOSTs drains and gates to minimise the inner spacing between transistors to around 7 μm. This avoids the CM power loss by asymmetrical connections caused by lines crossing. The spacing between both TLG microstrips is set to less than 4 μm and each trace is 3 μm wide (Fig. 5b). The TLG layout is made with stacked metals, i.e. a M2/M3 stack starting from the gate to reach the opposite drain with a M4/M7 stack to minimise the series losses with reduced traces dimensions and form factor. The lumped RLC model for the dual TLG lines is shown in Fig. 6b. This model is extracted using S-parameter curves fitting with a four-port Momentum layout. The model is split into dual sections for improved accuracy, with k f representing the coupling factor between the respective line inductances [START_REF] Kopcsay | A comprehensive 2-D inductance modeling approach for VLSI interconnects: frequency-dependent extraction and compact circuit model synthesis[END_REF].

The TLd is placed on the bottom side of the gm-core to provide a direct access to the MOS gates in M1/M2 layers as shown in Fig. 5b. The TL is laid out in M8 layer over M1 GND and sized to 60 μm total length and 8 μm width. The relative RLC model (see Fig. 6c) is extracted using S-parameter curves fitting with a twoport Momentum layout. The TLd line presents a characteristic impedance of 44 Ω.

While the RLC modelling helps define the initial design, EM simulation (done with Momentum) remains unavoidable to take into account the mutual capacitances and magnetic coupling between the different design components. The EM simulation should be extensive and includes all elements without compromising the meshing accuracy or increasing excessively the simulation time. In this design, EM simulations are done on two main parts: the first comprises the metal connections above the MOSTs, TLd, TLG, along with the surrounding ground plane as shown in Fig. 6d, and the second regards the top TLs (TLtop) connection to the ground-signal-ground measurement pads as shown in Fig. 7a. The TLtop model was included as two separate parts in the oscillator schematic (Fig. 7a) due to negligible mutual coupling between both lines with grounded sidewalls and bottom-M1 waveguide structure [START_REF] Elgamel | Efficient shield insertion for inductive noise reduction in nanometer technologies[END_REF]. This also reduces the simulation time needed for optimisation due to the large layout area that would have taken a single EM-model for a dual TLtop layout. Each TLtop is sketched with a 12 μm wide M8 layer over 165 μm long track from the MOST to the signal pad. Shielding sidewalls are formed with M8-to-M1 metal stack and set to 13 μm width and 6 μm spacing away from the signal track. The TL gives around 40 Ω characteristic impedance.

Starting with λ/2 length at 2f 0 for TLtop lines provides 50-Ω CM impedance for Z TOP (looking upwards in Fig. 7a). With an inductive CM impedance presented by the entire bottom part of the oscillator (Z BOT in Fig. 7b), a matching issue comes up with a nonoptimum CM power transfer situation. A solution would consist of shortening TLtop length to allow Z TOP rotation on the Smith chart (see Fig. 7b) in order to place it on the conjugate side of Z BOT for optimum CM matching. The counterpart of this design change is the partial contribution of TLtop to the oscillation frequency (f 0 ) by creating a parallel differential path to the gm-core's one. Managing both constraints in terms of CM impedance matching and oscillation frequency sets TLtop length to slightly less than λ/4 (or ∼165 μm) for a final simulated frequency at 218 GHz and -6.1 dBm output power under 1.5 V DC supply voltage. The partial CM matching in this case is shown in Fig. 7b. This oscillator version is termed (OSC1). Another version (OSC2) of this oscillator design has a shortened TLtop length of 70 μm to assess the effect on the frequency shift. TLtop is modified without any matching optimisation effort and with an expected power performance penalty. The simulated oscillation frequency for OSC2 showed around 250 GHz with about -9 dBm output power from 1.5 V DC supply voltage. The main goal of this VCO design is to illustrate the fundamental frequency enhancement effect while targeting a very small layout form factor. Additional effort is needed to increase the output power via on-chip or spatial combining.

Measurement results

Micrographs of the implemented oscillators OSC1 and OSC2 are shown in Figs. 8a andb, respectively. Both oscillators were directly probed on chip. The effective area is 200 μm × 250 μm and 200 μm × 200 μm for OCS1 and OSC2, respectively.

The testing setups are given in Figs. 9a andb for power and frequency measurements, respectively. Erickson PM5 power meter from Virginia Diodes Inc. (VDI) was used for power measurements, while VDI WR-3.4 SAX down-converter in addition to Keysight N9030A signal analyser was employed for frequency and phase-noise measurements. Both oscillators were tested in free-running mode without any frequency locking.

Measured oscillation frequencies for both oscillators showed an upward shift compared to simulations. OSC1 has an oscillation frequency from 248 to 251 GHz for a supply voltage (VDD) varied from 1.2 to 1.5 V as shown in Fig. 10a. This is compared to 218 GHz from simulations. The deviation in frequency can be due to conservative modelling (mainly with the MOST parasitics) and to the difficulty of separating shared or overlapping connections as well as too close mutual capacitances or magnetic coupling that would have been accounted multiple times. OSC2 presents about Power measurements are corrected with the different setup losses summarised in Table 1. The resulting output power figures after losses extraction are depicted for both oscillators in Figs. 11a andb as a function of VDD. OSC1 maximum output power reached -8.1 dBm at VDD = 1.5 V. OSC2 showed -14.8 dBm output power under 1.5 V. DC power consumptions for both oscillators versus VDD are also given in Figs. 11a andb with maximum 76 and 80 mW for OSC1 and OSC2, respectively. Output power deviations between simulated and measured results could also be explained by the MOST modelling as part of the combined wires resistance was based on extrapolating the Calibrefound resistance to include the skin effect. Other reasons could be due to some over or under estimated shared connections or intermetal layers vias or deviation from transistor parameters postfabrication.

Phase noise measurements are carried out using the same measurement setup with the N9030A signal analyser. OSC1 gives at 10 MHz offset around -104 dBc/Hz and -106.8 dBc/Hz under 1.2 and 1.5 V supply voltage, respectively. OSC2's phase noise is around -106 dBc/Hz at 10 MHz offset for both 1.2 and 1.5 V VDD. Measured phase noise curves for OSC1 and OSC2 at VDD = 1.5 V are given in Figs. 12a andb, respectively. A plateau is present up to a few MHz offset due to the unlocked oscillation frequency. A slope between -30 and -20 dB/decade starts appearing from 10 MHz offset.

This work is compared with the state of the art of THz oscillators implemented in 65-nm CMOS technology in the 250-300 GHz range in Table 2. The presented single-stage oscillators in our case can be distinguished by the low DC power consumption but also by the lowest form factor. The proposed architecture can be evolved into a coupled-array source to increase the output power while keeping a low DC-power level as well as a small form factor advantage. This is an essential characteristic for large-scale THz systems integration in low-cost CMOS processes.

Conclusion

In this paper, a push-push oscillator architecture based on differential gate equalisation to enhance the oscillation frequency while achieving an ultra-compact design form factor has been presented. The frequency enhancement is derived as a function of the equivalent RLC components of the oscillator's design parts.

The design concept has been proven on two single-stage push-push oscillators implemented in standard digital 65-nm CMOS technology. Measured oscillation frequencies showed around 250 and 285 GHz with output power levels around -8.1 and -14.8 dBm, respectively. Phase noise levels are at -106.8 and -106 dBc/Hz at 10 MHz offset. The achieved design area below 200 × 250 μm 2 represents a key feature for potential large-scale integrated coupled-oscillator arrays in low-cost and digital CMOS technologies. b Eight single-ended second-harmonic (2f 0 ) VCOs.

c Eight coupled differential self-feeding second-harmonic (2f 0 ) VCOs.

d Two mutually coupled single-ended triple-push (3f 0 ) oscillator cores.

e Quad-core-coupled triple-push (3f 0 ) VCOs.

f Four coupled differential fourth-harmonic (4f 0 ) VCOs.

  , gives Z MOS = 4 - j12 at 300 GHz, which can be modelled by an equivalent capacitance C M ≃ 48 fF. The effect of TLG can hardly render infinite input impedance from the gate level. However, making the transformed impedance (Z′ MOS ) inductive would relax the CM impedance matching. Setting TLG to an electrical length of roughly 30° and 50-Ω characteristic impedance is equivalent to a series inductance of L G ≃15 pH. Considering either TLG or L G would sensitively give the same impedance rotation starting from Z MOS on the Smith chart shown in Fig. 2c. This also means that the parasitic capacitance for TLG, i.e. C G ≃6 fF (calculated with Z c = (L G /C G ) 1/2 = 50 Ω) can be neglected with respect to C M . The path through TLG to the MOST can thus be considered as a delay cell formed basically by L G and C M .

Fig. 1

 1 Fig. 1 Push-push oscillator (a) With a traditional cross-coupled active gm-core, (b) With added TLG

Fig. 3

 3 Fig. 3 Oscillator design represented (a) Oscillator design with added TLd, (b) Equivalent dual-section π-model for TLd assumed lossless, (c) New delay cell with added elements L D and C D representing a half-length of TLd with respect to the middle AC-ground, (d) Delay cell sketched with the equivalent capacitance C eq of (L D , C D , C M )

Fig. 4 Fig. 5

 45 Fig. 4 FEF contours as a function of of L G /L D and C D /C M ratios varied from 0.1 to 10 each

Fig. 6

 6 Fig. 6 Extracted lumped RLC model (a) MOST routing, (b) Dual-section model for TLG, (c) TLd, (d) EM-simulated Momentum layout including TLd, TLG, the top MOSTs' connections, and the surrounding GND plane

Fig. 7

 7 Fig. 7 ADS simulation showing (a) The used schematic with the different EM-simulated-based models, (b) Partial CM matching between the bottom oscillator impedance (Z BOT ) and the top part impedance (Z TOP ) at 218 GHz frequency

Fig. 8 Fig. 9

 89 Fig. 8 Implemented oscillator micrographs (a) OSC1, (b) OSC2

Fig. 10

 10 Fig. 10 Measured oscillation frequency (a) Versus VDD for both oscillators, (b) OSC1 at VDD = 1.5 V, (c) OSC2 at VDD = 1.5 V

Fig. 11 Fig. 12

 1112 Fig. 11 Measured output power (after losses extraction) and DC power consumption versus VDD (a) OSC1, (b) OSC2

Table 1 Setup

 1 power losses (in dB)

	Circuit	Pad I325_DL2CD (probe)	WR3/WR10 taper (3 cm)	WR10 WG (5 cm)
	OSC1 (250	0.5	3.0	0.35	0.375
	GHz)
	OSC2 (285	0.5	3.4	0.385	0.39
	GHz)

Table 2

 2 Published 250-300 GHz VCOs in 65-nm CMOS Ref.

		Freq., GHz	P out , dBm	P DC , mW	Phase noise, dBc/Hz	Area, µm × µm
	[7] a	239	-4.8	18.5	-98.43 to -110.9 @ 10 MHz	400 × 450
	[9] b	256	4.1	227	-94 @ 1 MHz	670 × 650
	[10] c	260	0.5	800	-78.3 @ 1 MHz	1500 × 1500
	[11] d	288	-1.5	275	-87 @ 1 MHz	500 × 570
	[13] e	298	0.9	235	-79 @ 1 MHz	290 × 316
	[14] f	290	-1.2	325	-78 @ 1 MHz	600 × 600
	this work (OSC1)	248	-9	76	-106.8 @ 10 MHz	200 × 250
	this work (OSC2)	284	-15	80	-106.47 @ 10 MHz	200 × 200

a Single transformer-based push-push (2f 0 ) VCO.

Acknowledgment

The authors would like to thank Keysight for assistance with testing.