Grzegorz Zięba

Carl Smith

Michael G Fox

Stan Yavno

Eva Záhorská

Miroslaw Przybylski

Gérard Masson

Julien Cucherousset

Hugo Verreycken

Hein Van Kleef

| Carl Smith

Mirosław Przybylski

| Hein

| Gordon

H Copp

Gordon H Copp
email: gordon.copp@cefas.co.uk

Red operculum spots, body size, maturation and evidence for a satellite male phenotype in non-native European populations of pumpkinseed Lepomis gibbosus

Keywords: alternative mating strategy, carotenoid, GLMM, sexual selection, visual signal, Centrarchidae

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

| INTRODUCTION

Visual signals play a key role in the environmental biology of teleost fishes [START_REF] Marshall | The visual ecology of reef fish colours[END_REF][START_REF] Wootton | Ecology of Teleost Fishes[END_REF]. Pigmentation is an important component of visual signals in fishes, both within and among species. The expression of colour, which often comprises multicomponent signals, involves pigment-based and structural colours [START_REF] Price | Pigments, patterns, and fish behavior[END_REF][START_REF] Wootton | Reproductive Biology of Teleost Fishes[END_REF]. Carotenoid-based colouration, which typically appears as yellow, orange or red, has been of particular interest because these are costly to produce and thereby serve as a reliable indication of individual condition and foraging ability [START_REF] Olson | Costly sexual signals: are carotenoids rare, risky or required?[END_REF], with a potential role in mating system evolution [START_REF] Wootton | Reproductive Biology of Teleost Fishes[END_REF]. Different carotenoid pigments are responsible for conveying different information to potential mates or rivals [START_REF] Mclennan | The umwelt of the three-spined stickleback[END_REF]. Colour traits can function in sexual selection; selection that acts on differences in fitness among individuals that arise due to the number and identity of their mates. Because sexual selection tends to act more strongly on males than females [START_REF] Andersson | Sexual Selection[END_REF], it is males that usually express the brightest and most elaborate visual traits, including nuptial colouration. Though a broad simplification, there is a tendency in many taxa for a preference to mate with more colourful and ornamented individuals [START_REF] Andersson | Sexual Selection[END_REF][START_REF] Price | Pigments, patterns, and fish behavior[END_REF].

Understanding the adaptive value of mate choice (sensu [START_REF] Kokko | The evolution of mate choices and mating biases[END_REF] based on visual signals is a controversial field of research, but may arise through fitness benefits associated with mate choice, or an arbitrary mechanism [START_REF] Andersson | Sexual Selection[END_REF][START_REF] Wootton | Reproductive Biology of Teleost Fishes[END_REF].

Pigmentation in North American sunfishes (Centrarchidae), in common with most other teleosts, is more extensive in males than in females [START_REF] Scott | Freshwater Fishes of Canada[END_REF]. In the pumpkinseed Lepomis gibbosus males have a striking, red operculum spot that is associated with male aggression [START_REF] Stacey | Body color pattern and the aggressive behavior of male pumpkinseed sunfish (Lepomis gibbosus) during the reproductive season[END_REF]. In aquarium studies, breeding males were shown to attack fish models that intruded into their nesting area, and models with red colour on the operculum spot or iris received more attacks and aggressive displays than those without colour [START_REF] Stacey | Body color pattern and the aggressive behavior of male pumpkinseed sunfish (Lepomis gibbosus) during the reproductive season[END_REF]. It is recognised that dominance relationships among sunfishes are related to colouration, with dominant, brightly-coloured parental (territorial) pumpkinseed males building and defending nests, and performing courtship behaviour to attract females to their nest for spawning. Male sunfishes also adopt alternative mating strategies (sensu [START_REF] Wootton | Introduction: tactics and strategies in fish reproduction[END_REF], whereby subdominant 'cuckolder' males attempt to intrude into the nest of a parental male that is in the act of spawning [START_REF] Gross | Cuckoldry in sunfishes (Lepomis: Centrarchidae)[END_REF][START_REF] Gross | Sneakers, satellites and parentals: polymorphic mating strategies in North American sunfishes[END_REF]. Cuckolder males are early-maturing, do not perform courtship and avoid fitness costs associated with parental care [START_REF] Smith | The costs of parental care in teleost fishes[END_REF]. Small cuckolders, typically possessing disproportionately large testes, are referred to as 'sneakers' due to their rapid entry into the nest of a spawning territorial male to release sperm at the moment of spawning. Large cuckolders, which have testes of intermediate size relative to parental and sneaker males, are termed 'satellites' and participate in spawning by mimicking female characteristics and entering the nest of a territorial male as a potential mate [START_REF] Gross | Sneakers, satellites and parentals: polymorphic mating strategies in North American sunfishes[END_REF][START_REF] Garner | Alternative male reproductive tactics drive asymmetrical hybridization between sunfishes (Lepomis spp.)[END_REF]. Despite aggressive defence of their nest by parental males, cuckolder males frequently intrude successfully during spawning events [START_REF] Rios-Cardenas | The use of Molecular Genetic Analyses to Examine Mating Patterns a Paternal Care in the Pumpkinseed Sunfish[END_REF]. In the close congener, bluegill Lepomis macrochirus, a model taxon for research on alternative reproductive strategies [START_REF] Taborsky | Alternative reproductive tactics in fish[END_REF], [START_REF] Gross | Sneakers, satellites and parentals: polymorphic mating strategies in North American sunfishes[END_REF] characterised two cuckolder phenotypes, based on a population in Lake Opinicon, Canada. Sneaker bluegill are small and have relatively large testes relative to parental males, with a mean gonado-somatic index (GSI) of 4.6 (contrasting with parental males with a mean GSI of 1.1). In contrast, satellite male bluegill are larger than sneakers, with GSI values (mean 3.3) that are intermediate between those of parental males and sneakers. In applying these criteria to pumpkinseed from the same lake, [START_REF] Gross | Sneakers, satellites and parentals: polymorphic mating strategies in North American sunfishes[END_REF] identified the sneaker phenotype, but not that of the satellite. In a subsequent paternity study, [START_REF] Neff | Temporal variation in cuckoldry and paternity in two sunfish species (Lepomis spp.) with alternative reproductive tactics[END_REF] noted that the satellite strategy had yet to be reported for pumpkinseed. In a subsequent genetic study of bluegill-pumpkinseed hybridisation, [START_REF] Garner | Alternative male reproductive tactics drive asymmetrical hybridization between sunfishes (Lepomis spp.)[END_REF] did not distinguish between satellites and sneakers, referring to them collectively as cuckolders. Thus, there are no reports that the satellite phenotype (sensu [START_REF] Gross | Sneakers, satellites and parentals: polymorphic mating strategies in North American sunfishes[END_REF] is expressed in pumpkinseed and hereafter we use the term 'cuckolder' to refer to the alternative male mating strategy in pumpkinseed to the parental strategy.

The approach, onset and end of the reproductive period in fishes can often be identified by variation in GSI, the relationship between gonad weight and somatic weight [START_REF] Wootton | Ecology of Teleost Fishes[END_REF]. Male mating polymorphisms, particularly the cuckolder strategy, are associated with sperm competition, which occurs when the spermatozoa of different males compete for the same eggs [START_REF] Parker | Sperm competition and its evolutionary consequences in insects[END_REF]. Males expressing alternative mating strategies tend to experience a higher "risk" of sperm competition; i.e. an elevated probability that their sperm will compete with the sperm of other males, because they usually mate in the presence of at least one competitor, typically the territorial male and often other cuckolder males [START_REF] Wedell | Sperm competition, male prudence and sperm-limited females[END_REF][START_REF] Wootton | Reproductive Biology of Teleost Fishes[END_REF]. An outcome is that cuckolders are predicted to show adaptations for sperm competition, including relatively (or even absolutely) larger testis size, larger ejaculates, and faster swimming spermatozoa [START_REF] Montgomerie | Testes, sperm, and sperm competition[END_REF]. Consequently, individual GSI can be used to identify male mating strategy.

Pumpkinseed have been widely established outside their native North American range.

Most introductions have been in Europe, where the species was introduced in the late 19 th and early 20 th centuries, establishing populations in at least 28 countries [START_REF] Copp | Growth and life history traits of introduced pumpkinseed (Lepomis gibbosus) in Europe, and the relevance to its potential invasiveness[END_REF]. Most of the research undertaken on pumpkinseed in Europe has focused on female growth and life-history traits (reviewed in [START_REF] Fox | Old world versus new world: life-history alterations in a successful invader introduced across Europe[END_REF], with relatively few studies on males (e.g. [START_REF] Almeida | Cuckoldry features of introduced pumpkinseed sunfish (Lepomis gibbosus) in contrasting environmental conditions in southern Europe[END_REF][START_REF] Valente | Seasonal gonadal development and age-related maturity patterns of introduced pumpkinseed (Lepomis gibbosus Linnaeus, 1758) in a heated thermal reservoir and an adjacent river reach[END_REF]. In its native range, [START_REF] Neff | Temporal variation in cuckoldry and paternity in two sunfish species (Lepomis spp.) with alternative reproductive tactics[END_REF] observed temporal variability in the patterns of pumpkinseed paternity and suggested that these reflect temporal variation in mating opportunities, parental male defence ability, or cuckolder densities. Given this temporal variability, spatial variations are also feasible, possibly extending to visual signals, such as the operculum spot, which is assumed to be important for males for attracting mates [START_REF] Stacey | Body color pattern and the aggressive behavior of male pumpkinseed sunfish (Lepomis gibbosus) during the reproductive season[END_REF].

The aim of the present study was to determine whether the presence and size of the red operculum spot serves as a signal of sex and/or mating strategy in pumpkinseed. To do this we used estimates of GSI to categorise males as nesters or cuckolders and fitted models to test whether spot presence and size varied between the sexes and between males adopting different reproductive strategies. We predicted that parental males, which experience intrasexual contests and intersexual mate choice, would be more likely than cuckolder males and females to express a red operculum spot and that these males would also express larger red spots.

| Methods

Previously unpublished data for 1118 pumpkinseed were obtained from established nonnative populations in Belgium, England, France, the Netherlands and Poland (Table 1), with all data elements being available for 924 specimens of which 465 were males, 387 were females, and 266 were juveniles (Table 2). The fish were collected by a variety of sampling gear (Table 1), depending upon location characteristics and gear availability.

Note that sampling of the two Netherlands populations was incorrectly reported as by funnel trap in [START_REF] Cucherousset | Life-history traits and potential invasiveness of introduced pumpkinseed Lepomis gibbosus populations in northwestern Europe[END_REF]. Where possible, the sampling aimed to obtain ≈ 100 specimens for growth and life-history trait analysis (e.g. [START_REF] Villeneuve | Interpopulation variation in the growth and life history traits of the introduced sunfish, pumpkinseed Lepomis gibbosus, in Southern England[END_REF][START_REF] Cucherousset | Life-history traits and potential invasiveness of introduced pumpkinseed Lepomis gibbosus populations in northwestern Europe[END_REF] and so was not quantitative. In the three English sites, sampling was by catch-per-unit-effort (CPUE) with minnow traps (see [START_REF] Fox | Seasonal reproductive allocation, localscale variation and environmental influences on life history traits of introduced pumpkinseed (Lepomis gibbosus) in Southern England[END_REF]. The data on female growth and life-history traits, but not operculum spot occurrence or size, have been previously published for all these populations [START_REF] Villeneuve | Interpopulation variation in the growth and life history traits of the introduced sunfish, pumpkinseed Lepomis gibbosus, in Southern England[END_REF][START_REF] Cucherousset | Life-history traits and potential invasiveness of introduced pumpkinseed Lepomis gibbosus populations in northwestern Europe[END_REF][START_REF] Fox | Seasonal reproductive allocation, localscale variation and environmental influences on life history traits of introduced pumpkinseed (Lepomis gibbosus) in Southern England[END_REF][START_REF] Valente | Seasonal gonadal development and age-related maturity patterns of introduced pumpkinseed (Lepomis gibbosus Linnaeus, 1758) in a heated thermal reservoir and an adjacent river reach[END_REF], except the two populations from northwest Poland (Table 1): i) a former clay pit (Brodowski Pond) of ≈ 0.9 ha (maximum depth of 7 m) situated in Szczecin with a water temperature that ranges from 4 °C under thick ice cover in winter to 25 °C in summer; and ii) an artificial channel (4 km long, mean depth of 4.5 m), situated on the lower River Oder near the city of Gryfino that receives heated water effluent from a power plant, resulting in daily water temperature fluctuations of 5-15 °C, with maximum temperatures of 30 °C in summer and 15 °C in winter. There is currently no published study on the North American source populations of pumpkinseed in Europe.

For all populations, pumpkinseed were collected at the beginning of the spawning season (between May and June), immediately killed with a lethal dose of anaesthetic following licenced procedures in the country of collection and cooled to freezing for subsequent processing in the laboratory. After thawing, specimens were photographed (left side only) using a Nikon CoolPix 5000 camera for morphological analysis. The presence of operculum spots (denoted as present or absent) was recorded, including the rare occurrence of deviations from red (i.e. yellow or orange, or reddish shades thereof).

Note that the reflectance spectra, hue or intensity of operculum spots were not measured, nor was the spectral sensitivity of pumpkinseed in the red region of the light spectrum known. The presence of an operculum spot was recorded for both sides of each specimen assessed, but only the left side of the fish (the standard for fish measurements) was used for quantifying the size of the red operculum spot, which was estimated as the surface area of the spot (ImageJ, 1.48v). Fish were measured for standard length (SL, nearest 1 mm) and mass (nearest 1 g) and dissected to determine sex and gonad mass (in mg) when present. GSI was calculated as: 100 ´ gonad mass ÷ total body mass) [START_REF] Wootton | Reproductive Biology of Teleost Fishes[END_REF].

| Statistical analysis

To categorise males as parentals or cuckolders, a frequency distributions of male GSI was plotted by country. The resulting multimodal GSI frequency distributions were decomposed to normal distributions using the Bhattacharya method [START_REF] Bhattacharya | A simple method of resolution of a distribution into Gaussian components[END_REF], with a meaningful separation of normal distributions identified with a separation index (SI) > 2 [START_REF] Gayanilo | A draft guide to the Complete ELEFAN[END_REF]. GSI values below the threshold value were categorised as parental and those equal to or above the threshold value were categorised as cuckolders. Data were pooled by country, rather than population, to maximise the precision with which threshold values could be identified.

To test whether parental males were more likely to express a red operculum spot than other males and females, the presence of a red spot in pumpkinseed was modelled using a Bernoulli Generalised Linear Mixed Model (GLMM). Fish SL and mass were collinear, and mass was subsequently dropped from the analysis. Only a single fish that was classified as a juvenile expressed a red spot, and juveniles were subsequently dropped from the analysis. The distribution for the model response variable, link function and predictor function were specified as:

Spotij ~ Binomial(pij, nij) E(Spotij) ~ nij x pij and var(Spotij) = nij x pij x (1 -pij) hij = b1 + b2 x SLij + b3 x strategyij + populationi logit (pij) = hij populationj ~ N(0, s 2 population)
Where Spotij is the probability of fish i in population j having a red operculum spot in nij independent individuals, which was assumed to follow a binomial distribution with an expected probability (E) of expressing a spot of mean piNi and variance pi x (1-pi), with a logit link function. The logit function ensures the fitted probability of a red spot falls between 0 and 1. The variable strategyij is a categorical covariate with three levels, corresponding with fish mating strategy; female, territorial or cuckolder. The model also contained a linear effect for fish SL (SLij). Exploratory analyses demonstrated differences among populations in the size and expression of red operculum spots. To accommodate this effect in the model, the random intercept populationj was included to introduce a correlation structure between observations for different fish from the same population, with variance spopulation distributed normally and equal to 0.

To test whether parental males expressed larger red operculum spots than other males and females, the size of red spots was modelled on a zero-truncated subset of the data using a Gaussian GLMM, specified as:

log10Spotsizeij ~ N(µij, s 2 ) E(log10Spotsizeij) = µij and var(log10Spotsizeij) = s 2 µij = b1 + b2 x log10SLij x b3 x strategyij + populationi populationj ~ N(0, s 2 population)
Where log10Spotsizeij is the log10 size of red spot on fish i in population j and log10SLij is the log10 standard length of fish i in population j. Expected log10Spotsizeij is µij with variance s 2 . Analysis was conducted using the lme4 library [START_REF] Bates | lme4: Linear mixed-effects models using Eigen and S4[END_REF] in R (R Development Core Team, 2017). Model residuals were examined to ensure assumptions were met and that effects were adequately accounted for by the model.

| RESULTS

Male GSI distributions were multimodal within countries (Figure 1). With the exception of data from the Netherlands, the overlap of normally-distributed GSI groups was low; SI was always >2 and thereby provided meaningful separation of modal groups, even for the Netherlands specimens (separation index = 2.68). Decomposition of GSI distributions showed the threshold GSI value that distinguished parental and cuckolder males varied among countries, from 2.3 (Poland) to 4.2 (England) (Table 2,Figure1). In Belgian, French and Polish populations, separation of male strategies was ambiguous, with Bhattacharya decomposition revealing two thresholds within these data (Figure 1). In each case, the higher threshold was used as the cut-off between male strategies. The probability of expressing a red operculum spot was predicted by male SL, with a significantly greater likelihood of a red spot in larger males (Table 3; Figure 2). The probability of a red operculum spot was also predicted by sex, with both cuckolder and territorial males significantly more likely to express a red spot than females (Table 3; Figure 2). Among pumpkinseed that expressed a red spot, there was a significantly stronger interaction between SL and spot size in males than in females (Table 4; Figure 3), though there was no difference between male strategies classified by GSI (Figure 3). Two specimens were observed to possess a yellow operculum spot, a mature female from Schoapedobbe Pond and a male from Brière Marsh (Table 1).

| DISCUSSION

Our analysis of introduced pumpkinseed populations across Europe demonstrated that the expression of a red operculum spot was more likely in mature males than females, (Table 3), whereas only a single sexually immature juvenile expressed a red spot. Further, a significantly stronger interaction was observed (Table 4) between red spot size and body size in parental and cuckolder males than in females. Taken together, these findings indicate that in pumpkinseed, the presence and size of a red operculum spot signals sexual maturity and correlates strongly with body size in males. This signal, however, does not distinguish between the male mating strategies of parental and cuckolder.

Colour signals can have multiple functions, potentially conveying information to rivals, mates or both. In pumpkinseed, the function of the red operculum spot is equivocal. The present findings are consistent with previous studies in demonstrating that a red spot has a potential role in male-male interactions by eliciting an aggressive response by males [START_REF] Stacey | Body color pattern and the aggressive behavior of male pumpkinseed sunfish (Lepomis gibbosus) during the reproductive season[END_REF]. Notably, during aggressive intra-sexual contests, male sunfishes flare their opercula [START_REF] Colgan | Dynamics of aggression in male pumpkinseed sunfish (Lepomis gibbosus) over the reproductive phase[END_REF], displaying the red spot to a rival facing them head-on. Our results support a role for the red operculum spot in male contests, with the expression of a red spot primarily associated with sexually mature males and correlating with body size. However, this finding contradicts our prediction that parental males, which experience intrasexual contests and intersexual mate choice, should be more likely than cuckolder males to express a red operculum spot.

The expression of a red carotenoid-based colour signal may convey information about male condition, which could be functional in the context of male-male contests.

Carotenoids are acquired solely in the diet by vertebrates and these compounds have a number of critical physiological functions, making carotenoid-based signals a potentially reliable indication of individual condition and fighting ability [START_REF] Olson | Costly sexual signals: are carotenoids rare, risky or required?[END_REF].

Whether additional information is conveyed by the size, reflectance spectra and intensity of the red colour spot, or whether the red operculum spot comprises one component of a multimodal signal in pumpkinseed, remains to be investigated. Experimental studies would be particularly effective in differentiating these aspects of signal evolution in pumpkinseed.

Surprisingly, we detected no significant difference between parental and cuckolder males, either in the probability of expressing a red spot, or in the strength of the relationship between red spot size and male size. Alternative mating strategies in fish, which are exclusively shown by males in teleost fishes, are reflected most strikingly in reproductive behaviour, but also in external appearance and physiology. In contrast to [START_REF] Gross | Sneakers, satellites and parentals: polymorphic mating strategies in North American sunfishes[END_REF] interpretation of his sunfish model to a single native population of pumpkinseed, our analysis of non-native populations in Europe suggests that, while a clear distinction can be seen between the parental and cuckolder male reproductive strategies in relative testes size, the GSI value threshold that distinguishes parental and cuckolder males varies among populations, ranging from 2.3 (Poland) to 4.2 (England).

A possible alternative explanation for the observed range of GSI values may relate to variation in pumpkinseed population density, which has consequences for the reproductive success of cuckolder males. In water bodies where the populations occurred at high density, traits associated with sperm competition, including large testes size, are favoured, thereby potentially selecting for relatively larger gonads (high GSI) in both parental and cuckolder males, with a consequent upward shift in the threshold between them. In the three ponds in England for which quantitative density data were available [START_REF] Fox | Seasonal reproductive allocation, localscale variation and environmental influences on life history traits of introduced pumpkinseed (Lepomis gibbosus) in Southern England[END_REF], pumpkinseed was virtually the only species present in the two ponds where the species have been observed in high density and males in these ponds had higher mean and maximum male GSI values than observed in the low-density pond. In water bodies where populations occur at low densities, where the risk of sperm competition for parental males is low, relative testes size is predicted to be smaller and the GSI threshold between parental and cuckolder males lower [START_REF] Wootton | Reproductive Biology of Teleost Fishes[END_REF].

A caveat to these conclusions was the finding that many females also expressed red spots in their opercula, albeit with lower probability than males (Figure 2). The size of red operculum spots in females also conveyed less information about size than in males (Figure 3). If red operculum spots function primarily in male-male contests, then it is unclear why females would display this trait. For example, the red operculum spot may have multiple functions, potentially signalling male size and fighting ability, but also undergoing intersexual selection in the context of mate choice. If the red operculum spot is selected through intersexual selection, then it is possible that females express the trait for the same reason as males. Female ornamentation is surprisingly common in teleosts, and in the context of the pumpkinseed mating system, in which males are the sole providers of parental care, the evolution of male mate choice and female ornamentation is possible [START_REF] Wootton | Reproductive Biology of Teleost Fishes[END_REF].

The expression of red operculum spots in males is adaptive solely through intra-sexual selection, whereas in females the expression of red spots may be through a genetic correlated response. In threespine stickleback (Gasterosteus aculeatus), red throat colouration is known to function in intra-and intersexual selection on males and is a reliable indicator of body condition and parasite resistance [START_REF] Milinski | Female sticklebacks use male coloration in mate choice and hence avoid parasitized males[END_REF][START_REF] Barber | Indirect fitness consequences of mate choice in sticklebacks: offspring of brighter males grow slowly but resist parasitic infections[END_REF]. The red ornament, which has a strong genetic component, is also found in females and its expression in females is consistent with pleiotropy [START_REF] Yong | Genetic architecture of conspicuous red ornaments in female threespine stickleback[END_REF]. Pleiotropy may also mediate the expression of a red operculum spot in female pumpkinseed, and possibly also in cuckolder males, though this conclusion relies on there being a shared genetic basis to red spots in males and females and no adaptive value in expressing a red spot in females and cuckolders, and these assumptions remain to be tested.

The presence of two male GSI thresholds in pumpkinseed populations in Belgium, France and Poland provides evidence that, contrary to [START_REF] Gross | Sneakers, satellites and parentals: polymorphic mating strategies in North American sunfishes[END_REF], both the satellite and sneaker strategies may be present in some non-native European pumpkinseed populations (Figure 1). Alternatively, these findings may simply reflect the greater heterogeneity in age at maturity within some European populations [START_REF] Copp | Growth and life history traits of introduced pumpkinseed (Lepomis gibbosus) in Europe, and the relevance to its potential invasiveness[END_REF][START_REF] Fox | Old world versus new world: life-history alterations in a successful invader introduced across Europe[END_REF], which is a possible outcome of multiple sources of introduction or adaptation to conditions outside the natural range of the species. Behavioural studies, complemented with morphological analyses, are needed to clarify the existence of a satellite strategy in European pumpkinseed populations. Similarly, variation in predation among populations has the capacity to modify the selective landscape for red operculum spots, since sexual selection for visual signals can be limited through natural selection by predators in teleosts [START_REF] Endler | Natural selection on color pattern in Poecilia reticulata[END_REF][START_REF] Wootton | Reproductive Biology of Teleost Fishes[END_REF]. Although outside the scope of the present study, the effect of selective predation on pumpkinseed in response to the expression of the red operculum spot is not known and warrants further research.

In conclusion, the present study demonstrates that the presence of a red operculum spot in pumpkinseed is associated with sexual maturation, primarily in males. The size of red spot correlates with body size in parental and cuckolder males, although there was no difference in the presence or relative size of the red operculum spot between male mating strategies. The observed GSI distributions suggest the possible existence of the satellite mating strategy in some European populations, but alternatively this may be due to variations in pumpkinseed population density or to selective predation. The function, if any, of a red operculum spot in females is not clear but is consistent with pleiotropy.

FIGURE 1 .

 1 FIGURE 1. Observed gonado-somatic index (GSI): frequency plots for pumpkinseed among country of origin with Gaussian distribution curves and separation indices (SI), generated using Bhattacharyya's (1967) method. Threshold GSI values for separating frequency distributions for male mating strategies are indicated with arrows.

FIGURE 2 .FIGURE 3 .

 23 FIGURE 2. Mean fitted probability (solid line) of pumpkinseed expressing a red operculum spot as a function of standard length (SL in mm) with 95% confidence intervals (shaded area) for females, cuckolder males and parental males. Data were modelled with a Bernoulli GLMM. Fish from different populations were fitted as random intercepts in the model.

TABLE 1 .

 1 Country, population, location (WGS84 latitude and longitude coordinates) and capture method for male pumpkinseed Lepomis gibbosus, including number of specimens (n), mean and SD of standard length (mm) and the mean, minimum (min.) and maximum (max.) of the gonado-somatic index (GSI) for each population by water body.

	GSI

TABLE 2 .

 2 Mean (SD) pumpkinseed body mass, standard length (SL) and gonado-somatic 499 index (GSI; n/a = not available) among countries of origin and reproductive status. Male 500 reproductive strategy (either cuckolder or parental) was identified by separation of GSI 501 frequency distributions byBhattacharya's (1967) method. 502

	Country of origin	Reproductive status	n	Mass (g)	SL (mm)	GSI
	Belgium	Female	61	8.33	(9.89) 61.3 (16.0)	n/a
		Juvenile	10	3.85	(1.63) 52.5 (6.5)	n/a
		Cuckolder	21	4.81	(2.59) 55.9 (7.3)	3.60 (0.86)
		Parental	24	14.86 (10.22) 75.2 (16.4)	1.26 (0.59)
	France	Female	66	17.81 (11.68) 81.1 (18.5)	n/a
		Juvenile	16	1.28	(0.75) 38.9 (6.1)	n/a
		Cuckolder	15	5.03	(2.95) 55.6 (8.9)	4.90 (1.11)
		Parental	91	22.58 (14.21) 88.2 (17.1)	2.00 (0.83)
	Netherlands Female	62	23.77 (23.99) 86.7 (30.2)	n/a
		Juvenile	1	7.00	66.6	n/a
		Cuckolder	33	7.25	(1.16) 66.3 (4.5)	4.25 (1.32)
		Parental	51	36.35 (26.58) 101.3 (30.4)	2.10 (0.63)
	Poland	Female	18	20.25 (11.76) 76.4 (14.5)	n/a
		Juvenile	0	n/a	n/a	n/a
		Cuckolder	10	11.01	(7.24) 64.2 (10.3)	4.19 (1.05)
		Parental	39	23.79 (12.42) 80.4 (14.4)	1.26 (0.62)
	England	Female	180	14.46 (14.64) 76.5 (19.9)	n/a
		Juvenile	45	2.56	(3.14) 46.2 (13.4)	n/a
		Cuckolder	16	5.69	(7.08) 59.4 (14.4)	6.20 (2.23)
		Parental	165	12.29 (10.04) 76.1 (16.2)	n/a
	503

ACKNOWLEDGEMENTS

This study, which derives from work initiated as part of a N.A.T.O. Collaborative Linkage Grant (LST-CLG No. 979499), supported by consecutive research contracts (SF0238 and SF0248) from the UK Department of Environment, Food and Rural Affairs (to GHC), which were complemented by a Marie Curie post-doctoral fellowship (PIEF-GA-2008-219707) and a subsequent research grant from the National Science Centre, Poland (decision No DEC-2011/01/D/NZ8/01807) to GZ. We thank two anonymous reviewers for their constructive and perceptive comments.

Author contribution

The study was conceived by GHC and MGF. Specimens were collected and processed by GHC, GZ, MGF, SY, EZ, GM, JC, HV and HHvK. Data were analysed by GHC, SY, GZ, MP and CS. The initial draft of the manuscript was prepared by SY, EZ, GZ and GHC, then extensively revised by GZ, CS and MP, with co-author contributions.

2


Ivan Duchemin
email: ivan.duchemin@cea.fr

Xavier Blase

Cubic-scaling all-electron GW calculations with a separable density-fitting space-time approach

We present an implementation of the GW space-time approach that allows cubicscaling all-electron calculations with standard Gaussian basis sets without exploiting any localization nor sparsity considerations. The independent-electron susceptibility is constructed in a time representation over a non-uniform distribution of real-space locations {r k } optimized within a separable resolution-of-the-identity framework to reproduce standard Coulomb-fitting calculations with meV accuracy. The compactness of the obtained {r k } distribution leads to a crossover with the standard Coulomb-fitting scheme for system sizes below a few hundred electrons. The needed analytic continuation follows a recent approach that requires the continuation of the screened Coulomb potential rather than the much more structured self-energy. The present scheme is benchmarked over large molecular sets and scaling properties are demonstrated on a family of defected hexagonal boron-nitride flakes containing up to 6000 electrons.

Introduction

The GW approximation [START_REF] Hedin | New Method for Calculating the One-Particle Green's Function with Application to the Electron-Gas Problem[END_REF][START_REF] Strinati | Dynamical Correlation Effects on the Quasiparticle Bloch States of a Covalent Crystal[END_REF][START_REF] Hybertsen | Electron correlation in semiconductors and insulators: Band gaps and quasiparticle energies[END_REF][START_REF] Godby | Self-energy operators and exchangecorrelation potentials in semiconductors[END_REF][5][6][7][8][9][10][11] to the exchange-correlation self-energy has become a standard approach in solid-state physics to explore the electronic properties of metallic or semiconducting materials. Its accuracy was indeed proven superior to standard DFT calculations relying on the Kohn-Sham ansatz for the electronic energy levels (for large benchmark calculations on inorganic crystals, see e.g. Refs. 12,13). Further, following early applications in the late 90s, [14][START_REF] Van Der Horst | Ab Initio Calculation of the Electronic and Optical Excitations in Polythiophene: Effects of Intra-and Interchain Screening[END_REF][START_REF] Rohlfing | Optical Excitations in Conjugated Polymers[END_REF] the GW formalism is nowadays widely used as well for the study of gas phase or dense ordered or disordered organic molecular systems. The development of codes exploiting standard Gaussian atomic basis sets allowed in particular the comparison of all-electron GW calculations with higher level quantum-chemistry techniques (e.g. coupled-cluster) performed with the very same running parameters (geometry, atomic basis sets, resolution-of-the-identity, etc.) [START_REF] Bruneval | Benchmarking the Starting Points of the GW Approximation for Molecules[END_REF][START_REF] Krause | Coupled-cluster reference values for the GW27 and GW100 test sets for the assessment of GW methods[END_REF][START_REF] Rangel | Evaluating the GW Approximation with CCSD(T) for Charged Excitations Across the Oligoacenes[END_REF][START_REF] Knight | Accurate Ionization Potentials and Electron Affinities of Acceptor Molecules III: A Benchmark of GW Methods[END_REF] The scaling of the number of operations needed to perform GW calculations with respect to the system size is typically O(N 4 ) within traditional planewave implementations. This scaling can be preserved with localized basis sets provided that resolution-of-the-identity (RI) techniques [START_REF] Whitten | Coulombic potential energy integrals and approximations[END_REF][START_REF] Baerends | Self-consistent molecular Hartree-Fock-Slater calculations I. The computational procedure[END_REF][START_REF] Dunlap | On some approximations in applications of Xα theory[END_REF][START_REF] Vahtras | Integral approximations for LCAO-SCF calculations[END_REF][START_REF] Klopper | Explicitly correlated second-order Møller-Plesset methods with auxiliary basis sets[END_REF][START_REF] Ren | Resolution-of-identity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions[END_REF][START_REF] Duchemin | Hybrid and Constrained Resolution-of-Identity Techniques for Coulomb Integrals[END_REF] are used to avoid calculating response functions, such as the susceptibility, in the product space associated with valence-to-virtual molecular orbital products. While such a moderate scaling already allows calculations on systems containing well over a hundred atoms on supercomputers, [START_REF] Duchemin | Short-Range to Long-Range Charge-Transfer Excitations in the Zincbacteriochlorin-Bacteriochlorin Complex: A Bethe-Salpeter Study[END_REF][START_REF] Govoni | Large Scale GW Calculations[END_REF][START_REF] Li | Correlated electron-hole mechanism for molecular doping in organic semiconductors[END_REF][START_REF] Ben | Large-scale GW calculations on pre-exascale HPC systems[END_REF] attempts to deliver GW calculations with a lower scaling appeared with the seminal space-time approach by Rojas, Needs, Godby in 1995 [START_REF] Rojas | Space-Time Method for Ab Initio Calcu-lations of Self-Energies and Dielectric Response Functions of Solids[END_REF] and are now blooming. [START_REF] Foerster | An O(N 3 ) implementation of Hedin's GW approximation for molecules[END_REF][START_REF] Vlček | Stochastic GW Calculations for Molecules[END_REF][START_REF] Neuhauser | Expeditious Stochastic Calculation of Random-Phase Approximation Energies for Thousands of Electrons in Three Dimensions[END_REF][START_REF] Liu | Cubic scaling GW : Towards fast quasiparticle calculations[END_REF][START_REF] Wilhelm | Toward GW Calculations on Thousands of Atoms[END_REF][START_REF] Gao | Accelerating Time-Dependent Density Functional Theory and GW Calculations for Molecules and Nanoclusters with Symmetry Adapted Interpolative Separable Density Fitting[END_REF][START_REF] Kim | Complex-time shredded propagator method for large-scale GW calculations[END_REF][START_REF] Kutepov | Self-consistent GW method: O(N) algorithm for polarizability and self energy[END_REF][START_REF] Förster | Low-Order Scaling G 0 W 0 by Pair Atomic Density Fitting[END_REF][START_REF] Wilhelm | Low-Scaling GW with Benchmark Accuracy and Application to Phosphorene Nanosheets[END_REF] This space-time formalism [START_REF] Rojas | Space-Time Method for Ab Initio Calcu-lations of Self-Energies and Dielectric Response Functions of Solids[END_REF] stands as the first cubic-scaling GW approach relying on the separability of the independent-electron susceptibility χ 0 as the product of two Green's functions when expressed over a real-space grid, adopting further a time representation. This factorisation allows decoupling the summation over occupied and virtual molecular orbital contributions, leading to a cubic scaling scheme instead of the traditional quartic scaling calculation of χ 0 . Such a reduced scaling does not rely on any localization nor sparsity considerations associated with e.g. 3-center integrals within the local direct overlap metric [START_REF] Foerster | An O(N 3 ) implementation of Hedin's GW approximation for molecules[END_REF][START_REF] Wilhelm | Toward GW Calculations on Thousands of Atoms[END_REF][START_REF] Förster | Low-Order Scaling G 0 W 0 by Pair Atomic Density Fitting[END_REF] or a range-truncated Coulomb metric. [START_REF] Wilhelm | Low-Scaling GW with Benchmark Accuracy and Application to Phosphorene Nanosheets[END_REF] The imaginary-time formulation at the heart of the GW space-time approach is identical to the Laplace transform idea already in use in quantum chemistry for e.g. MP2 calculations. [START_REF] Almlöf | Elimination of energy denominators in Møller-Plesset perturbation theory by a Laplace transform approach[END_REF][START_REF] Häser | Laplace transform techniques in Møller-Plesset perturbation theory[END_REF] On the contrary, the use of a real space-grid was more naturally rooted in the pseudopotential planewave community, allowing by Fourier transform of the planewave basis to obtain relatively sparse uniform real-space grids. The space-time GW approach was more recently adapted to a full potential projector-augmented wave methodology, [START_REF] Liu | Cubic scaling GW : Towards fast quasiparticle calculations[END_REF] building on an earlier application to calculating RPA correlation energies with cubic scaling. [START_REF] Kaltak | Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in Si[END_REF] In the case of all-electron calculations, the size of the real-space grid may seem an a priori bottleneck. However, real-space quadrature strategies have been already developed with much success in quantum chemistry for accelerating the calculation of 2-electron Coulomb integrals, including the chain-of-sphere (COSX) semi-numerical approach to exchange integrals [START_REF] Neese | Efficient, approximate and parallel Hartree-Fock and hybrid DFT calculations. A 'chain-of-spheres' algorithm for the Hartree-Fock exchange[END_REF][START_REF] Izsák | An overlap fitted chain of spheres exchange method[END_REF][START_REF] Izsák | Robust fitting techniques in the chain of spheres approximation to the Fock exchange: The role of the complementary space[END_REF] or the general tensor hypercontraction mathematical framework in its specific least-square grid optimization implementation (LS-THC). [START_REF] Parrish | Tensor hypercontraction. II. Least-squares renormalization[END_REF][START_REF] Hohenstein | Communication: Tensor hypercontraction. III. Least-squares tensor hypercontraction for the determination of correlated wavefunctions[END_REF][START_REF] Schumacher | Tensor Hypercontraction Second-Order Møller-Plesset Perturbation Theory: Grid Optimization and Reaction Energies[END_REF] More recently, the interpolative separable density fitting (ISDF) approach [START_REF] Lu | Compression of the electron repulsion integral tensor in tensor hypercontraction format with cubic scaling cost[END_REF][START_REF] Lu | Fast algorithm for periodic density fitting for Bloch waves[END_REF] represents a versatile strategy to combine the standard quantum chemistry resolution-of-the-identity (RI) techniques with a separable representation of the coefficients of molecular orbital products over auxiliary basis sets. The ISDF approach is now developing in the pseudopotential planewave or real-space grid community, [START_REF] Gao | Accelerating Time-Dependent Density Functional Theory and GW Calculations for Molecules and Nanoclusters with Symmetry Adapted Interpolative Separable Density Fitting[END_REF][START_REF] Hu | Interpolative Separable Density Fitting Decomposition for Accelerating Hybrid Density Functional Calculations with Applications to Defects in Silicon[END_REF][START_REF] Dong | Interpolative Separable Density Fitting through Centroidal Voronoi Tessellation with Applications to Hybrid Functional Electronic Structure Calculations[END_REF][START_REF] Hu | Chapter Accelerating Optical Absorption Spectra and Exciton Energy Computation via Interpolative Separable Density Fitting[END_REF][START_REF] Hu | Accelerating Excitation Energy Computation in Molecules and Solids within Linear-Response Time-Dependent Density Functional Theory via Interpolative Separable Density Fitting Decomposition[END_REF] including a recent GW implementation. [START_REF] Gao | Accelerating Time-Dependent Density Functional Theory and GW Calculations for Molecules and Nanoclusters with Symmetry Adapted Interpolative Separable Density Fitting[END_REF] Similarly, building on the expertise with resolution-of-the-identity (RI) techniques and/or real-space quadratures for Coulomb integrals, the ISDF scheme is also being explored by the quantum chemistry community working with localized (e.g. Gaussian) basis sets for explicitly correlated all-electron calculations such as QMC or Möller-Plesset techniques. [START_REF] Malone | Overcoming the Memory Bottleneck in Auxiliary Field Quantum Monte Carlo Simulations with Interpolative Separable Density Fitting[END_REF][START_REF] Lee | Systematically Improvable Tensor Hypercontraction: Interpolative Separable Density-Fitting for Molecules Applied to Exact Exchange, Second-and Third-Order Møller-Plesset Perturbation Theory[END_REF] In a recent study, we presented an alternative to the ISDF formalism applied to allelectron Hartree-Fock, MP2 and RPA calculations with Gaussian basis sets. [START_REF] Duchemin | Separable resolution-of-the-identity with all-electron Gaussian bases: Application to cubic-scaling RPA[END_REF] In this scheme, standard auxiliary {P µ } basis sets (e.g. cc-pVXZ-RI [START_REF] Weigend | Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations[END_REF] or def2-XZVP-RI [START_REF] Weigend | RI-MP2: optimized auxiliary basis sets and demonstration of efficiency[END_REF] ) are provided as an input, as in any standard resolution-of-the-identity (RI) calculation, but the fitting procedure takes as an intermediate the expression of wave function and densities over compact non-uniform real-space grids {r k }. The corresponding fitting weights result from solving a quadrature equation that aims at reproducing the results of a standard Coulomb-fitting (RI-V) calculation. Adopting the space-time approach, cubic-scaling calculations of the independent-electron susceptibility at imaginary frequencies, and resulting RPA correlation energy, could be achieved. As compared to the corresponding RI-V calculation, an accuracy of a few µHartree/electron was demonstrated for RPA correlation energies with {r k } distributions typically 4 times larger than the input auxiliary basis, allowing a crossover with the standard quartic-scaling RI-V RPA calculations for systems of the size of pentacene.

In the present work, we extend this formalism to all-electron cubic-scaling GW calculations. In contrast with the RPA formalism, where only imaginary-frequency susceptibilities are needed, we further exploit a recently developed [START_REF] Friedrich | Tetrahedron integration method for strongly varying functions: Application to the GT self-energy[END_REF][START_REF] Duchemin | Robust Analytic-Continuation Approach to Many-Body GW Calculations[END_REF] analytic continuation scheme that brings to the real-frequency axis the dynamically screened Coulomb potential W rather than the much more structured GW self-energy. As compared to standard RI-V GW calculations, we demonstrate an accuracy at the meV level for the quasiparticle energies of large molecular sets. Finally, cubic scaling is evidenced using a family of defected hexagonal boron-nitride flakes with increasing radius containing up to 6000 electrons, with a crossover with the standard quartic scaling RI-V GW calculations for systems containing a very few hundred electrons.

Theory

The GW formalism with resolution-of-the-identity

We start by describing the standard resolution-of-the-identity (RI) framework for GW calculations. A more detailed discussion on RI techniques applied to MBPT can be found in Ref. 62. We just recall here that the essence of RI approximations, developed in particular to tackle the calculation of 2-electron 4-center Coulomb integrals with localized atomic orbital (AO) basis sets, [START_REF] Whitten | Coulombic potential energy integrals and approximations[END_REF][START_REF] Baerends | Self-consistent molecular Hartree-Fock-Slater calculations I. The computational procedure[END_REF][START_REF] Dunlap | On some approximations in applications of Xα theory[END_REF][START_REF] Vahtras | Integral approximations for LCAO-SCF calculations[END_REF][START_REF] Klopper | Explicitly correlated second-order Møller-Plesset methods with auxiliary basis sets[END_REF][START_REF] Ren | Resolution-of-identity approach to Hartree-Fock, hybrid density functionals, RPA, MP2 and GW with numeric atom-centered orbital basis functions[END_REF][START_REF] Duchemin | Hybrid and Constrained Resolution-of-Identity Techniques for Coulomb Integrals[END_REF] amounts to expressing the product of 2 molecular orbitals (MOs) over an auxiliary basis set {P µ }, namely:

φ n (r)φ m (r) = µ F µ (φ n φ m )P µ (r) (1)
where we work with finite size systems allowing real molecular orbitals. For localized atomicorbitals (AO) basis calculations, the auxiliary basis is typically 2-3 times larger than the AO basis set used to expand the MOs. As an example, the accurate RI-V Coulomb-fitting approach [START_REF] Vahtras | Integral approximations for LCAO-SCF calculations[END_REF] defines the coefficients F µ as:

F V µ (φ n φ m ) = ν [V -1 ] µν (P ν |φ n φ m ) (2)
with (P ν |φ n φ m ) the 3-center Coulomb integrals

(P ν |φ n φ m ) = drdr P ν (r)φ n (r )φ m (r ) |r -r |
and V µν the Coulomb matrix elements in the auxiliary basis. Coming now to the GW formalism, we start with the expression of the independent-electron susceptibility along the imaginary-frequency axis:

χ 0 (r, r ; iω) = 2 ja φ * j (r)φ a (r)φ * a (r )φ j (r ) iω -(ε a -ε i ) + c.c. (3)
where (j) and (a) index occupied and virtual molecular orbitals (MOs), respectively, and where the factor (2) indicates a closed shell system. Expanding MO products over an auxiliary basis in the case of real-valued MOs leads to:

χ 0 (r, r ; iω) = µν P µ (r) • [χ RI 0 (iω)] µν • P ν (r ) (4)
with

[χ RI 0 (iω)] µν = 2 ja F µ (φ j φ a )F ν (φ j φ a ) iω -(ε a -ε j ) + c.c. (5)
In the standard RI framework, it is the quantity [χ RI 0 (iω)] µν that is calculated in the auxiliary {P µ } basis. The following steps start with the definition of the GW correlation self-energy as a convoluted integral along the real-energy axis

Σ C (r, r ; E) = i 2π ∞ -∞ dω e iω0 + G(r, r ; E + ω) W (r, r ; ω) (6)
with W = (W -V ) and where G, W and V are the time-ordered 1-body Green's function, the screened and bare Coulomb potentials, respectively. In the contour-deformation approach, [START_REF] Godby | Self-energy operators and exchangecorrelation potentials in semiconductors[END_REF]5 this expression is transformed into an integral along the imaginary-energy axis, plus the contribution of a few residues involving the screened Coulomb potential W calculated at real energies:

Σ GW C (r, r ; E) = -1 2π ∞ -∞ dω G(r, r ; E + iω) W (r, r ; iω) (7)
- i φ i (r)φ i (r ) W (r, r ; ε i -E)θ(ε i -E) + a φ a (r)φ a (r ) W (r, r ; E -ε a )θ(E -ε a )
Expressing the Green's function in a quasiparticle form

G(r, r ; iω) = n φ n (r)φ n (r ) iω -ε n + iη × sgn(ε n -µ) (8)
with η = 0 + and µ the chemical potential, it appears that the expectation value of the GW correlation self-energy φ n |Σ GW C |φ n operator only requires integrals of the kind:

φ n φ m |W (z)|φ n φ m = µν F µ (φ m φ n )F ν (φ m φ n ) P µ |W (z)|P ν (9)
with z along the imaginary or real-energy axes. The needed P µ |W (z)|P µ matrix elements of the screened Coulomb potential are obtained from a Dyson-like equation projected into the auxiliary basis

P µ | W (z)|P ν = P µ |V |P ν + ζρ P µ |V |P ζ • [χ RI 0 (z)] ζρ • P ρ |W (z)|P ν ( 10
)
where the random phase approximation (RPA) approximation is used.

2.2

The space-time approach from a separable resolution-of-identity framework With the size of the auxiliary basis scaling linearly with system size, straightforward calculation of the [χ RI 0 (iω)] µν matrix elements from equation 5 requires O(N 4 ) steps. In other words, for each (P µ , P ν ) pair a double summation over occupied and unoccupied MOs is required. Following Almlöf and Häser, [START_REF] Almlöf | Elimination of energy denominators in Møller-Plesset perturbation theory by a Laplace transform approach[END_REF][START_REF] Häser | Laplace transform techniques in Møller-Plesset perturbation theory[END_REF] the Laplace transform:

1 iω -(ε a -ε j ) + c.c. = -2 +∞ 0 dτ cos(ωτ )e -(εa-ε j )τ (11)
where the time integral converges since (ε jε a ) < 0, allows to disentangle occupied and virtual energy levels in the denominator, leading to an imaginary time formulation

[χ RI 0 (iτ )] µν = -2i ja F µ (φ j φ a )F ν (φ j φ a )e ε j τ e -εaτ (12)
where the (i) factor is introduced to match the standard definition of the independentelectron susceptibility in the time domain. However, occupied and virtual MOs are still entangled in the F µ/ν (φ j φ a ) weight factors. This is precisely the goal of the separable RI introduced in Ref. 94 in the context of RPA total energies, with the expansion :

F RS µ (φ j φ a ) = k M µk φ j (r k )φ a (r k ) (13)
where the φ j and φ a MOs are factorized, leading to the wording separable-RI that we label RI-RS where RS stands for real-space. The present scheme targets a standard calculation with input molecular orbitals (MO) Gaussian basis and its associated auxiliary basis sets, and the {r k } distribution is an intermediate representation designed to reproduce the accuracy of a standard Coulomb-fitting calculation with the input basis sets. This is described here below and in Ref. 94 for Hartre-Fock, RPA and MP2 calculations.

The separable real-space RI (RI-RS) leads to expressing the independent-electron susceptibility matrix elements in the auxiliary basis as

[χ 0 (iτ )] µν RI-RS = kk M µk • χ 0 (r k , r k ; iτ ) • M νk (14)
where for any (r, r ) pair of points in real-space

χ 0 (r, r ; iτ ) = -iG(r, r ; iτ )G(r , r; -iτ ) (15)
G(r, r ; iτ ) = i occ j φ j (r)φ j (r )e ε j τ (τ > 0) (16) = -i vir a φ a (r)φ a (r )e εaτ (τ < 0) ( 17
)
where G is the time-ordered one-body Green's function. Considering here that the MOs are spin-orbitals, one can recover the factor two in Equation 12 for spin-restricted systems. Eq. 15

is the standard equation for the original space-time approach [START_REF] Rojas | Space-Time Method for Ab Initio Calcu-lations of Self-Energies and Dielectric Response Functions of Solids[END_REF] with the summations over occupied and virtual MOs completely decoupled, leading to a strictly cubic-scaling number of operations, independently of any localization nor sparsity properties. After construction of the imaginary-time independent-electron susceptibility in the {P µ } Gaussian auxiliary basis following Eq. 14, the [χ RI 0 (iω)] µν are obtained by Fourier transform at imaginary frequencies.

Following Eq. 10, the screened Coulomb potential P µ | W (z)|P ν is finally obtained along the imaginary frequency axis. Details about the time and frequency grids will be given below when discussing the analytic continuation of W to the real-axis. The overall flow of calculations is presented in Fig. 1. While we use in the present scheme the analytic continuation of W , the standard analytic continuation of Σ to the real axis can also be used once the χ 0 (P, Q; iω) susceptibilities are obtained using the real-space imaginary-time approach.

Construction of the {r k } distributions

A crucial aspect of the present scheme is the size of the {r k }-set that controls the prefactor associated with the present cubic-scaling scheme. Following our implementation of cubicscaling RPA calculations in an all-electron space-time approach, [START_REF] Duchemin | Separable resolution-of-the-identity with all-electron Gaussian bases: Application to cubic-scaling RPA[END_REF] the central idea is not to use a generic real-space grid (such as Becke grids 99 adopted to express densities in DFT codes) but to optimize for each atomic species a reduced set of {r k } points sufficient to reach the accuracy of the standard Coulomb-fitting RI-V approximation in conjunction with the chosen auxiliary basis set. Such task is performed by minimizing the difference between the F RS and F V fitting procedure as defined in Eqs.

F RS µ (αα ) -F V µ (αα ) P µ 2 V ( 18
)
where the {α} are the Gaussian basis functions used to expand the MOs.

For a given atom, the initial set of points are constructed as a superposition of high symmetry subsets of Lebedev grids up to order 9, associated with different sphere radii (see Supporting Information Ref. 94). The optimization process starts by minimizing the penalty function of Eq. 18 adjusting first the radii. This is similar to the optimization strategy adopted in the grid-based formulation of LS-THC 85 but fitting the codensities coefficients rather than the 4-center Coulomb integrals. In a second step, all constraints are raised and every point is allowed to move independently. This non-linear minimization process is performed using a basin-hopping mechanism coupled to a L-BFGS (limited memory Broyden-Fletcher-Goldfarb-Shanno) algorithm. We emphasize that such a step is done once for all for a given element and the chosen basis sets. Experimenting with such a strategy for the def2-TZVP / def2-TZVP-RI associated basis sets leads to 100 {r k } points for H and He, and 336, 436 and 536 points for elements in the second, third, and fourth row of the periodic table, respectively. Such grid sizes allow an agreement at the meV level between subsequent quasiparticle energies calculated with the present real-space approach and the standard Coulomb-fitting RI-V scheme. Except for the first row, this is typically 3.5 to 4.5 times larger than the number of elements in the def2-TZVP-RI set. Better optimization schemes and reducing the initial number of Lebedev subsets may lead to reduce these {r k } distribution sizes. However, as shown below, the present approach already provides an excellent accuracy-to-cost ratio and appears to be very robust, with no outliers as tested over large molecular sets.

In a second step, the {r k } distribution for the molecular system is built as the superposition of the isolated atoms {r k } distributions and only the weights {M µk }, as defined in Eq. 13, need to be calculated for each considered molecular system. Such a step only requires O(N 3 ) operations since the least-square estimator matrix [M ] µk is obtained in a matrix multiplication/inversion formulation from the target F V µ coefficients (see Ref. 94).

As such, the weights {M µk } are univocally defined once the {r k } grid and F V µ factors are set-up. We provide in the Supporting Information a graph confirming the cubic-scaling of the {M µk } construction that amounts to about 25% of the total CPU time, including the calculation of the target F V µ , for non-self-consistent G 0 W 0 calculations. This 2-step process, namely the optimization of the {r k }-distribution on isolated atoms, dramatically simplifies the minimization process while preserving excellent accuracy as demonstrated below.

Concerning previous RPA calculations using the present separable RI with Laplace transform scheme, [START_REF] Duchemin | Separable resolution-of-the-identity with all-electron Gaussian bases: Application to cubic-scaling RPA[END_REF] the crucial observation was that indeed {r k }-sets typically 4 times larger than the used auxiliary {P µ } basis set were sufficient to reproduce the accuracy of standard RI-V calculations to within a few µHartree/electron for the exchange, RPA and MP2 total energies. Namely, replacing the standard RI-V F V µ coefficients by their F RS µ approximants preserved an excellent accuracy, with a number of points sufficiently small to offer a crossover with the standard quartic-scaling RI-V RPA approach for systems of the size of pentacene.

We will perform here below the corresponding accuracy check for the quantity of interest here, namely the GW quasiparticle energies, showing that meV accuracy can be achieved with a crossover between the present separable RI-RS scheme and standard RI-V calculations for systems containing less than a few hundred electrons. This crossover is independent of the compactness and dimensionality of the studied systems since sparsity and localization are not exploited.

Analytic continuation, frequency and time grids

An important aspect of the space-time approach is the required analytic continuation from the imaginary to the real-frequency axis. With the calculation of the susceptibility χ 0 (iτ ) at imaginary times, the imaginary-frequency χ 0 (iω) analog can be obtained by Fourier transform. From such quantities, the screened Coulomb potential W (iω) and self-energy Σ(iω) can be obtained at imaginary frequencies and efficiently continued analytically to the real energy axis as performed in many codes.

An alternative to the analytic continuation of the self-energy was proposed by Christoph

Friedrich in the context of GT calculations on solid iron, [START_REF] Friedrich | Tetrahedron integration method for strongly varying functions: Application to the GT self-energy[END_REF] and by ourselves in the present case of GW calculations on molecular systems with extensive benchmark accuracy checks. [START_REF] Duchemin | Robust Analytic-Continuation Approach to Many-Body GW Calculations[END_REF] The central idea is to adopt the contour deformation scheme where the quantity needed along the real-axis is no longer the self-energy directly, but the screened Coulomb potential (see second and third lines of Eq. 7). Since the self-energy contains (N W × N G ) poles, where N G and N W are respectively the number of poles of the Green's function and screened Coulomb potential, the screened Coulomb potential is much less structured than the self-energy itself, leading to a much more robust analytic-continuation scheme. Difficult test cases drawn from the GW 100 test sets, [START_REF] Van Setten | Benchmarking G 0 W 0 for Molecular Systems[END_REF] such as the M gO or BN dimers, were shown to be very accurately treated with the calculations of the screened Coulomb potential W (iω) for no more than 12 frequencies along the imaginary-axis. In particular, theses frequencies are constructed so as to minimize the error over the imaginary axis integration contribution to Equation 7. We address the reader to Ref. 98 for a detailed presentation of this scheme. We keep the number of imaginary frequencies to nω=12 that was shown in this former study to lead to sub-meV accuracy, as compared to the contour deformation scheme, for GW calculations on frontier orbitals using this "robust" analytic continuation scheme. Once the imaginary frequencies are set, the corresponding imaginary-time grid is constructed following Ref. 94, where the present space-time approach was explored for RPA total energy calculations. The selected times {τ p , p = 1, nτ } are optimized for the chosen set of imaginary frequencies z k (k=1, nω) through the minimization process: arg min

ω p k ,τp   k ln(Emax) ln(E min ) du p ω p k e -τpe u - 1 e u + iz k - 1 e u -iz k 2   ( 19
)
where ω p k is the weight associated with a given τ p time for a targeted z k frequency. The energies E min and E max are the energy gap and the maximum (ε aε i ) value, respectively. The 1/(e u ±iz k ) factors represent the pole structure of the independent-electron susceptibility along the imaginary axis. The p ω p k e -τpe u approximant translates the fact that e -a|τ | (a > 0) is the Fourier transform of 2a/(a 2 +ω 2 ) within a prefactor. Following Refs. 94,98, the log scale is used so as to allow a regular sampling of the error oscillations at energies between E min and E max . The problem can be then solved in a traditional least square approach using a uniform sampling in u. Such a formulation conserves an excellent accuracy, as demonstrated in Table 1, with a number of grid points comparable to that commonly adopted with the more elaborated minimax approach. [START_REF] Liu | Cubic scaling GW : Towards fast quasiparticle calculations[END_REF][START_REF] Kaltak | Cubic scaling algorithm for the random phase approximation: Self-interstitials and vacancies in Si[END_REF]100 Similarly to the minimax formulation, the grids points τ p have been pre-tabulated with the E max /E min ratio as a single parameter, so as to minimize the computational effort of the setup. On the other hand, the ω p k coefficients can be conveniently recalculated on the fly as the result of a simple linear least square equation. In association with nω=12 imaginary frequencies, nτ =18 times are selected to reach sub-meV accuracy on the quasiparticle energies. We provide in Table 1 a typical test of accuracy, selecting the def2-TZVP G 0 W 0 @PBE0 HOMO and LUMO energies of acridine, the first element of the molecular set of Ref. 41 studied in full details in the next section. In this Table, RI-RS without Laplace transform only differs from the standard Coulomb-fitting (RI-V) by the construction of the F µ fitting coefficients. We observe in particular that the dependence on the nω grid falls well below the meV for nω ≥ 10. The real-space approach with Laplace-transform (RI-RS+LT) depends further on the imaginary-time grid. However, for a given nω imaginary-frequency grid, a time-grid with nτ = 1.5×nω introduces negligible errors, comforting overall our choice of nω = 12 and nτ =18 running parameters.

Results

Validation and accuracy

We benchmark the accuracy of the present scheme using the recent set of 24 intermediate size molecules with acceptor character proposed in Ref. 41. Our calculations are performed at the def2-TZVP G 0 W 0 @PBE0 level associated with the corresponding def2-TZVP-RI auxiliary basis. [START_REF] Weigend | RI-MP2: optimized auxiliary basis sets and demonstration of efficiency[END_REF] Our goal here is not to carry calculations in the complete-basis set limit, but rather to assess the accuracy of the present space-time approach, as compared to the standard Coulomb-fitting (RI-V) scheme, using a reasonable basis set. The real-space {r k } sets were thus optimized for the def2-TZVP and def2-TZVP-RI basis sets association, following the scheme described above and summarized in Eqn. 18. As discussed above, the size of the {r k } atomic distributions amounts to 136 for H and 336 for second row elements.

We provide in Table 2 the def2-TZVP G 0 W 0 @PBE0 highest occupied (HOMO) and lowest unoccupied (LUMO) molecular orbitals energies calculated using the standard Coulombfitting (RI-V) scheme and the separable real-space (RI-RS) approach in conjunction with the Laplace transform scheme. Both calculations are performed with the "robust" analytic continuation (AC) scheme. [START_REF] Duchemin | Robust Analytic-Continuation Approach to Many-Body GW Calculations[END_REF] Direct contour-deformation calculations without any analytic continuation, namely calculating directly the needed residues of W along the real-axis within the standard quartic scaling RI-V formalism, shows that the errors introduced by the AC are well below the meV for the HOMO and LUMO energy levels of the molecules contained in this set. The analysis of the results evidences that for such {r K } distributions, the error on the quasiparticle energies remains below the meV. Such an accuracy may be tuned by increasing/decreasing the size of the {r k } distributions, but the present accuracy-to-size trade-off is already excellent in practice.

We further perform benchmark def2-TZVP G 0 W 0 @PBE0 calculations on the GW 100 test set [START_REF] Krause | Coupled-cluster reference values for the GW27 and GW100 test sets for the assessment of GW methods[END_REF][START_REF] Van Setten | Benchmarking G 0 W 0 for Molecular Systems[END_REF][START_REF] Maggio | A Plane Wave Perspective for Small Molecules[END_REF][START_REF] Govoni | Comparison of Methods and Accuracy of Results Obtained with the WEST Code[END_REF][START_REF] Gao | Accelerating Time-Dependent Density Functional Theory and GW Calculations for Molecules and Nanoclusters with Symmetry Adapted Interpolative Separable Density Fitting[END_REF][START_REF] Duchemin | Robust Analytic-Continuation Approach to Many-Body GW Calculations[END_REF]101 that contains elements from the third and fourth periods of the periodic table, including transition metal complexes. We exclude the 5 systems containing 5-th period elements for which the def2-TZVP basis set requires the use of an effective core potential, namely Xe, Rb 2 , I 2 , vinyl iodide (C 2 H 3 I) and aluminum iodide (AlI 3 ). Similarly to the previous test set, the error induced by the real-space RI-RS with Laplace transform technique, as compared to a standard RI-V calculations, is below the meV for most molecules as reported in Fig. 2. Three systems (Kr, CuCN, SF 4 ) show an error on the HOMO slightly larger than 1 meV (in absolute value), while all errors on the LUMO value are below the Table 2: HOMO and LUMO energies at the def2-TZVP G 0 W 0 @PBE0 level for the molecular set of Ref. 41. Molecules are arranged by chemical families in the order of Ref. 41. The standard Coulomb-fitting scheme (RI-V) performed with the def2-TZVP-RI auxiliary basis set [START_REF] Weigend | RI-MP2: optimized auxiliary basis sets and demonstration of efficiency[END_REF] is compared to the present real-space Laplace-transform (RI-RS) calculations. Negative and positive maximum errors, the mean absolute (MAE) and mean signed (MSE) 2: HOMO/LUMO def2-TZVP G 0 W 0 @PBE quasi-particle energy discrepancy analysis over the GW 100 test, excluding the 5 systems containing 5-th period elements (see text).

The error is that of the present real-space Laplace-transform (RI-RS) approach with respect to the standard Coulomb-fitting (RI-V) approach performed with the corresponding def2-TZVP-RI auxiliary basis set. [START_REF] Weigend | RI-MP2: optimized auxiliary basis sets and demonstration of efficiency[END_REF] The molecules of the set are sorted according to the maximum period (row) involved within the periodic table.

meV. As emphasized above, optimizing further the distribution of {r k } points may bring all errors below the meV, but the purpose of the present study is to show that the present scheme, as it stands, already brings very consistently the error at the meV level. Overall, the mean absolute (MAE) errors amount to 0.21 meV and 0.09 meV for the HOMO/LUMO, respectively. All data can be found in the Supporting Information. The present data confirm, in the specific case of GW calculations, previous studies reporting on the excellent accuracy-to-cost ratio associated with grid-based techniques for the evaluation of exact exchange or explicit correlation energies in the context of all-electron atomic-orbital basis sets calculations. [START_REF] Neese | Efficient, approximate and parallel Hartree-Fock and hybrid DFT calculations. A 'chain-of-spheres' algorithm for the Hartree-Fock exchange[END_REF][START_REF] Izsák | An overlap fitted chain of spheres exchange method[END_REF][START_REF] Izsák | Robust fitting techniques in the chain of spheres approximation to the Fock exchange: The role of the complementary space[END_REF][START_REF] Parrish | Tensor hypercontraction. II. Least-squares renormalization[END_REF][START_REF] Hohenstein | Communication: Tensor hypercontraction. III. Least-squares tensor hypercontraction for the determination of correlated wavefunctions[END_REF][START_REF] Schumacher | Tensor Hypercontraction Second-Order Møller-Plesset Perturbation Theory: Grid Optimization and Reaction Energies[END_REF][START_REF] Duchemin | Separable resolution-of-the-identity with all-electron Gaussian bases: Application to cubic-scaling RPA[END_REF]102,103

Scaling analysis

We finally address the issue of scaling with respect to the system size through the example of finite-size hexagonal boron-nitride (h-BN) "flakes" containing a central point-defect. The study of the optical emission mediated by defects in h-BN is an important technological research area, with the prospect of having at hand stable, room-temperature, polarized and ultrabright single-photon sources, together with a scientific challenge when it comes to identify the defects and mechanisms responsible for such sharp emission lines in the visible range. [104][105][106][107][108] We select as a test case the C B V N (nitrogen vacancy plus carbon substitution to neighbouring boron) defect that has been recently identified as a possible candidate for emission at about 2 eV. 109 Our goal here is not to confirm the likeliness of such a defect, but rather to start exploring whether many-body calculations with a typical defect can be performed using finite-size clusters, rather than the traditional supercell approach using periodic boundary conditions (PBC). Indeed, the use of PBC complicates the calculations of charged excitations due to the electrostatic interaction between cells, and the Coulomb potential must be truncated to avoid spurious contributions even in the limit of large supercells. As such, the modeling of the opto-electronic properties of defects in h-BN at the many-body GW and Bethe-Salpeter level remains scarce due in particular to the cost of performing the required large-scale GW calculations. 109,110 The edge of the h-BN flakes are passivated by hydrogen atoms to avoid dangling bonds and the HOMO-LUMO gap is clearly controlled by very localized defect states yielding energy levels within the gap of pristine h-BN (see Inset Fig. 3 for the LUMO). The size of the studied flakes correspond to average radii ranging from 21.5 to 56.4 Å, containing from 137 to 941 C, B or N atoms, that is from 167 to 1019 atoms including passivating H atoms.

The average diameter is defined as D = 2 N at S at /π, where N at is the number of B or N atoms, and S at = 3 √ 3d 2 BN /4 is the effective surface per B or N atom in the hexagonal lattice, with d BN the BN bond length. Structural relaxation at the PBE0 6-311G* level indicates that the ground-state for these systems is not spin-polarized.

Before discussing scaling properties, we briefly comment on the evolution of the HOMO-LUMO energy gap obtained at the 6-311G* G 0 W 0 @PBE0 level 111 as a function of system Figure 3: Plot of the 6-311G* G 0 W 0 @PBE0 HOMO-LUMO gap as a function of the inverse flake average diameter. The Inset represents the Kohn-Sham LUMO associated with our largest flake (1019 atoms). size (see Fig. 3). Our goal here is not to obtain converged values with respect to basis completeness, but to study the evolution of the gap with system size using a minimal triplezeta plus polarization basis, keeping in mind that the 6-311G* HOMO-LUMO gaps for our defected flakes are typically 70 meV larger than that obtained with the larger def2-TZVP basis set. The decrease of the gap with increasing diameter can be attributed to polarization effects, namely the fact that upon calculating the ionization potential or electronic affinity, as measured by a photo-emission experiment, the added charge localized on the defect generates a long-range Coulomb field that polarizes the surrounding atoms. Such a polarization, properly described within the GW formalism, stabilizes the added hole or electron, closing the gap. In the case of finite size systems, this polarization is incomplete as compared to an infinite sheet, leading to a HOMO-LUMO gap that is too large. Performing a fit of the GW gap up to second order in (1/D), the linear contribution is found to be negligible, leading to a simple quadratic dependence. Such a quadratic behaviour stems from the reaction field generated by the 2D density of dipoles induced by a charge added or removed on/from the LUMO/HOMO levels. 112 The extrapolated gap at infinite radius amounts to 4.96 eV, still 60 meV away from the gap of the largest flake studied. As expected, the G 0 W 0 gap is much larger than the PBE0 Kohn-Sham gap of 3.07 eV obtained for the largest system. Fitting the data associated with the 4 smallest flakes, the extrapolated value remains within 20 meV of the extrapolated value with the fit containing all points. This indicates the stability of the extrapolation scheme and suggests that an accurate asymptotic value may be obtained with calculations performed on systems containing a rather limited number of atoms. While the isolated defect limit GW quasiparticle gap needs extrapolating to infinite sizes, preliminary results indicate that optical excitations, which are neutral excitations of the system, converge much faster as a function of system size.

We now plot (log scale) in Fig. 4 the total CPU time (namely the sum of all cores CPU time, over the complete run) associated with the G 0 W 0 calculations reported in Fig. 3. We compare in particular the standard RI-V (Coulomb fitting) calculations (filled blue triangles) performed with the universal Coulomb fitting auxiliary basis of Ref. 113 and the present realspace Laplace-transform approach (filled black circles) with a real-space {r k } distribution optimized as described above for this specific auxiliary basis set. Timings are reported in Tables S2 andS3 of the Supporting Information where the specific contributions from the RI-RS set-up, susceptibility calculations, Dyson-equation inversion and self-energy calculations are further provided.

These calculations confirm that the present space-time scheme offers a cubic scaling with system size (see black dot-dashed fit) with an (extrapolated) crossover with the standard quartic scaling RI-V scheme taking place for about 350 electrons. We further explore limit of small system sizes with a larger def2-TZVP basis associated with its def2-TZVP-RI auxiliary basis. We can see that despite the larger basis sets, the overhead of the 128 core distribution and pre-computation phases still hinders slightly the fit by perfect cubic/quartic lines in the small system size limit. Nonetheless, these later calculations confirms a crossover that consistently takes place at about 350 electrons ( 50-60 B/N atoms). A similar crossover was observed in the case of cc-pVTZ RPA calculations where we used a similar real-space Figure 4: Scaling properties (log scale) for the standard Coulomb-fitting (RI-V) G 0 W 0 calculations compared with the present real-space Laplace-transform (RI-RS+LT) scheme. Calculations have been performed with the def2-TZVP and 6-311G* basis sets. For the standard RI-V scheme, the auxiliary def2-TZVP-RI [START_REF] Weigend | RI-MP2: optimized auxiliary basis sets and demonstration of efficiency[END_REF] and universal Coulomb fitting 113 basis sets, respectively, were used. Calculations were performed on a set of hexagonal boron-nitride flakes with up to 6000 electrons. Calculations have been performed on AMD Rome Epyc nodes with 128 cores/node and 1.85 Gb/core. The black (RI-RS+LT) and blue (RI-V) dot-dashed lines are cubic and quartic fits, respectively. An unconstrained fit yields a scaling exponent of 3.07 for the (RI-RS+LT) scheme.

Laplace-transform approach to build the independent-electron susceptibility (see Ref. 94).

All calculations have been performed on a supercomputer built of 128 cores 2,6 GHz AMD Rome nodes with 1.85 Gb/core memory. We only used fully filled nodes in order to maintain consistency between the timings presented here, meaning that the smallest system calculations have been distributed on a minimum of 128 cores. Under this constraint, we selected CPU grid sizes that roughly match the minimum memory requirement for each calculation, as detailed in Table S2 of the Supporting Information. Our real-space Laplacetransform (RI-RS+LT) calculations require much less cores than the standard RI-V scheme, partly due to the corresponding O(N 2 ) memory footprint, with the χ 0 (r k , r k ; iω) being the largest objects stored in memory. On the other hand, the memory requirement of our standard RI-V implementation grows as O(N 3 ), dominated in this case by the storage of (occupied)×(virtual) co-density auxiliary fits F V µ (φ i φ a ). Let us emphasize that within the RI-RS+LT approach, each 3-center integral is computed and immediately discarded during the RI setup.

Conclusions

We have presented an all-electron space-time GW formalism relying on a separable resolutionof-the-identity (RI) formalism, offering cubic-scaling GW calculations that do not exploit any sparsity nor localization considerations. This allows a crossover with the quartic scaling Coulomb-fitting RI-V GW calculations for systems containing a very few hundred electrons, independently of the dimensionality of the studied system. As compared to the interpolative separable density fitting (ISDF) scheme, the present approach preserves the use of standard auxiliary Gaussian basis sets that are taken as an input, and not constructed by the ISDF algorithm. The needed distribution of {r k } points are optimized to recover at the meV level the results of a standard Coulomb-fitting (RI-V) GW calculation. Precalculated grids to be associated with a larger collection of standard basis sets, beyond the 6-311G* and def2-TZVP sets adopted in this study, comes now as a prerequisite for a broader use of the present scheme. Scaling with system size could be further reduced possibly by exploiting stochastic techniques [START_REF] Neuhauser | Expeditious Stochastic Calculation of Random-Phase Approximation Energies for Thousands of Electrons in Three Dimensions[END_REF] or the decay properties of the space-time Green's functions 114 at long-range in the case of very large systems. The performances as they stand today clearly illustrate however the interest of real-space quadrature, as developed by the quantum chemistry community for all-electron atomic-orbital basis sets calculations, and the progress performed by the GW community to install this family of many-body perturbation techniques as a valuable tool for moderately correlated systems, with an excellent trade-off between accuracy and CPU time, allowing to tackle finite size or periodic systems, metallic or semiconducting, containing up to several hundred atoms.

(

Graphical TOC Entry

Scaling with system size for 6-311G* G 0 W 0 calculations performed with the cubic scaling realspace Laplace transform scheme (RI-RS+LT) as compared to the quartic scaling traditional Coulomb fitting (RI-V) approach for a family of defected boron-nitride flakes. We further give the total CPU time associated with the main cubic-scaling operations (Fig. S1 and Table S3), namely the construction of the grid weights M µ,k (RI setup), the construction of the independent-electron susceptibilities (χ 0 ), the inversion of the Dyson equation to obtain the screened Coulomb potential (W ), and finally the evaluation of the self-energy (Σ GW ). Taking the largest system, for which the 6-311G* space-time G 0 W 0 calculation represents a total CPU time of about 14336 hours, the sum of these 4 steps amounts to 13984 hours, namely 98% of the total. For this largest flake, the RI set-up represents about 25% and the χ 0 calculations about 32% of the total CPU time.

Table S3: Number of atoms, electrons, cores and timings in hh:mm (total CPU time) for the principal computational steps of the 6-311G* RI-RS+LT G 0 W 0 @PBE0 calculations.

Figure 1 :

 1 Figure 1: Schematic representation of the steps involved in the present cubic scaling allelectron space-time approach. The optimized set of real-space positions {r k } is typically 4 times as large as the corresponding input Gaussian auxiliary basis {P µ }. The number of imaginary times and frequencies is set by nτ and nω. In the contour deformation approach (see Inset) only the screened Coulomb potential W needs to be continued from the imaginary to the real-energy axis, avoiding the continuation of the much more structured self-energy (see Ref. 98).

Figure S1 :

 S1 Figure S1: Scaling behavior for the 4 principal steps of our cubic-scaling G 0 W 0 implementation. Dashed lines are guides for the eyes using an exact cubic behavior (slope equal to 3).

Table 1 :

 1 Acridine def2-TZVP G 0 W 0 @PBE0 HOMO and LUMO (in eV) convergence against the imaginary frequency grid size. We keep the ratio nτ = 1.5 × nω. The real-space with Laplace transform (RI-RS + LT) depends on both frequency and time grids. All calculations are performed with the auxiliary def2-TZVP-RI basis set, while RI-RS calculations use an extra {r k } distribution optimized for the corresponding def2-TZVP/def2-TZVP-RI basis sets association.

	nω	RI-V	RI-RS RI-RS + LT
			HOMO
	6 -7.55777 -7.55777	-7.55787
	8 -7.56073 -7.56072	-7.56071
	10 -7.56113 -7.56112	-7.56112
	12 -7.56112 -7.56111	-7.56111
	14 -7.56112 -7.56111	--
			LUMO
	6 -0.75723 -0.75719	-0.75713
	8 -0.76020 -0.76016	-0.76018
	10 -0.76046 -0.76042	-0.76042
	12 -0.76044 -0.76040	-0.76040
	14 -0.76044 -0.76040	--

  errors are indicated in meV. Values leading to the largest error are in bold.

		HOMO	LUMO
		RI-V [eV]	RI-RS [eV] RI-V [eV]	RI-RS [eV]
	anthracene	-7.0787	-7.0787	-0.4233	-0.4234
	acridine	-7.5611	-7.5611	-0.7604	-0.7604
	phenazine Figure	-7.9755	-7.9754	-1.1796	-1.1799
	azulene	-7.1271	-7.1272	-0.5595	-0.5597
	benzoquinone (BQ)	-9.7275	-9.7273	-1.6004	-1.6006
	naphthalenedione	-9.3297	-9.3296	-1.5520	-1.5521
	dichlone	-9.4112	-9.4111	-1.9812	-1.9814
	F4-BQ	-10.5855	-10.5853	-2.3221	-2.3220
	Cl4-BQ	-9.7245	-9.7246	-2.5192	-2.5191
	nitrobenzene	-9.7868	-9.7867	-0.5485	-0.5486
	F4-benzenedicarbonitrile	-10.2336	-10.2334	-1.7207	-1.7210
	dinitro-benzonitrile	-10.7596	-10.7595	-1.8683	-1.8684
	nitro-benzonitrile	-10.2038	-10.2038	-1.3992	-1.3993
	benzonitrile	-9.5108	-9.5107	0.1831	0.1828
	fumaronitrile	-10.9712	-10.9710	-1.0740	-1.0740
	mDCNB	-10.0173	-10.0175	-0.7227	-0.7228
	TCNE	-11.4676	-11.4674	-3.2543	-3.2544
	TCNQ	-9.1373	-9.1373	-3.5795	-3.5795
	maleic-anhydride	-10.7783	-10.7782	-1.0029	-1.0031
	phthalimide	-9.6757	-9.6755	-0.6380	-0.6381
	phthalic-anhydride	-10.1111	-10.1111	-0.9060	-0.9061
	Cl4-isobenzofuranedione	-9.5809	-9.5807	-1.7124	-1.7127
	NDCA	-8.7061	-8.7058	-1.3302	-1.3304
	BODIPY	-7.8008	-7.8008	-1.6315	-1.6316
	Max. err.	+0.32 meV / -0.23 meV +0.07 meV / -0.33 meV
	MAE	0.12 meV	0.14 meV
	MSE	0.08 meV	0.12 meV

  5) Farid, B.; Daling, R.; Lenstra, D.; van Haeringen, W. GW approach to the calculation of electron self-energies in semiconductors. Phys.Rev. B 1988, 38, 7530-7534.

	(6) Aryasetiawan, F.; Gunnarsson, O. The GW method. Rep. Prog. Phys. 1998, 61, 237-
	312.
	(7) Farid, B. In Electron Correlation in the Solid State -Chapter 3 ; March, N., Ed.;
	Imperial College Press, London, 1999.
	(8) Onida, G.; Reining, L.; Rubio, A. Electronic excitations: density-functional versus
	many-body Green's-function approaches. Rev. Mod. Phys. 2002, 74, 601-659.
	(9) Ping, Y.; Rocca, D.; Galli, G. Electronic excitations in light absorbers for photo-
	electrochemical energy conversion: first principles calculations based on many body
	perturbation theory. Chem. Soc. Rev. 2013, 42, 2437-2469.
	(10) Martin, R.; Reining, L.; Ceperley, D. Interacting Electrons: Theory and Computational
	Approaches; Cambridge University Press, 2016.
	Phys. Rev. Lett. 2006, 96, 226402.
	(13) Shishkin, M.; Marsman, M.; Kresse, G. Accurate Quasiparticle Spectra from Self-
	Consistent GW Calculations with Vertex Corrections. Phys. Rev. Lett. 2007, 99,
	246403.
	(14) Ethridge, E. C.; Fry, J. L.; Zaider, M. Quasiparticle spectra of trans-polyacetylene.

(11) Golze, D.; Dvorak, M.; Rinke, P. The GW Compendium: A Practical Guide to Theoretical Photoemission Spectroscopy. Front. Chem. 2019, 7, 377. (12) van Schilfgaarde, M.; Kotani, T.; Faleev, S. Quasiparticle Self-Consistent GW Theory. Phys. Rev. B 1996, 53, 3662-3668.

Table S1 -

 S1 GW100 G 0 W 0 @PBE0 data, continued from previous page

		HOMO	LUMO
		RI-V [eV]	RI-RS [eV] RI-V [eV]	RI-RS [eV]
	Arsine	-10.1110	-10.1111	3.0944	3.0946
	Benzene	-9.0143	-9.0140	1.5996	1.5997
	Berylliummonoxide	-12.9594	-12.9595	-11.6442	-11.6442
	Borane	-13.0064	-13.0064	0.4548	0.4547
	Boronnitride	-11.3781	-11.3782	-3.5945	-3.5945
	Bromine	-10.1504	-10.1495	-0.8253	-0.8247
	Buthane	-11.5753	-11.5754	2.9304	2.9302
	Carbondioxide	-13.4373	-13.4371	2.9403	2.9403
	Carbondisulfide	-9.7561	-9.7562	0.2424	0.2424
	Carbonmonoxide	-13.9578	-13.9576	1.0780	1.0779
	Carbonoxyselenide	-10.1634	-10.1633	1.3249	1.3248
	Carbonoxysulfide	-10.9454	-10.9453	1.7287	1.7287
	Carbontetrabromide	-9.9729	-9.9719	-0.4853	-0.4847
	Carbontetrachloride	-11.1414	-11.1409	0.7781	0.7780
	Carbontetrafluoride	-15.7426	-15.7419	4.9648	4.9653
	Chlorine	-11.1315	-11.1316	-0.1991	-0.1991
	Coppercyanide	-10.0069	-10.0052	-1.0328	-1.0327
	Copperdimer	-7.3478	-7.3479	-0.3626	-0.3624
	Cyclooctatetraene	-8.0888	-8.0886	0.5793	0.5791
	Cyclopentadiene	-8.3697	-8.3699	1.5707	1.5707
	Cyclopropane	-10.6135	-10.6138	3.5310	3.5311
	Cytosine	-8.4546	-8.4546	0.8030	0.8029
	Diborane6	-11.9952	-11.9951	1.2461	1.2460
	Dipotassium	-3.9849	-3.9849	-0.4377	-0.4377
	Disilane	-10.3480	-10.3480	2.2701	2.2702
	Ethane	-12.4466	-12.4468	3.1998	3.1998
	Ethanol	-10.3587	-10.3592	2.9515	2.9515
	Ethoxyethane	-9.5363	-9.5364	3.0017	3.0014
	Ethybenzene	-8.5805	-8.5805	1.5685	1.5687
	Ethylene	-10.3456	-10.3455	2.5831	2.5831
	Fluorine	-15.3320	-15.3322	0.1373	0.1372
	Fluoroborane	-10.7671	-10.7671	1.5093	1.5092
	Formaldehyde	-10.5488	-10.5487	1.5744	1.5745
	Formic Acid	-11.0648	-11.0647	2.5916	2.5915
	Galliummonochloride	-9.5438	-9.5438	0.3533	0.3532
	Germane	-12.1362	-12.1362	3.2973	3.2973
	Guanine	-7.7447	-7.7449	1.3749	1.3750
	Helium	-23.7926	-23.7926	22.2760	22.2760
	Hexafluorobenzene	-9.6562	-9.6563	0.9102	0.9101
	Hydrazene	-9.4265	-9.4264	2.6692	2.6691
				Continued on next page

Table S1 -

 S1 GW100 G 0 W 0 @PBE0 data, continued from previous page

		HOMO	LUMO
		RI-V [eV]	RI-RS [eV] RI-V [eV]	RI-RS [eV]
	Hydrogen	-15.9491	-15.9491	4.4635	4.4635
	Hydrogenazide	-10.4807	-10.4811	1.8777	1.8776
	Hydrogenchloride	-12.2631	-12.2632	2.8900	2.8900
	Hydrogencyanide	-13.3566	-13.3569	3.1549	3.1550
	Hydrogenfluoride	-15.5491	-15.5491	3.2673	3.2673
	Hydrogenperoxide	-11.2224	-11.2222	3.1763	3.1764
	Hydrogensulfide	-10.0084	-10.0084	3.2376	3.2376
	Krypton	-13.4876	-13.4865	10.4118	10.4122
	Lithiumdimer	-5.0389	-5.0389	-0.1772	-0.1771
	Lithiumfluoride	-10.4604	-10.4603	0.1233	0.1232
	Lithiumhydride	-7.3073	-7.3072	0.1841	0.1841
	Magnesiumchloride	-11.2469	-11.2469	-0.1141	-0.1142
	Magnesiumfluoride	-12.9920	-12.9921	0.0598	0.0599
	Magnesiummonoxide	-7.3227	-7.3226	-1.5708	-1.5708
	Methane	-14.0449	-14.0452	3.5594	3.5594
	Methanol	-10.7266	-10.7268	3.1640	3.1639
	Neon	-20.7658	-20.7657	20.8885	20.8885
	Nitrogen	-15.2455	-15.2455	2.9042	2.9041
	Ozon	-12.3624	-12.3620	-1.8652	-1.8654
	Pentasilane	-8.9797	-8.9796	0.7289	0.7289
	Phenol	-8.3968	-8.3966	1.4663	1.4663
	Phenol (v2)	-8.2916	-8.2912	1.5195	1.5196
	Phosphine	-10.2794	-10.2793	3.1459	3.1460
	Phosphorusdimer	-10.1723	-10.1723	-0.2715	-0.2715
	Phosphorusmononitride	-11.4976	-11.4975	0.2857	0.2857
	Potassiumbromide	-7.6723	-7.6718	-0.2585	-0.2584
	Potassiumhydride	-5.4022	-5.4022	0.0491	0.0492
	Propane	-11.8751	-11.8750	3.0092	3.0093
	Pyridine	-9.3883	-9.3877	1.0332	1.0330
	Silane	-12.4656	-12.4657	3.2217	3.2217
	Sodiumchloride	-8.5723	-8.5725	-0.4087	-0.4087
	Sodiumdimer	-4.9076	-4.9076	-0.3663	-0.3663
	Sodiumhexamer	-4.2910	-4.2910	-0.6794	-0.6794
	Sodiumtetramer	-4.1678	-4.1677	-0.7052	-0.7052
	Sulferdioxide	-12.0048	-12.0045	-0.4411	-0.4412
	Sulfertetrafluoride	-12.2741	-12.2753	0.9815	0.9813
	Tetracarbon	-11.0137	-11.0139	-2.5154	-2.5154
	Thymine	-8.8343	-8.8341	0.5947	0.5946
	Titaniumfluoride	-14.6675	-14.6670	-0.2178	-0.2178
	Tuloene	-8.6321	-8.6316	1.5334	1.5333
				Continued on next page

Table S1 -

 S1 GW100 G 0 W 0 @PBE0 data, continued from previous page Timings for the 6-311G* G 0 W 0 @PBE0 calculationsWe present in TableS2additional details on the size, timing (walltime) and number of cores used concerning the defected h-BN systems. Calculations have been performed on 128-cores 2.6 GHz AMD Rome nodes with 1.85 Gb/core. For consistency between different systems, all calculations used fully filled nodes. Notice that real-space Laplace-transform (RI-RS+LT) calculations required much less cores than the standard RI-V scheme as the signature of reduced memory requirements.

		HOMO	LUMO
		RI-V [eV]	RI-RS [eV] RI-V [eV]	RI-RS [eV]
	Uracil	-9.2341	-9.2340	0.5442	0.5440
	Urea	-9.6647	-9.6650	2.4063	2.4063
	Vynilbromide	-8.9477	-8.9475	1.9401	1.9401
	Vynilbromidev2	-9.4866	-9.4861	1.8184	1.8183
	Vynilchloride	-9.7883	-9.7883	2.0109	2.0108
	Vynilfluoride	-10.2372	-10.2377	2.7246	2.7247
	Water	-12.1646	-12.1646	3.0779	3.0779
	Max. err.	+1.67 meV / -1.16 meV +0.66 meV / -0.30 meV
	MAE	0.21 meV	0.09 meV
	MSE	0.08 meV	0.001 meV
	2

Table S2 :

 S2 Number of B,C,N and H atoms, number of electrons, number of cores and timings in hh:mm:ss (walltime) for the 6-311G* G 0 W 0 @PBE0 calculations on defected hexagonal boron-nitride flakes. The memory indicated represent the maximum memory available for the run. B/C/N N H N el N CP U walltime mem.(Tb) N CP U walltime mem.(Tb)

				RI-V			RI-RS+LT
	N 137	30 852	384	00:15:37	0.72	128	00:19:26	0.24
	203	36 1254	768	00:32:21	1.44	256	00:31:21	0.48
	371	48 2274 3072 01:21:12	5.76	640	01:08:41	1.20
	448	54 2742 6144 01:37:25	11.52	1152 01:19:18	2.16
	674	65 4108 23040 02:46:32	43.20	2560 01:59:11	4.80
	941	78 5724	-	-	-	5120 02:48:00	9.60

  CPU time N B/C/N N H N el N CP U RI setup

					χ 0	W	Σ GW
	137	30 852	128	11:37	10:58	2:09	14:52
	203	36 1254	256	33:56	39:40	6:45	49:31
	371	48 2274	640	171:42	202:50	39:37	301:34
	448	54 2742 1152	338:16	466:55	88:04	592:49
	674	65 4108 2560	1219:20 1555:57 241:13 1866:17
	941	78 5724 5120	3633:19 4564:49 715:42 5070:11

Acknowledgement

The authors are indebted to Pascal Pochet for discussions concerning defects stability in periodic and finite size systems and to Gabriele D'Avino for explaining the quadratic behaviour of the polarization energy with inverse diameter. Calculations have been performed thanks to an allocation on the French CEA-TGCC Joliot-Curie supercomputer comprizing 2292 AMD Rome Epyc nodes with 128 cores/node and 2 Go/core. This work received support from the French Agence Nationale de la Recherche (ANR) under contract ANR-20-CE29-0005.

Supporting Information Available

We provide in the Supporting Information the details of the GW 100 test set data (Table S1)

and additional information about the systems size, timings and number of cores used for the 6-311G* G 0 W 0 @PBE0 calculations on the hexagonal boron-nitride flakes (Tables S2-S3 and Fig. S1).

Supporting Information for:

Cubic-scaling all-electron GW calculations with a separable density-fitting space-time approach Ivan Duchemin * , † and Xavier Blase ‡ †Univ. Grenoble Alpes, CEA, IRIG-MEM-L Sim, 38054 Grenoble, France ‡Univ. Grenoble Alpes, CNRS, Inst NEEL, F-38042 Grenoble, France E-mail: ivan.duchemin@cea.fr 1 Details of the def2-TZVP G 0 W 0 calculations on the GW 100 test set.

Table S1: HOMO and LUMO energies at the def2-TZVP G 0 W 0 @PBE0 level for the GW 100 molecular test set. S1 The standard Coulomb-fitting scheme (RI-V) performed with the def2-TZVP-RI auxiliary basis set S2 is compared to the present real-space Laplace-transform (RI-RS) calculations. Negative and positive maximum errors, the mean absolute (MAE) and mean signed (MSE) errors are indicated in meV. Values leading to the largest error are in bold.