Eric Heitz

Stephen Hill

Morgan Mcguire

Combining Analytic Direct Illumination and Stochastic Shadows

Introduction

This research is motivated by the lack of shadowing in certain stateof-the-art real-time shading techniques, such as the area-lighting framework developed by [START_REF] Heitz | Realtime polygonal-light shading with linearly transformed cosines[END_REF]. While it accurately models the effects of complex materials and lighting, the lack of shadowing ultimately produces unrealistic images (Fig. 2(a)) and this problem cannot be solved robustly by using a shadow map or a fake soft shadow (Fig. 2(b,c)). Given the recent progress in real-time path-tracing [START_REF] Schied | Spatiotemporal varianceguided filtering: Real-time reconstruction for path-traced global illumination[END_REF][START_REF] Mara | An efficient denoising algorithm for global illumination[END_REF], it is tempting to believe that current real-time techniques might soon be replaced by fully stochastic approaches combined with efficient denoising algorithms. However, these stochastic techniques either suffer from noise or overblurred shading details due to the aggressive denoising that needs to be used for reconstruction (Fig. 2(d)). In summary, on the one hand we have real-time shading techniques that produce clean and sharp shading details but miss shadowing effects, and on the other we have stochastic techniques that incorporate all of the effects but yield noisy or overblurred results.

In this paper, we show how to get the best of both worlds by completing real-time shading techniques with ray-traced shadows in a way that is mathematically correct. Our main insight is the formulation of shadowed illumination as the product of the unshadowed illumination and the illumination-weighted shadow. Based on this idea, we develop the following contributions:

• A ratio estimator of shadowed illumination, in which shadowing is the only source of variance. The estimator produces noise-free results in unshadowed regions of the image.

• A denoising formulation that allows for denoising shadows separately from shading. Thanks to this, we can denoise our results without compromising the sharpness of shading details as only the shadows are overblurred (Fig. 2(e)). Note that any existing denoiser can benefit from this formulation.

• A real-time implementation of our estimator and a custom realtime denoiser renders high-quality results using only 1-4 shadow rays per light per pixel.

Mathematical Overview of our Method

Shadowed direct illumination The shadowed direct illumination at a point is the integral of the cosine-weighted Bidirectional Reflectance Distribution Function (BRDF), the incident lighting, and the visibility towards the light:

S = Ω BRDF × Light × Visibility. (1)
Real-time unshadowed direct illumination Many real-time techniques provide efficient ways for computing the direct illumination with various surface materials and light types but without accounting for visibility, i.e. they compute the unshadowed direct illumination:

U = Ω BRDF × Light. ( 2
)
For instance, [START_REF] Heitz | Realtime polygonal-light shading with linearly transformed cosines[END_REF] provide a closed-form evaluation for integrating the BRDF against polygonal area lights. Another typical example is the precomputation of the integral of the BRDF and an environment map using a prefiltered representation, which is a common practice in games [START_REF] Kautz | Unified approach to prefiltered environment maps[END_REF][START_REF] Karis | Real shading in unreal engine 4[END_REF]]. These techniques produce shadowless results, as shown in Figure 2(a).

Problem Adding mathematically correct shadows on top of an analytic computation of U is complicated because in general the visibility is not separable:

S = Ω BRDF × Light × Ω Visibility. (3)
Assuming separability of the visibility is mathematically wrong and might produce unrealistic results, as in Figure 2(c). Mathematically correct visibility is thus a stumbling block for real-time rendering and there is currently no alternative for obtaining exact and robust shadows other than turning to fully stochastic methods.

Our solution Our insight is to split the shadowed illumination into two parts:

S = Ω BRDF × Light unshadowed illumination U analytic or precomputed × Ω BRDF × Light × Visibility Ω BRDF × Light illumination-weighted shadow W stochastic + denoising . (4)
We refer to the first part as the unshadowed illumination U , which we compute via existing analytic or precomputed real-time techniques as in Equation ( 2). The second part is the illuminationweighted shadow W , which we compute via stochastic raytracing followed by denoising. This formulation is inspired by the ratio estimator, a variance-reduction technique from the field of statistics [START_REF] Hartley | Unbiased ratio estimators[END_REF][START_REF] Mickey | Some finite population unbiased ratio and regression estimators[END_REF].

Advantages This formulation provides two main advantages. First, only the illumination-weighted shadow W is computed stochastically. Consequently, only shadowed regions of the final image appear noisy and lit regions are rendered noise-free even at 1 sample per pixel. Second, the unshadowed illumination U containing all of the sharp shading details (textures, normals, etc.) is analytic and does not need to be denoised. Hence, we denoise only the stochastic illumination-weighted shadow W such that agressive denoising overblurs only shadows while sharp shading details are preserved in the final image.

Previous Work

Unshadowed Direct Illumination

Our technique is built on top of existing real-time techniques that compute direct illumination without shadows in an efficient manner, i.e. they provide an efficient and noise-free solution to Equation (2) and can be used in our formulation of the illumination shown in Equation (4).

Analytic illumination techniques originated from the study of form factors that yield the irradiance arriving at a point -i.e. direct illumination with a diffuse BRDF -from various geometric shapes [START_REF] Baum | Improving radiosity solutions through the use of analytically determined form-factors[END_REF]]. The advantage of diffuse BRDFs is that they do not depend on the view direction, making it possible to derive the irradiance from complex emitters such as linearly varying area lights [START_REF] Chen | A closed-form solution for the irradiance due to linearly-varying luminaires[END_REF]]. The concept of using analytic form factors has been extended to non-diffuse BRDFs with area lights [START_REF] Arvo | Applications of irradiance tensors to the simulation of non-lambertian phenomena[END_REF][START_REF] Snyder | Area light sources for real-time graphics[END_REF]. For a long time these computations were too costly, so real-time rendering applications used coarser approximations [START_REF] Drobot | Physically based area lights[END_REF]], but they were recently made affordable thanks to accurate approximations [START_REF] Lecocq | Accurate analytic approximations for real-time specular area lighting[END_REF]].

The derivation of new directional distributions that offer analytic integration properties for area lighting is an active research topic [START_REF] Heitz | Realtime polygonal-light shading with linearly transformed cosines[END_REF][START_REF] Heitz | Linear-light shading with linearly transformed cosines[END_REF][START_REF] Dupuy | A spherical cap preserving parameterization for spherical distributions[END_REF].

Precomputed illumination techniques are used for image-based emitters such as scene lighting represented as an environment map. Tyically, diffuse BRDFs use precomputed irradiance maps [Ramamoorthi and Hanrahan 2001] and glossy BRDF models use preconvolved environment maps [START_REF] Kautz | Unified approach to prefiltered environment maps[END_REF][START_REF] Ramamoorthi | Frequency space environment map rendering[END_REF][START_REF] Karis | Real shading in unreal engine 4[END_REF]]. The same approach can also be used for refraction [START_REF] De Rousiers | Real-time rough refraction[END_REF]].

Variance-Reduction Techniques

Ratio estimators Ratio estimators are variance-reduction techniques often used in statistics [START_REF] Hartley | Unbiased ratio estimators[END_REF][START_REF] Mickey | Some finite population unbiased ratio and regression estimators[END_REF]]. The goal is to estimate the expectation of a random variable Y by taking advantage of the information provided by a positively correlated random variable X whose expectation is known. A classic estimator consists in averaging observations yn from Y :

E [Y ] = lim N →∞ 1 N N n=1 yn. (5)
If X and Y are positively correlated and if the expectation E [X] of X is known, it is possible to reduce the variance using:

E [Y ] = E [X] lim N →∞ N n=1 yn N n=1 xn . ( 6
)
We obtained Equation (4) by applying this ratio estimator using the unshadowed direct illumination for X and the shadowed direct illumination for Y .

A ratio estimator has been employed previously in rendering by [START_REF] Stachowiak | Stochastic screen-space reflections[END_REF], who uses a precomputed integral of the BRDF alone in order to absorb the variance due to imperfect importance sampling of the BRDF. To our knowledge, the idea of using ratio estimators when the integral of one or several terms of the integrand is known has not been used further in rendering applications.

Control variates method

The main concurrent approach that takes advantage of analytic integrals is the control variates method, which has been used in rendering [START_REF] Hery | Physically based lighting at pixar[END_REF]. We compare against in Section 3. Control variates are a variance reduction technique that consists of using an easily integrable approximation of the function to integrate, together with an error-correction term that is computed stochastically:

S = c Ω BRDF×Light control variate (analytic) - Ω BRDF×Light×(c -Visibility) error (stochastic) , (7)
where c is a correlation coefficient that needs to be chosen carefully to ensure effective variance reduction. Computing the optimal value for c requires either supplemental computations or dedicated optimizations and data structures [START_REF] Fan | Optimizing control variate estimators for rendering[END_REF][START_REF] Clarberg | Exploiting visibility correlation in direct illumination[END_REF]. In this paper, we are interested in estimators that work as a replacement for a fully stochastic evaluation without precomputations or supplemental data structures. Hence, the optimal value for c is unknown, so we use c = 1 by default for our comparisons. In Section 3, we show that a control variate used with a default value does not perform well. The main advantage of our formulation in Equation ( 4) is that it does not require finding such a coefficient and hence it does not require supplemental computations nor data structures to be effective. Furthermore, we show that our formulation allows us to make robust separate approximations or denoising for the illumination and the shadows, which is not possible with control variates.

Ratio Estimator

In this section, we introduce our formulation of the direct illumination, build a ratio estimator based on this formulation, and discuss its properties and performance in terms of variance reduction. For brevity we will use a more compact notation: R denotes the cosineweighted BRDF, L the incident lighting, and V the visibility, and the shadowed and unshadowed direct illumination are

S = Ω R(ω) L(ω) V (ω) dω, (8)
U = Ω R(ω) L(ω) dω. (9)

Previous Estimators of the Shadowed Illumination

Full-stochastic estimator A classic estimator of S is

S sto N = 1 N N n=1 R(ωn) L(ωn) V (ωn) p(ωn) , (10)
where the directions ωn are importance sampled from a Probability Density Function (PDF) p. We call it the full-stochastic estimator, since all of the terms are evaluated randomly.

Control-variate estimator Following Equation ( 7), with a correlation coefficient c = 1 as explained in the previous work section, we obtain the control-variate estimator

S CV N = U - 1 N N n=1 R(ωn) L(ωn) (1 -V (ωn)) p(ωn) , (11)
which uses an analytic expression for the unshadowed illumination U and a stochastic estimator for the additive error term.

Split Formulation of the Shadowed Illumination

Following Equation ( 4), we formulate direct illumination as

S = U × W, (12)
W = S U , ( 13
)
where W is the illumination-weighted shadow. Intuitively, it is the average visibility of the rays weighted by their contribution to the illumination integral. This is shown more clearly if we write

W = Ω R(ω) L(ω) V (ω) dω Ω R(ω) L(ω) dω = Ω R(ω) L(ω) Ω R(ω) L(ω) dω V (ω) dω, (14)
which emphasizes that RL/ RL acts as a normalized filter over the binary visibility values V . As a result, the values of the weighted-shadow image are always distributed in [0, 1] and the more a ray contributes to the illumination -i.e. it has both high BRDF and light values -the more it is accounted for in the weighted-shadow image.

Interpretation of the illumination-weighted shadow

The value W can be seen as the probability that a ray sampled perfectly from the normalized product of the BRDF and the light is unshadowed.

Ratio Estimator of the Shadowed Illumination

The formulation of Equation ( 12) yields our ratio estimator:

S ratio N = U × S sto N U sto N , (15)
S sto N = 1 N N n=1 R(ωn) L(ωn) V (ωn) p(ωn) , U sto N = 1 N N n=1 R(ωn) L(ωn) p(ωn) ,
where S sto N and U sto N are stochastic estimators of the shadowed and unshadowed illumination that are computed with the same samples ωn generated from PDF p, i.e. they use the same random numbers. It is important that the R(ωn) L(ωn)/p(ωn) values are the same in S sto N and U sto N to ensure that the result of their division is a normalized average of visibility values V (ωn) that remains in [0, 1].

Note that any other sampling strategy can be used to evaluate S sto N and U sto N as long as they use the same samples. In our implementation, we use Multiple Importance Sampling (MIS) to compute all of the estimators with lower variance, i.e. we generate samples from both the BRDF and the light and we weight them with the balance heuristic [START_REF] Veach | Optimally combining sampling techniques for monte carlo rendering[END_REF].

The ratio estimator is biased The expectation of the ratio estimator is not the shadowed illumination. Indeed, for a fixed number of samples N , the expectation and the fraction do not commute:

E S sto N U sto N = S U ⇒ E S ratio N = S, (16)
and the analytic U does not cancel out with the expectation of U sto N such that only S would remain.

The ratio estimator is consistent The ratio estimator converges toward the right result as the number of samples increases because the limit and the fraction commute:

lim N →∞ S sto N U sto N = S U ⇒ lim N →∞ S ratio N = U × S U = S, (17)
and the analytic U effectively cancels out with the limit of U sto N such that only S remains.

full-stochastic

CV ratio (ours) reference

Figure 3: Variance of the estimators at 4spp. Note that the controlvariate (CV) estimator produces negative values that are clamped to 0 (black pixels). Because of this, the shadow appears gray but its average is effectively 0. Furthermore, in the last row, the controlvariate result appears less noisy, but the black pixels actually span more variance in the negative domain.

The ratio estimator has low variance In our experience, the variance of our ratio estimator is generally lower or similar to the variance of both the full-stochastic and control-variate estimators. However, we could not find a mathematical proof of this. To support this observation, we investigate the problem through numerical experiments. A comparison of the noisy results produced by the different estimators is available in Figures 3 and6, and our supplemental material provides a numerical convergence analysis with an equivalent 1D problem. Our results show that the full-stochastic estimator is generally good in shadowed regions but noisy in lit regions, while the control-variate estimator has the opposite property. Our ratio estimator is as good as the both estimators in their respective good configurations: i.e. it yields noise-free results in lit regions like the control-variate estimator and does not increase the variance in shadowed regions. In penumbra regions, the three estimators are roughly equivalent.

Approximate Ratio Estimators

Grayscale approximation of illumination-weighted shadows

We must compute the estimator for each RGB channel to be exact. This means storing three RGB images: analytic, stochastic shadowed, and stochastic unshadowed. To save memory and band-width, especially in denoising, it is possible to compute and process stochastic shadowed luminance and stochastic unshadowed luminance as before, but only a grayscale weighted-shadow image (Figure 4). The approximation diverges from the exact result if the scene exhibits strong correlation between directions and lightsource colors, as in Figure 4. However, the approximate result remains plausible and noise-free, and for many applications superior to using no shadow at all.

no shadow gray shadow (approx.) RGB shadow (exact)

Figure 4: Grayscale approximation of the illumination-weighted shadows. In order to save memory and bandwidth it is possible to compute a grayscale (rather than RGB) weighted-shadow image.

Approximate real-time unshadowed illumination So far, we assumed that the exact unshadowed illumination was available analytically or through a precomputed representation. However, in many cases, only an approximate result is available in real-time.

Our method allows us to extend even approximate real-time shading techniques with robust shadows, as shown in Figure 5. For instance, the integral of the BRDF and environmental lighting is hard to precompute for complicated BRDFs. The common real-time approximation for this is to use an environment map prefiltered by an approximate isotropic filter [START_REF] Karis | Real shading in unreal engine 4[END_REF]]. The approximation introduced by this technique is thus the replacement of the exact BRDF R by an approximate BRDF R:

Ω R L = Ω R L, (18)
which is what we desire to obtain a sharp and clean real-time result. However, for the stochastic computation of the weighted shadow, it is simpler to use the exact material R. Indeed, we usually know how to importance sample the exact R but we might not know how to sample the approximate R that has been used to prefilter the environment map, especially if the prefiltering has been further approximated and optimized [START_REF] Manson | Fast filtering of reflection probes[END_REF]. In this case, the denominator of the weighted shadow no longer cancels with the unshadowed illumination in our formulation:

Ω R L × Ω R L V Ω R L = Ω R L V. ( 19
)
Fortunately, this approximation is visually plausible and artifactfree. Indeed, the weighted shadow Ω R L V / Ω R L is a value in [0, 1] that modulates the approximate unshadowed illumination Ω R L. What we obtain in Figure 5 is thus the approximate unshadowed illumination modulated such that the exact shadow appears. This result is plausible if the real-time unshadowed illumination approximation is plausible.

Denoising the Ratio Estimator

We now explain how to apply our ratio estimator with a generic denoise(). In Section 4, we propose a novel denoiser for shadows.

unshadowed illumination shadowed illumination real-time approx. reference ratio approx. reference [START_REF] Karis | Real shading in unreal engine 4[END_REF] :

Ω R L Ω R L Ω R L × Ω R L V Ω R L Ω R L V
Figure 5: Approximate real-time unshadowed illumination. We use our formulation to incorporate the exact shadow with a realtime approximation of the unshadowed illumination.

Computations We compute three images: the analytic unshadowed illumination image U , the stochastic shadowed illumination SN , and the stochastic unshadowed illumination UN :

U = Ω R(ω) L(ω) dω, (20)
SN = N n=1 R(ωn) L(ωn) V (ωn) p(ωn) , (21)
UN = N n=1 R(ωn) L(ωn) p(ωn) , (22)
The analytic image U can be computed using an exact solution if it is available or with an approximate real-time technique otherwise, as explained in Section 3.4. We obtain the stochastic illuminationweighted shadow WN with the division

WN = SN UN , (23)
and our stochastic (noisy) result is the product

result = U × WN . ( 24
)
Denoising before dividing An important point is that one should denoise the stochastic shadowed and unshadowed illuminations before computing their division:

denoised result = U × denoise[SN ] denoise[UN ] , (25)
instead of denoising WN directly.

Mathematical justification An alternative could be to denoise the ratio instead, i.e. performing denoising after the division. To get an intuition, let us consider denoising a set of pixels that have exactly the same BRDF, light, and visibility, and should therefore evaluate to the same result. In this case, the denoising operation between these pixels should not introduce bias and we expect convergence towards the exact result as the number of denoised pixels increases. Similar to Equation ( 16), denoising after the division does not make the denoising operation converge toward the right result because the estimator is biased:

denoise SN UN -→ S U . ( 26
)
Similar to Equation ( 17), denoising before and dividing after makes the denoising operation converge toward the right result because the estimator is consistent:

denoise[SN ] denoise[UN ] -→ S U . ( 27
)
Practical justification Besides being the right mathematical approach, denoising before dividing (as in Equation ( 27)) also presents the advantage of being robust to pixels with low denominators or undefined values. Let us consider the case where the illumination weight evaluates to 0, i.e.

SN = UN = 0 =⇒ WN = SN UN = 0 0 = NaN. ( 28
)
If we denoised the weighted-shadow values WN = SN /UN , this NaN value would propagate and cause visual artifacts. In contrast, by denoising SN and UN before dividing as we suggest, they get mixed with their neighboring pixel values and become positive (as long as they are positive for at least one pixel in the neighborhood) and the division is well-defined. In general, with this approach, pixels with small denominators UN do not contribute to themselves or their neighborhoods, which is the expected behavior. Indeed, the weighted-shadow is by definition a weighted sum of shadow values.

The smaller the weights, the less relevant the shadow values.

Robustness for denoising Figures 6 and 7 apply a naive bilateral filter over the different estimators. This simple denoising model demonstrates the robustness of the estimators against aggressive filtering; we later show a more sophisticated denoiser. Naive denoising works well for all the estimators where neighboring pixels expect similar results. However, if the image contains highfrequency shading details, the methods behave differently. Denoising the full-stochastic estimator blurs these details. Denoising the control-variate estimator introduces visual artifacts because some pixels might have negative values of high magnitude. These negative values can be spread safely to similar neighboring pixels but they should not be spread to different pixels that do not expect so much variance in the negative domain. The consequence in this example is that denoising the control-variate estimator makes the shadow disappear. In contrast, denoising our weighted-shadow image, as explained above, only makes the shadow appear blurrier and shading details are preserved.

Using better denoisers

Recent denoisers, such as that of [START_REF] Schied | Spatiotemporal varianceguided filtering: Real-time reconstruction for path-traced global illumination[END_REF], produce significantly less overblurring of shading details and achieve better results compared to our simple bilateral filter in Figures 6 and7. We did not choose such a state-of-the-art denoiser for our comparisons because we wanted to validate the robustness of our formulation against rough denoising techniques, i.e. in the worst case.

Furthermore, the objective of recent research focusing on denoising is to produce the best possible denoise [SN ] as a final result. Our denoising formulation U × denoise[S N ] denoise[U N ] works on top of that. It uses the information provided by the analytic U to prevent shading noise and shading overblur. Thanks to this, our formulation ensures that, for any given denoise() function, our result is always better. Hence, even the most efficient direct-illumination denoisers can only be improved by our formulation. Furthermore, our formulation with an existing denoiser comes at negligible cost because the analytic computation of U is usually cheap and the denoising of UN is exactly the same as SN , so most computations can be factorized.

In summary, we directly benefit from existing (and future) efficient direct-illumination denoisers and, reciprocally, our formulation provides an improved way to use them that is simple to implement and comes at negligible overhead.

Our Real-Time Algorithm

In this section, we describe our real-time implementation dedicated to our primary motivating case of efficiently shadowing area lights.

U UN SN E denoise[S N ] denoise[U N ] U × denoise[S N ] denoise[U N ]

Main algorithm

Our main algorithm operates in five full-screen passes per frame, as illustrated in Figure 8:

Temporally integrate the final image via reverse reprojection

Filter passes 3 and 4 use an underlying Gaussian kernel with the bilateral weights driven by plane differences between the G-buffer position and normal at the center of the kernel and at the tap location. This is a standard denoising practice.

We set the standard deviation of the Gaussian kernel proportional to the noise estimate E. This allows the filter to grow in noisy regions produced by large penumbrae without overblurring where shadows harden, near surfaces and under small light sources. It also means that our filter incurs almost no cost outside of the penumbrae. These are innovations over previous real-time denoising approaches.

We take advantage of the need to denoise two signals by computing them in the same pass with multiple render targets, which amortizes the cost of the G-buffer memory access for bilateral filter coefficients. We also directly write the final image in the vertical pass without explicitly computing the illumination-weighted shadow image WN = SN /UN , to reduce output bandwidth requirements. The visualizations shown in our result figures are for exposition and do not exist in memory in the actual implementation.

The final temporal integration is the game industry's standard practice for antialiasing [START_REF] Salvi | Anti-aliasing: Are we there yet?[END_REF]]. It reduces all forms of image aliasing, including any residual low frequency noise from shadow undersampling.

Noise Estimate

Total-variation noise estimation The choice of estimator E is key. At first, one might consider statistical variance of each neighborhood, but this is actually a poor choice because it ignores ordering. For example, the variance of a comb function (impulse train) and a step function can be the same, but in our application the comb is probably a noisy shadow that should be blurred and the step function a well-sampled hard-shadow edge that should be preserved. So, we instead designed an estimator that considers pixel ordering.

We compute E at each pixel x from the total variation of the color gradient of WN = S N U N estimated by n path integrals across a set of uniformly-angled lines { } through x:

E[x] = 1 n n i i ||d 2 WN [ i(t)]|| dt 2 dt. (29)
When both SN and UN are close to zero, arbitrarily let WN = 1.

E[x] is a robust measure for the amount of "wiggle" in the 2D neighborhood. We use an asterisk shape of n = 4 lines with a random per-pixel offset to avoid bias. The idea of total variation as a noise estimate was first introduced by Rudin et al. [1992], and has previously been used for denoising by solving a nonlinear optimization problem to minimize it. We instead use this measure to drive a joint-bilateral filter and work with derivatives. As desired, our method preserves both gradients and is much more efficient because it admits a closed-form implementation within a pixel shader, whereas an iterative solver would create execution divergence.

Because we estimate E from randomly chosen paths, it is slightly noisy itself. Rather than increasing n towards convergence, we simply box-filter the result. See the supplement for a GLSL implementation of the noise estimator.

SN E denoise[SN ] Final Image

Figure 9: Detail of denoising terms from the teapot scene.

Validation of our noise estimate

The definition of E in Equation (29) ensures that only penumbrae and not stochastic sampling of shading or materials are considered noise; after all, we have analytic shading results and this term is only computing shadowing. This point is subtle, but the implications for the difference from fullscreen denoising are profound in both potential quality and performance, because the filter cost is dominated by bandwidth required, which is proportional to the radius. 4.3 Real-Time Area Lights [START_REF] Heitz | Realtime polygonal-light shading with linearly transformed cosines[END_REF] give an efficient analytic solution suitable for use as U , for unshadowed direct illumination from polygonal area sources. We render a G-buffer and then compute U , UN , and SN in a single deferred shading pass. The terms used to compute the stochastic estimates are three textures that store the shadow ray origins (with start distance offset in the alpha channel), the shadow ray directions (with end distance in the alpha channel), and the radiance from the source along that ray. For a budget of N rays per pixel, these textures are N times larger than full-screen images.

Following the deferred shading pass, we cast all of the shadow rays in parallel. (Our implementation uses OptiX for peak performance; FireRays and Embree are viable alternatives). Then, we add the corresponding radiance contribution at each pixel to SN if the shadow ray is unoccluded.

The stochastic images UN and SN depend on the choice of rays, i.e., points on the light to match with the fixed surface point at each pixel. To compute these, we apply weighted multiple importance sampling to the light surface and the BRDF at each pixel, with a fixed ray budget per light per pixel. For each ray, we choose BRDF (vs. light) sampling with probability proportional to the analytic reflected illumination U . For BRDF samples, we choose glossy (vs. matte) lobe sampling with probability proportional to the glossy reflected illumination divided by the total reflected illumination (again from U ). For light samples, we use the solid angle sampling approach of Ureña et al. [2013]. To minimize clumping, we make all of these choices in a low-discrepancy fashion in screen space, so we use Halton sequences with 2, 3, 5, and 7 jittered by a 2D screen-space blue noise distribution as in [START_REF] Georgiev | Blue-noise dithered sampling[END_REF], for the four pseudo-random values required. See our supplement for the GPU implementation and normalization terms.

Results

Figure 8 shows the key terms, including the implicit conditional shadow image, for lights of increasing size casting shadows on cubes, using our area source application. Note the low-discrepancy pattern of the stochastic samples, and the preservation of texture and edge detail since only shadows and not the final image are denoised. The high-noise estimate is correctly tightly restricted to the noisy shadow edge for the left-most cube with the small light source and simultaneously correctly grows to fill the penumbra for the large source.

Figure 11 compares final images from two different state-of-the-art approaches: our area-light shadowing method, and a nearly converged reference image produced from 256 shadow rays per light per pixel. The state of the art for real-time rendering is analytic shading either with no shadows at all (left column) or with point shadows from the center of the light (right column). Point shadows can be produced by shadow maps or ray tracing; we use ray tracing here to avoid shadow map filtering issues. Our result (second column) produced with two shadow rays per light per pixel is closer to the reference image (third column) than these alternatives. The primary remaining error in our result is slight overblurring of some penumbrae due to aggressive denoising. Observe that material details, glossy highlights, and geometric edges are never blurred by our technique.

The first row shows a living room lit by two area sources that simulate indirect light. There is a slightly blue area light on the ceiling simulating reflected light that passed through the windows from the sky, and a slightly orange area light light off to the camera's left simulating indirect light from artificial sources inside the room.

The extremely soft shadows throughout the room resemble ambient occlusion. Compare this to the implausible result from point lights.

The second row shows a breakfast room lit by sunlight through the window on the right. There is also a dim environment probe in this room so that the backs of objects are not black. The shadows in our result correctly soften with distance from the window and the shading wraps around the curve of the chair and the pitcher.

The third row shows the Sponza atrium with the Lucy statue in it.

Lucy is backlit by a large area source the reflects in the glossy floor.

The "reflection" of the statue is actually only a shadow that masks the area light analytic reflection, but due to the contrast may read visually as indirect light. Again, the point light and unshadowed analytic results are poor approximations compared to ours. The varying smoothness of the marble and sharp edges of the light in reflection are preserved because we only denoise the shadow, not the light or material.

Table 1 gives timings for the GPU render passes of this scene. The essential point is that it produces reasonable area light shadows in milliseconds, versus seconds for convergence by stochastic methods. The GPU ray casting dominates our denoising and weighted shadow computations. GPU ray casting is both very fast now and on a strong acceleration trend (especially for shadow rays). As ray cast cost decreases in the future, the net cost the passes in Table 1 for high-quality area lights could soon be less than the cost of shadow map generation and filtering for medium-quality point lights.

The fourth row shows San Miguel at sunset. The sky is modeled as very large purple area light overhead and the sun is a smaller orange light off the left of the frame. In the video, the sun sets and the sky dims as the camera tracks. The sky shadows are extremely diffuse due to the size of that source, yet are stable and noise free in our result.

Conclusion and Future Work

We have shown that our formulation, based on the ratio estimator, allows us to combine real-time illumination techniques with raytraced shadows. We believe that a similar approach could be fruitful for other problems:

Adding shadows from dynamic objects to precomputed irradiance maps Irradiance (a.k.a. light) maps contain precomputed illumination that is typically precomputed via radiosity or path tracing, or computed at runtime with very infrequent updates. As a result, dynamic objects such as animated characters do not correctly shadow irradiance-mapped surfaces.

To correct this, consider a light map rendered with no final-bounce shadowing to be a converged, shadow-less illumination image U , and compute stochastic final-bounce illumination with (SN ) and without (UN ) shadows. Then, compute and apply the denoised weighted shadow image to incorporate dynamic shadowing.

BRDF-aware denoising Separately denoising glossy and matte illumination-weighted shadowing terms with separate noise esti-mates and thus kernel sizes would avoid overblurring glossy shadows at the cost of doubling the number of terms to be denoised.

Low-variance estimators for multiple-bounce path tracing, photon mapping, or virtual point lights The weighted "shadow" terminology follows from our explanation of the algorithm in terms of direct illumination. But nothing in our derivation precludes its application to indirect illumination, or requires the error in the U image to be brighter than the true result. We can apply the technique with U as any noise-free approximation of the illumination, UN as a stochastic estimator of U , and SN as a stochastic estimator of the true result, so long as the error has only low-frequency terms that are robust to denoising.

For a path tracer, we can compute U and UN inexpensively using only direct illumination and let SN capture a noisy estimate of the net global illumination. In this case, W will primarily contain magnitudes greater than 1 as its role is increasing illumination.

Figure 2 :

 2 Figure 2: Comparison of real-time options for incorporating shadows in the direct illumination of an area light. Options (d) and (e) are rendered at 2spp and denoised with the same bilateral filter.

Figure 6 :

 6 Figure6: Our ratio estimator. All results are rendered using 4 shadow rays per pixel and the denoiser is a bilateral filter. In contrast to a classic full-stochastic estimator, our ratio estimator renders noise-free results in unshadowed regions of the image. Furthermore, the denoising formulation of the ratio estimator only overblurs shadows and shading details remain sharp.

Figure 7 :

 7 Figure 7: Robustness for denoising. Estimators denoised by a σ = 4-pixel Gaussian. Denoising the estimators on a smooth surface produces acceptable results. With this normal map, denoising the full-stochastic estimator blurs the shading details, denoising the controlvariate estimator removes the shadow because it spreads negative values from the error correction terms, and denoising the ratio estimator produces the best result.

Figure 8 :

 8 Figure 8: Our real-time algorithm: computations, intermediate terms and result. The figures are at 1920×1080 each; zoom in to see stochastic sampling noise.

  1. Compute the noise estimate E from WN = S N U N 2. Denoise E itself via a 3 × 3 box filter 3. Horizontally joint-bilateral filter SN and UN simultaneously 4. Vertically joint-bilateral filter SN and UN simultaneously and produce a single output: final image U × denoise[S N ] denoise[U N ]

Figure 9

 9 Figure 9 shows detail near the top knob from the teapot scene of the numerator term undergoing denoising. Observe in the center that the E estimate accurately identifies the penumbra. The denoise[SN ] subimage is the denoised numerator. It shows that no blurring (and thus, minimal cost) occurs outside of the penumbra. The umbra on the back of the teapot near the bottom of the figure and the unshadowed highlight remain untouched. Yet, the highlight is not noisy in the denoised weighted shadow or the final image Figure 2(e) because that stochastic noise is identical in the denominator denoise[UN ] and thus cancels. See Figure 10 for more examples of accurate noise estimation.

Figure 10 :Figure 11 :

 1011 Figure 10: Validation of our noise estimate. Our noise estimate detects regions where the stochastic weighted-shadow image WN is noisy and is insensitive if only SN is noisy. Thanks to this property, we use large denoising kernels only in shadowed regions of the image. The figures are at 1920×1080 each; zoom in to see stochastic sampling noise.

Table 1 :

 1 Pass timings with a GeForce 1080 at 1920×1080 on Sponza, using 16-bit float precision for all buffers and N shadow rays per pixel.

	Shade Pass	Time Denoise Pass	Time
	Generate rays	N × 0.36 ms Estimate noise	2.12 ms
	OptiX ray cast N × 3.68 ms Horizontal bilateral 0.70 ms
	Accumulate	N × 0.10 ms Vertical bilateral	0.70 ms


Linfeng Ding
email: dingli@uni-mainz.de

Stephan Buhre

Clemens Kunisch

Boris Kaus

Pressure dependence of density and structural relaxation of glass near the glass transition region

Keywords: Glass, High pressure, Piston-cylinder apparatus, Structural relaxation, Fictive temperature, Fictive pressure

A simplified and effective pressure cell together with an experimental procedure has been developed and applied to compress samples of SCHOTT N-BK7

 glass under static high pressures in a piston-cylinder apparatus. Results from the density and volume recovery measurements show that, the glass samples were effectively densified in piston-cylinder apparatus with the density at room temperature increasing linearly with frozen-in pressure. To explain the experimental data, we developed a mathematical model based on a suggestion by Gupta (1988) with two internal parameters, named fictive temperature ( ) and fictive pressure ( ), which fits the experimental data well.

Introduction

It is well known that the properties of glass at room temperature are influenced by temperature, particularly around the glass transition region [START_REF] Ritland | Limitations of the fictive temperature concept[END_REF] . Glass is not in equilibrium, but instead, the structure and properties of glass are time-dependent while it aims to reach equilibrium [START_REF] Scherer | Volume relaxation far from equilibrium[END_REF][START_REF] Welch | Dynamics of glass relaxation at room temperature[END_REF] . The temperature and time-dependence of properties as well as the structural relaxation of glass in the glass transition region have been well-characterized experimentally and can be explained with several mathematical models [START_REF] Ritland | Limitations of the fictive temperature concept[END_REF][START_REF] Scherer | Volume relaxation far from equilibrium[END_REF][START_REF] Tool | Relation between inelastic deformability and thermal expansion of glass in its annealing range[END_REF][START_REF] Narayanaswamy | A model of structural relaxation in glass[END_REF][START_REF] Mauro | Nonequilibrium viscosity of glass[END_REF] . Yet, those studies were all performed for 1 bar, and did not consider the effect of increases in pressure.

However, pressure also changes the properties and the behavior of glass even at room temperature. Experimental examples are hardness measurement [START_REF] Taylor | Plastic deformation of optical glass[END_REF][START_REF] Shang | Creep behavior of soda-lime glass in the 100-500 k temperature range by indentation creep test[END_REF] and scribing of glass. In the glass transition range, the viscosity strongly depends on density [START_REF] Cook | Pressure and temperature dependent viscosity of two glass forming liquids: Glycerol and dibutyl phthalate[END_REF] and thus on pressure. For example, Bridgman and Simon [START_REF] Bridgman | Effects of very high pressures on glass[END_REF] compressed thin samples of oxide glasses (SiO 2 and B 2 O 3 ) and several mixed silicate glasses samples between two flat dies (uniaxial) made of carboloy at room temperature, and demonstrated that the vitreous SiO 2 had a density increase of almost 7.5% at 20 GPa, whereas the vitreous B 2 O 3 has an asymptotic increase of density near 6% at pressures close to 20 GPa. An experimental study on a borosilicate glass with pressure of up to 0.66 GPa and temperatures up to 300°C, showed that the density change was very small with only 0.05% at 0.66 GPa, 200°C and 24 hours [START_REF] Anderson | Effect of pressure on glass structure[END_REF] . A compilation of experimental data shows that there are significant differences in the manner in which different glass-types react to increases in pressure [START_REF] Rouxel | Poisson's ratio and the densification of glass under high pressure[END_REF] . Jin et al. [START_REF] Jin | Pressure effect on the structural relaxation and glass transition in metallic glasses[END_REF] reported a 6 K/GPa change in of a metallic glass by DSC measurements, and activation volume was suggested to describe the pressure dependence of glass relaxation in a free volume model.

Several phenomenological models [START_REF] Tribone | Pressure-jump volume-relaxation studies of polystyrene in the glass-transition region[END_REF][START_REF] Hodge | Enthalpy relaxation and recovery in amorphous materials[END_REF][START_REF] Simon | Enthalpy recovery of a glass-forming liquid constrained in a nanoporous matrix: Negative pressure effects[END_REF][START_REF] Grassia | Modeling volume relaxation of amorphous polymers: Modification of the equation for the relaxation time in the kahr model[END_REF] describe the time dependent frozen-in state and departure from equilibrium () [START_REF] Williams | Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function[END_REF] of glass forming materials depending on temperature history [START_REF] Tool | Relation between inelastic deformability and thermal expansion of glass in its annealing range[END_REF][START_REF] Narayanaswamy | A model of structural relaxation in glass[END_REF][START_REF] Moynihan | Structural relaxation in vitreous materials[END_REF] and also pressure history [START_REF] Tribone | Pressure-jump volume-relaxation studies of polystyrene in the glass-transition region[END_REF][START_REF] Hodge | Enthalpy relaxation and recovery in amorphous materials[END_REF][START_REF] Simon | Enthalpy recovery of a glass-forming liquid constrained in a nanoporous matrix: Negative pressure effects[END_REF][START_REF] Grassia | Modeling volume relaxation of amorphous polymers: Modification of the equation for the relaxation time in the kahr model[END_REF] . Those models are based on the Tool-Narayanaswamy-Moynihan [START_REF] Tool | Relation between inelastic deformability and thermal expansion of glass in its annealing range[END_REF][START_REF] Narayanaswamy | A model of structural relaxation in glass[END_REF][START_REF] Moynihan | Structural relaxation in vitreous materials[END_REF] model and use one internal parameter named fictive temperature ( ). On the other hand, Davies and Jones 20 pointed out that fictive pressure ( ) could also be an internal parameter. Gupta [START_REF] Gupta | Fictive pressure effects in structural relaxation[END_REF] introduced the to his structural relaxation model as an internal parameter together with , but has not compared his two-internal-parameter model to experimental data.

Here, we present new experimental data to study the density and structural relaxation of glass under elevated pressure of up to 1.5 GPa and elevated temperature, which were obtained using a piston-cylinder apparatus. We describe the new experimental procedure developed for this as well as a two-internal-parameter model that can overall fit the data.

Experimental procedure

Samples

The SCHOTT N-BK7  glass samples had been originally cooled at a rate of -0.4

K/h at ambient pressure and machined to cylinders with diameters of 7.400.02 mm and heights of 6.000.02 mm for the experiments. The density of those slow-cooled samples at room temperature is 2.50820.0006 g/cm [START_REF] Welch | Dynamics of glass relaxation at room temperature[END_REF] . The glass transition temperature T g =561 °C is taken as the temperature at which the viscosity is 10 12 Pa•s. Specially, the sample relatively fast-cooled at a rate of -120 K/h at ambient pressure in the lab was used for comparison to the samples after high pressure treatment, as all the samples in high pressure treatment were cooled at a rate of -120 K/h. The density of this fast-cooled sample at room temperature is 2.49570.0007 g/cm 3 .

High pressure treatment

Test equipment: A self-assembled piston-cylinder apparatus [START_REF] Ziaja | Metamorphism and melting of picritic crust in the early earth[END_REF] with an improved design based on the principle of Boyd and England [START_REF] Boyd | Apparatus for phase-equilibrium measurements at pressures up to 50-kilobars and temperatures up to 1750-degrees-c[END_REF] with uniaxial pressure up to 3 GPa was applied for the high pressure treatment. The uncertainty of pressure in the piston-cylinder apparatus was 0.05 GPa.

Pressure cell development:

A special pressure cell needed to be developed, to obtain glass samples with no cracks at the end of the experiments such that we obtain sufficient material for further measurements. If not, during subsequent measurements, the cracks in the sample may open and increase the apparent macroscopic volume of the sample, which erroneously would be interpreted as a decrease in density.

Experimental process: The experimental processes of high pressure treatment are shown in Fig. 1 (a). Glass samples are first subjected to high pressures and then heated up to 700 °C at 600 K/h. The temperature is chosen for the glass to be viscous enough for a hydrostatic pressure state to build up and to avoid crystallization. After keeping the temperature fixed for 24 hours (2 hours for the 1 GPa sample at the first experiment), the samples are relatively fast-cooled at -120 K/h to room temperature. Pressure release is done at a rate of -0.2 GPa/h to 0.5 GPa and further at a lower rate of -0.05 GPa/h to 0.1 MPa (-1.5 GPa/h to 0.1 MPa for 1.25 GPa sample due to the programming error,

purple line with open stars in Fig. 1 (a)). Glass samples reach stage ① after this process. Stress condition: The stress state in the glass sample in the piston-cylinder apparatus is influenced by gaps in the pressure cell that develop during assembling, especially the gap from sample to the other parts, the stiffness and Poisson's ratio of the insulator [START_REF] Rouxel | Poisson's ratio and the densification of glass under high pressure[END_REF] and by the thermal expansion coefficients of the materials used for the pressure cell. To make the glass sample more deformable than the surrounding materials, it is heated up to a temperature of 700 °C at which it behaves mainly viscous such that the glass sample obtains a hydrostatic stress state. During cooling, the glass changes from mainly viscous to frozen-in, during which the hydrostatic stress state may change into a uniaxial stress state. Therefore, the frozen-in stress state is not known precisely.

For the calculations, we assume a hydrostatic stress state always. Compared to this uncertainty in the overall stress state, the stress introduced to the sample by the frozen-in temperature gradient during cooling is small.

Density measurement

Test equipment: The Mettler Toledo-XP205 was used for the density measurement using the Archimedes principle. The density uncertainty was calculated by the statistical uncertainty from the 95% confidence interval of the student's t-test and the standard deviation.

Experimental process: Firstly, the mass of the samples was measured in air, then the mass of samples in the liquid was measured and the temperature of the liquid was recorded at the same time. Each sample was measured 5 times and the results reported are the average values rounded to four decimals.

Volume recovery measurement

Test equipment: The Bähr DIL 806 optical dilatometer is used to measure the length change of the sample diameter without mechanical contact. The temperature of the dilatometer has been calibrated by a pure quartz crystal with a repeatability of α-β transformation temperature of 0.39 °C in the same position.

Experimental process: The change in diameter was measured by applying the temperature history shown in Fig. 1 (b). The dilatometer was heated from room temperature to 500 °C at a rate of +1200 K/h, then at +120 K/h to 525 °C and holding for 72 hours. At the end of this isothermal hold the sample diameter does not change any more. This equilibrium state is denoted as stage ②. Afterwards, the samples were freely cooling down to the room temperature, at which they reached stage ③.

Results

High pressure treatment

The standard pressure cell (1#) for rock experiment resulted in a glass sample broken to many pieces. Thus, we designed and prepared the other four pressure cells seeking for a suitable pressure cell for the high pressure experiments to obtain glass samples with no crack and no blackening. The designed pressure cells and accordingly images of glass samples after high pressure treatment are shown in Fig. 2. For further measurements, we used the pressure cell 2#, because 2# was the only pressure cell which obtained glass sample without cracks and blackening. The further compression treatment at 0.75 GPa and 1.1 GPa resulted in only small amounts of cracks at the surface and no blackening inside (similar to the 0.5 GPa sample).

Moreover, the slight black film on the 1.1 GPa sample disappeared during the volume recovery measurement. Therefore, the pressure cell 2# is suitable for compressing glass.

Density after high pressure treatment

In the high pressure treatment, the glass samples were cooled at the fast cooling rate of -120 K/h to room temperature (Fig. 3, stage ①). Compared to the fast-cooled uncompressed sample with density of 2.49570.0007 g/cm [START_REF] Welch | Dynamics of glass relaxation at room temperature[END_REF] , their density increases approximately linearly with pressure, which can be fitted with:

(

where is the density of the sample at room temperature in g/cm [START_REF] Welch | Dynamics of glass relaxation at room temperature[END_REF] , P is the applied pressure in GPa. Specifically, the density increased by 4.49% after the 1.5 GPa treatment. The sample compressed at 1 GPa (red square in Fig. 3) broke into several parts which led to a very high uncertainty.

Volume recovery

The measured relative length changes in diameter during volume recovery of the uncompressed samples were smaller in magnitude than those of the compressed samples (Fig. 4), which go from an increase in diameter of 1.03 % after 0.5 GPa compression to 1.89 % after 1.5 GPa compression. Particularly, the two arrows in Fig. 4 indicate that the 1.1 GPa curve shows fluctuations after reaching equilibrium and lost data points after approximately 37 h (the sample was broken up which caused measurement error). The difference between these two arrows in relative length change is only approximately 0.02%, which leads to an uncertainty in specific volume of 0.0003 cm 3 /g. This sample does not reach 72 h, therefore does not contribute to the average specific volume.

During the volume recovery, the more the samples have been compressed in the high pressure treatment, the more the measured change in relative length increases.

Mathematical Model

To describe the measured volume recovery, we firstly applied the exponential decay model [START_REF] Welch | Dynamics of glass relaxation at room temperature[END_REF][START_REF] Williams | Non-symmetrical dielectric relaxation behaviour arising from a simple empirical decay function[END_REF][START_REF] Phillips | Stretched exponential relaxation in molecular and electronic glasses[END_REF][START_REF] Koontz | Prony series spectra of structural relaxation in n-bk7 for finite element modeling[END_REF] , which fits the data well but reached different relaxation times in the new equilibrium state at 525 °C and 1 bar. A new mathematical model is developed following the suggestion by Gupta 21 , in which the volume relaxation depends on two internal parameters, named fictive temperature ( ) [START_REF] Tool | Relation between inelastic deformability and thermal expansion of glass in its annealing range[END_REF] and fictive pressure ( ) [START_REF] Gupta | Fictive pressure effects in structural relaxation[END_REF] , to account for the thermal history and the pressure history.

The mathematical model assumes a true hydrostatic stress during the whole high pressure treatment. The specific volume depends on T, ,P, according to

(2) with =525 °C, =0.0001 GPa, and =0.4065 cm 3 /g. The rates of change in fictive temperature and fictive pressure

=(L TT L TP V )( )( (3)
=(L TP L PP V )( )( (4)
are taken from Gupta [START_REF] Gupta | Fictive pressure effects in structural relaxation[END_REF] . Here, is the specific heat capacity of the liquid (equilibrium) state and the frozen (glassy) state . Similarly, is the difference in thermal expansion coefficient between liquid and frozen glass, is the difference in compressibility between liquid glass and frozen glass with (K is the bulk modulus). , , and are kinetic coefficients [START_REF] Gupta | Fictive pressure effects in structural relaxation[END_REF] which depend on T, P, and . For the model in this report, the assumptions , and are used. This introduces three time constants: for changes in specific heat, for changes in compressibility, and for changes in thermal expansion coefficient. Those three different time constants

differ from the characteristic time for changes in the shear modulus by constant factors , , only.

Following the idea of Markovsky and Soules [START_REF] Markovsky | An efficient and stable algorithm for calculating fictive temperatures[END_REF] , the and at the end of a time step t are calculated from the and at the beginning of the same time step t.

(8) (9)
In the model of this report, the characteristic time constant for changes of the shear modulus (10) depends on real pressure and fictive pressure via an apparent activation volume [START_REF] Jin | Pressure effect on the structural relaxation and glass transition in metallic glasses[END_REF][START_REF] Karato | Deformation of earth materials: An introduction to the rheology of solid earth[END_REF] ( ) with a percentage factor of (0≤ ≤1) and . In case of , the characteristic time (T, ,P, ) reduces to the characteristic time in the Tool-Narayanaswamy 4, 5 model (

where is the fictive temperature and x is a fixed factor.

This mathematical model captures the influence of the pressure pre-history on the volume recovery measurements at 525 °C with a single set of parameters (see figure caption for employed constants in Fig. 5 (a)). The calculation for each sample comprises the initial cooling at fast or slow cooling rate, the high pressure treatment (if any), the heating-up stage and, finally, the isothermal volume recovery measurement.

Discussion

In our experiments, sample blackening occurred during pressurization of silicate glasses in all pressure cells using h-BN as pressure transmitting medium. The mechanism of blackening is not clear. Brooker et al. [START_REF] Brooker | Reduction in piston-cylinder experiments: The detection of carbon infiltration into platinum capsules[END_REF] reported that it may be caused by carbon infiltration [START_REF] Brooker | Reduction in piston-cylinder experiments: The detection of carbon infiltration into platinum capsules[END_REF] from the graphite furnace. The blackening of 0.9 GPa sample was not reversible at 525 °C and showed no absorption bands measured in the transmission spectrum. On the other hand, the sample from the 1.1 GPa experiment (no h-BN) became visibly clear after heating up again in air to 525 °C during the volume recovery, which may thus be due to the reduction of oxygen [START_REF] Ratso | Highly active nitrogen-doped few-layer graphene/carbon nanotube composite electrocatalyst for oxygen reduction reaction in alkaline media[END_REF] in the glass and different from the nonreversible blackening.

We note that the shapes of length change curves in Fig. 4 are not exactly the same although the overall change is quite regular. This is attributed to the measurement accuracy of the optical dilatometer. The volume recovery shows that, at sufficiently high temperature and holding time, the densification of glass under high pressure is reversible [START_REF] Anderson | Effect of pressure on glass structure[END_REF] . The isothermal hold of 72 h at 525 °C is sufficient for SCHOTT N-BK7

Fig. 1

 1 Fig. 1 (a) Process of high pressure treatment in piston-cylinder apparatus; (b) process of volume recovery measurements in dilatometer for all the samples. Arrow in (a) shows the time sequence of high pressure treatment. Stage ①: samples compressed by high pressure and before volume recovery measurement; Stage ②: samples reach

Fig. 2

 2 Fig. 2 Pressure cells development and images of glass samples after high pressure treatment for different pressure cells. 1-Pyrophyllite; 2-Steel plug; 3-Thermocouple; 4-Graphite; 5-NaCl; 6-Sample; 7-Pt; 8-Al 2 O 3 ; 9-Borosilicate glass; 10-h-BN.

Fig. 3

 3 Fig. 3 Density of glass samples at three different stages: Stage ① = at room temperature after high pressure treatment; Stage ② = at 525 °C, at the end of the

Fig. 4

 4 Fig. 4 Relative length change of samples versus time at volume recovery temperature of 525 °C.

Fig. 5

 5 Fig. 5 (a) Measured volume recovery data fitted by the new model with H/R=74000 K, x=0.6, =2, =3, =20, =20 GPa, =1. (b) Change of and of fast-cooled samples versus time in our mathematic model.

Fig. 6

 6 Fig. 6 and of fast-cooled samples at the beginning of the isothermal volume recovery at 525 °C as a function of fictive pressure in our mathematic model.



  Fig.7 (a)) after 1 GPa high pressure treatment (with higher at the beginning) are

Fig. 7

 7 Fig. 7 (a) Prediction of the volume recovery at a constant heating rate after various temperature and pressure treatments from our mathematic model. (b) Corresponding change of and in our mathematic model.

Acknowledgement

This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 642029 -ITN CREEP. The authors appreciate the constructive comments from two anonymous reviewers. The authors thank Ms. Silvia Weitzel for the help on the volume recovery measurements, Ms. Andrea Baltes on the density measurements, and Mr. Harald Schmück on the high pressure treatment.

Conclusions

The effects of high pressure on density and volume recovery of SCHOTT N-BK7  glass near the glass transition region has been measured and quantitatively described in this work.

(1) A simplified and effective pressure cell (2#) together with an experimental process is developed to compress glass at high pressures around the transition region using a piston-cylinder apparatus.

(2) Glass samples can be effectively densified in the piston cylinder apparatus with an approximately linear increase in density at increasing pressure.

(3) A mathematical model with two internal parameters is developed that describes the volume recovery of all samples with a single set of parameters.


Andrei Novikov

On the estimation of the expectation of a heavy-tail distribution

Keywords: Extreme value theory, Heavy-tail distributions, Hill estimator, Regular variation, Tail index

We use bias-reduced estimators of extreme quantiles from heavy-tail distributions, to build a new estimator of the expectation in the case of infinite second moment. The asymptotic normality of the proposed estimator is established under mild conditions.

Introduction

Let X 1 , X 2 , ... be independent and identically distributed non-negative random variables with expectation µ < ∞, variance σ 2 and cumulative distribution function (cdf) F. Suppose that the tail of F is regularly varying at infinity with tail index (-α) < 0, that is

lim t→∞ 1 -F (tx) 1 -F (t) = x -α , for any x > 0, (1)
Such heavy-tail distributions include distributions such as Pareto, Burr, Student, α-stable (0 < α < 2) , and log-gamma, which are known to be appropriate models for fitting large insurance claims, large fluctuations of prices, log-returns, etc. We are concerned with the construction of a bias-reduced asymptotically normal estimator for the expectation

µ := ∞ 0 xdF (x) ,
which could be rewritten, in terms of the quantile function (corresponding to the cdf F )

Q(s) := inf {x : F (x) ≥ s} , 0 < s < 1, as µ = 1 0 Q (1 -s) ds. (2)
For a given sample X 1 , ..., X n , let

Q n (s) := inf {x ∈ R : F n (x) ≥ s} , 0 < s ≤ 1,
denote the sample quantile function (classical non-parametric estimator of Q) associated to the empirical cdf defined on the real line by F n (x) := n -1 n i=1 I (X i ≤ x) , with I (•) being the indicator function. The natural (unbiased) estimator of µ is the sample mean

1 0 Q n (1 -s) ds = 1 n n i=1 X i =: X n . (3)
From the Central Limit Theorem (CLT), the sequence of r.v.'s √ n X nµ /σ, n ≥ 1 converges in distribution to the standard Gaussian r.v., provided that the second-order moment E X 2 1 is finite. This is a very restrictive condition in the context of heavy-tail distributions as the following considerations show. Assume that the r.v. X 1 follows the Pareto law with index α > 0, that is, 1 -F (x) = x -α for x ≥ 1. When α > 1, the expectation µ exists, but E X 2 1 is only finite for α ≥ 2. Hence, the range α ∈ (1, 2) is not covered by the CLT and thus we need to seek another approach to handle this situation. Making use of Weissman's estimator of high quantiles an alternative estimator for µ was proposed and its asymptotic normality was established for any α ∈ (1, 2). Let us define the following estimator for Q :

Q n (1 -s) :=    Q W n (1 -s) for 0 < s < k/n Q n (1 -s) for k/n ≤ s < 1,
where

Q W n (1 -s) := (k/n) 1/ α H n X n-k,n s -1/ α H n , s ↓ 0 (4)
is Weissman's estimator of high quantiles, with

α H n := k -1 k i=1 log X n-i+1,n -log X n-k,n -1 , (5)
being the well-known Hill estimator of the tail index α, and X 1,n ≤ ... ≤ X n,n denoting the order statistics pertaining to the sample X 1 , ..., X n . The number k represents the number of upper order statistics used in the computation of α H n , it is an integer sequence k = k n satisfying

1 < k < n, k → ∞ and k/n → 0 as n → ∞. (6)
By replacing

Q (1 -s) by Q n (1 -t) in formula (2)
, an alternative estimator for µ is as follows:

µ P n = µ P n (k) := 1 0 Q n (1 -s)ds = k/n 0 Q W n (1 -s) ds + 1 k/n Q n (1 -s)ds,
which, by a straightforward calculation, is equal to

µ P n := k n α H n α H n -1 X n-k,n + 1 n n i=k+1 X n-i+1,n , (7)
provided that α H n > 1. Moreover, the same author showed that, under suitable regularity assumptions, for any α ∈ (1, 2),

√ n µ P n -µ k/nX n-k,n d → N 0, σ 2 (α) , as n → ∞, (8)
where

σ 2 (α) := α/ (1 -α) 4 (2 -α) .
Throughout this paper, the standard notations Actually, this estimator is defined in the more general situation where the r.v. X is real (not necessarily non-negative) with lower and upper heavy tails. He simultaneously took into account the regular variations of both tails of G and the balance condition

lim t→∞ (1 -F (t)) / (1 -F (t) -F (-t)) = p ∈ [0, 1] .
In this paper, we only consider non-negative r.v.'s. Our motivation comes from the actuarial risk theory where insurance losses are represented by such r.v.'s. In this case, µ P n may be interpreted as an estimator of a risk measure called the net premium. Note that in our case, since r.v. X is non-negative, we have F (-x) = 0 for x ≥ 0, which yields p = 1 in the above balance condition.

Hill's estimator α H n plays a pivotal role in statistical inference on distribution tails. This estimator has been thoroughly studied, improved and even generalized to any real parameter α. Weak consistency of α H n was established assuming only that the underlying cdf F satisfies condition (1). The asymptotic normality of α H n has been established under the following stricter condition that characterizes Hall's model.

1 -F (x) = cx -α + dx -β + o x -β , as x → ∞, (9)
for some c > 0, d = 0 and β > α > 0. Note that [START_REF] Bouchard | Nonparametric frontier estimation by linear programming[END_REF], which is a special case of a more general second-order regular variation condition, is equivalent to

Q (1 -s) = c 1/α s -1/α 1 + α -1 c -β/α ds β/α-1 + o (1) , as s ↓ 0. ( 10
)
The constants α and β are called, respectively, first-order (tail index, shape parameter) and second-order parameters of cdf F.

Extreme value based estimators essentially rely on the number k of upper order statistics involved in estimate computation. Hill's estimator has, in general, a substantial variance for small values of k and a considerable bias for large values of k. Hence, one has to look for a k value, denoted by k * , that balances between these two vices. The choice of this optimal value k * represents a thorny issue in the process of estimating the tail index and related quantities. To solve this problem, several adaptive procedures are available. A theoretical optimal choice of k is obtained by minimizing the asymptotic mean squared error (RMSE) of α H n . Indeed, under condition (9), we have

k * := 2 -1 αβ 2 (β -α) -3 d -2 c 2β/α α 2β-α n 2β-2α 2β-α . (11)
Though Peng's estimator µ P n enjoys the asymptotic normality property, it still has a problem due to the fact that, it is based on Weissman's estimator Q W n known to be largely biased. Fortunately, many estimators with reduced biases are proposed in the literature as an alternative to Q W n .

In this paper, we use the bias-reduced estimator of the high quantile Q (1s) , recently proposed by who exploited the censored maximum likelihood (CML) based estimators α, β of the couple of regular variation parameters (α, β). The CML estimators α, β are defined as the solution of the two equations (under the constraint

β > α H n ) 1 k k i=1 1 G i (α, β) = 1 and 1 k k i=1 1 G i (α, β) log X n-i+1,n X n-k,n = β -1 , (12)
where

G i (α, β) = α β 1 + αβ α -β H (α) X n-i+1,n X n-k,n β-α - αβ α -β H (α) , (13)
with

H (α) = 1 α - 1 k k i=1 log X n-i+1,n X n-k,n .
The bias-reduced estimators

Q LP Y n (1 -s) , of the high quantiles Q (1 -s)
, are obtained by substituting α, β to (α, β) in [START_REF] Byström | Managing extreme risks in tranquil and volatile markets using conditional extreme value theory[END_REF]. That is

Q LP Y n (1 -s) := c 1/ α s -1/ α 1 + α -1 c -β/ α ds β/ α-1 , s ↓ 0, (14)
where

               c = α β α -β k n X α n-k,n 1 β - 1 k k i=1 log X n-i+1,n X n-k,n , d = α β β -α k n X β n-k,n 1 α - 1 k k i=1 log X n-i+1,n X n-k,n . (15)
The consistency and asymptotic normality of Q LP Y n (1s) are established by the same authors. Now we can define another estimator for the quantile function Q as follows:

Q n (1 -s) =    Q LP Y n (1 -s) for 0 < s < k/n Q n (1 -s) for k/n ≤ s < 1.
By replacing Q by Q n , in formula (2), we get

µ n = µ n (k) := 1 0 Q n (1 -s)ds = k/n 0 Q LP Y n (1 -s) ds + 1 k/n Q n (1 -s)ds. (16)
An elementary integral calculation leads to a new bias-reduced estimator for µ defined by the following formula:

µ n := (k/n) (n c/k) 1/ α α α -1 + d c -β/ α (k/n) β/ α-1 β -1 + 1 n n i=k+1 X n-i+1,n , (17)
provided that β > α > 1 so that µ n be finite.

The rest of the paper is organized as follows. In Section 2, we briefly discuss the third ordercondition of regular variation before establishing the asymptotic normality of µ n . Some concluding remarks notes made in Section 4. Finally, some of the main results used in Section 3 are gathered in the Appendix.

Main results

In the theory of extremes, a function, denoted by U and (sometimes) called tail quantile function, is used quite often. It is defined by

U (t) := (1/ (1 -F )) -1 (t) = Q (1 -1/t) , 1 < t < ∞.
In terms of this function, Hall's conditions ( 9) and ( 10) are equivalent to

U (t) = c 1/α t 1/α 1 + α -1 c -β/α dt 1-β/α + o (1) , t → ∞. (18)
This implies that

lim t→∞ log [U (tx) /U (t)] -α -1 log x A 1 (t) = x 1-β/α -1 1 -β/α , for any x > 0, (19)
where

A 1 (t) := dα -1 (1 -β/α) c -β/α t 1-β/α .
The function A 1 (t) , which tends to zero as t → ∞ (because β > α), determines the rate of convergence of log [U (tx) /U (t)] to its limit α -1 log x. Relation ( 19) is known as the second-order condition of regular variation.

Unfortunately, the second-order regular variation is not sufficient to find asymptotic distributions for the estimators defined by the systems ( 12) and [START_REF] Daouia | Functional convergence of quantile-type frontiers with application to parametric approximations[END_REF]. We strengthen it into a condition, called third-order condition of regular variation and given by [START_REF] Daouia | Kernel estimators of extreme level curves[END_REF], that specifies the rate of [START_REF] Daouia | On kernel smoothing for extremal quantile regression[END_REF].

lim t→∞ log [U (tx) /U (t)] -α -1 log x A 1 (t) - x 1-β/α -1 1 -β/α A 2 (t) = D (α, β, ρ) , (20)
where A 2 (t) → 0 as t → ∞, with constant sign near infinity and

D (α, β, ρ) := 1 ρ x 1-β/α+ρ -1 1 -β/α + ρ - x 1-β/α -1 1 -β/α ,
with ρ being a positive constant called third-order parameter. The asymptotic normality of α, β and c was established under the following extract conditions on the sample fraction k, as n → ∞,

(i) √ k |A 1 (n/k)| → ∞, (ii) √ kA 2 1 (n/k) → 0, (iii) √ kA 1 (n/k) A 2 (n/k) → 0. ( 21
)
As for d, it is asymptotically normal under the assumption

√ k |A 1 (n/k)| / log (n/k) → ∞ added to (ii) and (iii) .
Example 1 Consider the Fréchet cdf with shape parameter α > 0

F (x) = exp -x -α , x > 0. ( 22
)
The corresponding tail quantile function is defined by

U (t) = (-log(1 -1/t)) -1α , for t > 1.
Applying Taylor's expansion (to the third order) to U and identifying with [START_REF] Daouia | Nadarayas estimates for large quantiles and free disposal support curves[END_REF], yield β = 2α, c = 1 and d = -1/2. The condition (20) holds for A 1 (t) = t -1 /2α, A 2 (t) = (α -3) t -2 /12α 2 and ρ = 3α.

Note that, from a theoretical point of view, assumptions ( 6) and ( 21) are realistic, as the following example shows. Indeed, let us choose k = [n ǫ ] , 0 < ǫ < 1, then it easy to verify that these assumptions hold for any 2/3 < ǫ < 4/5. The notation [•] stands for the integer part of real numbers.

Our main result, namely the asymptotic normality of the bias-reduced estimator µ n , is formulated in the last of the following four theorems. In Theorem 1, we give an approximation of α in terms of Brownian bridges, which leads to its asymptotic normality stated in Theorem 2. We do the same thing to µ n in Theorem 3. But, this does not meet our needs to achieve the major object of this paper. Then, we need to approximate both α and µ n by linear functional of the same sequence of standard Brownian bridges B n (s) .

Theorem 1 Assume that the third order condition (20) holds with β/α =: λ > 1 and let k = k n be an integer sequence satisfying ( 6) and [START_REF] Daouia | Robustness and inference in nonparametric partial frontier modeling[END_REF]. Then there exists a sequence of Brownian bridges

{B n (s) , 0 ≤ s ≤ 1} such that √ k ( α -α) = α (η 1 W 1n + η 2 W 2n + η 3 W 3n ) + o p (1) , as n → ∞
where W 1n , W 2n and W 3n are sequences of centered Gaussian r.v.'s defined by

W 1n := n/kB n (1 -k/n) -n/k 1 0 s -1 B n (1 -ks/n) ds, W 2n := λ -1 -1 n/kB n (1 -k/n) + (λ -1) n/k 1 0 s λ-2 B n (1 -ks/n) ds, W 3n := (1 -λ) n/k 1 0 s λ-2 (log s) B n (1 -ks/n) ds + λ -2 n/kB n (1 -k/n) -n/k 1 0 s λ-2 B n (1 -ks/n) ds,
and

η 1 := λ 4 (λ -1) 4 , η 2 := λ 2 (2λ -1) (3λ -1) (λ -1) 5 η 3 := λ 3 (2λ -1) (λ -1) 4 2 .
Theorem 2 Under the assumptions of Theorem 1, we have

√ k ( α -α) d → N 0, α 2 β 4 / (α -β) 4 , as n → ∞. (23)
Theorem 3 Under the assumptions of Theorem 1, we have, as

n → ∞ √ n k/n (nc/k) 1/α { µ n -µ} = - α (α -1) 2 {η 1 W 1n + η 2 W 2n + η 3 W 3n } + W 4n + o p (1) ,
where W 1n , W 2n and W 3n are those of Theorem 1 and

W 4n := - 1 k/n B n (1 -s) dQ (1 -s) k/n (nc/k) 1/α .
Theorem 4 Under the assumptions of Theorem 1, we have

√ n k/n (nc/k) 1/α { µ n -µ} d → N 0, σ 2 (α, β) , as n → ∞, (24)
where

σ 2 (α, β) := α 2 β 4 (α -1) 4 (α -β) 4 + 2 2 -α + 2αβ 2 (α -1) 2 (α -β) 2 . ( 25
)
The following corollary to Theorem 4 provides a straightforward practical way to build confidence intervals for µ.

Corollary 1 Under the assumptions of Theorem 1, we have

√ n k/nσ α, β (n c/k) 1/ α { µ n -µ} d → N (0, 1) , as n → ∞,
where α, β and c are the estimates of α, β and c given in ( 12) and ( 15) respectively.

Proofs

3.1 Proof of Theorem 1.

First recall that λ = β/α > 1. Then, we have, as n → ∞

α -α = α η 1 (S 1 -1) + η 2 S 2 -λ -1 + η 3 S 3 -λ -2 + o p k -1/2 , (26)
where

S 1 := 1 k k i=1 log Y n-i+1,n Y n-k,n , S 2 := 1 k k i=1 Y n-i+1,n Y n-k,n 1-λ , (27)
and

S 3 := 1 k k i=1 Y n-i+1,n Y n-k,n 1-λ log Y n-i+1,n Y n-k , n , (28)
with Y 1,n ≤ ... ≤ Y n,n being the order statistics pertaining to a sample Y 1 , ..., Y n of independent and identically distributed non-negative random variables, defined on the same probability space as the

X ′ i s, with cdf G (y) = 1 -y -1 , for y > 1. (29)
A probability space (Ω, A,P) may be constructed, carrying an infinite sequence U 1 , U 2 , ... of independent (0, 1) -uniform r.v.'s and a sequence of Brownian bridges {B n (s) , 0 ≤ s ≤ 1} , n = 1, 2, ..., having, amongst others, the property stated in Lemma 1. Let U 1,n ≤ ... ≤ U n,n denote the order statistics pertaining to U 1 , ..., U n and define the empirical quantile function V n (s) as

V n (s) = U i,n for (i -1) /n < s ≤ i/n, i = 1, ..., n, and V n (0) = U 1,n .
Lemma 1 On this probability space, for every 0 ≤ τ < 1/2, we have, as n → ∞ sup

1/n≤s≤1-1/n | √ n (s -V n (s)) -B n (s)| (s (1 -s)) 1/2-τ = O p n -τ . ( 30
)
Without loss of generality, we assume that

Y i = G -1 (U i ) = (1 -U i ) -1 , i = 1, ..., n, and
Y i,n = G -1 (U i,n ) = (1 -U i,n ) -1 , i = 1, ..., n,
where G -1 denotes the quantile function pertaining to cdf G given by formula [START_REF] Dekkers | On the estimation of extreme-value index and large quantiles estimation[END_REF]. Then, this allows us to write

Y n-i+1,n = (1 -V n (1 -s)) -1 , for i -1 n < s ≤ i n , i = 1, ..., n.
Making use of the previous representation of the order statistics Y n-i+1,n , we may rewrite the three statistics in ( 27) and ( 28) into

S 1 = n k k/n 0 log 1 -V n (1 -s) 1 -U n-k,n -1 ds, S 2 = n k k/n 0 1 -V n (1 -s) 1 -U n-k,n -1+λ
ds, and

S 3 = n k k/n 0 1 -V n (1 -s) 1 -U n-k,n -1+λ log 1 -V n (1 -s) 1 -U n-k,n -1
ds.

Next, we show that, as

n → ∞ √ k (S 1 -1) = W 1n + o p (1) , √ k S 2 -λ -1 = W 2n + o p (1)
,

and √ k S 3 -λ -2 = W 3n + o p (1) ,
where W 1n , W 2n and W 3n are the Gaussian r.v.'s defined in Theorem 1. We will only consider the asymptotic distribution of S 3 . The proofs for S 1 and S 2 use similar arguments. By letting f (x) = x λ-1 log x, the statistic S 3 becomes

S 3 = -(n/k) k/n 0 f 1 -V n (1 -s) 1 -U n-k,n
ds.

An application of standard calculus gives

1 0 f (s) ds = -λ -2 . Therefore S 3 -λ -2 = -(n/k) k/n 0 f 1 -V n (1 -s) 1 -U n-k,n -f s k/n ds.
We divide the integral above in two parts, then we study the asymptotic behavior of each integral. Observe that

S 3 -λ -2 = -(n/k) 1/n 0 f 1 -V n (1 -s) 1 -U n-k,n -f s k/n ds -(n/k) k/n 1/n f 1 -V n (1 -s) 1 -U n-k,n -f s k/n ds =: -∆ n -Ω n .
Next, we show that √ k∆ n converges to 0 in probability. Indeed, we have 1-

V n (1 -s) = 1-U n,n , for 0 < s ≤ 1/n, it follows that ∆ n = (n/k) 1/n 0 f 1 -U n,n 1 -U n-k,n -f s k/n ds = k -1 f 1 -U n,n 1 -U n-k,n - 1/k 0 f (s) ds.
An elementary calculation gives

1/k 0 f (s) ds = λ -1 k -λ log k -1 -λ -1 , and we have n (1 -U n,n-k ) /k P → 1, as n → ∞, therefore ∆ n = {1 + o p (1)} k -λ log k -1 -λ -1 k -λ log k -1 -λ -1 .
Since λ > 1 and k → ∞, then k -λ+1/2 → 0 and k -λ+1/2 log k -1 → 0, it follows that √ k∆ n P → 0

as n → ∞. Consider now the second term Ω n which may be rewritten into

Ω n = (n/k) k/n 1/n f 1 -V n (1 -s) 1 -U n-k,n -f s 1 -U n-k,n ds
+ (n/k) k/n 1/n f s 1 -U n-k,n -f s k/n ds =: Ω n1 + Ω n2 .
Making use of Taylor's expansion of f, we get

f 1 -V n (1 -s) 1 -U n-k,n -f s 1 -U n-k,n = f ′ ϕ n (s) 1 -U n-k,n 1 -V n (1 -s) 1 -U n-k,n - s 1 -U n-k,n ,
and f s 1 -U n-k,n -f s k/n = f ′ (sψ n ) s 1 -U n-k,n - s k/n , where min {1 -V n (1 -s) , s} < ϕ n (s) < max {1 -V n (1 -s) , s} (31)
and

min 1 1 -U n-k,n , 1 k/n < ψ n < max 1 1 -U n-k,n , 1 k/n . ( 32
)
Observe now that Ω n1 and Ω n2 may be rewritten into

Ω n1 = (n/k) k/n 1/n f ′ s 1 -U n-k,n 1 -V n (1 -s) 1 -U n-k,n - s 1 -U n-k,n ds + Ω * n1 ,
and

Ω n2 = {1 + o p (1)} (n/k) k/n 1/n f ′ s k/n s 1 -U n-k,n - s k/n ds + Ω * n2 ,
where

Ω * n1 := (n/k) k/n 1/n f ′ ϕ n (s) 1 -U n-k,n -f ′ s 1 -U n-k,n × 1 -V n (1 -s) 1 -U n-k,n - s 1 -U n-k,n ds,
and

Ω * n2 := {1 + o p (1)} (n/k) k/n 1/n f ′ sψ n k/n -f ′ s k/n × s 1 -U n-k,n - s k/n ds.
From Lemma 3 (see the Appendix), both √ kΩ * n1 and √ kΩ * n2 converge to 0 in probability. Since

n (1 -U n-k,n ) /k P → 1, then Ω n1 = {1 + o p (1)} k/n 1/n f ′ s k/n [1 -s -V n (1 -s)] ds + o p k -1/2 ,
and

Ω n2 = -{1 + o p (1)} k/n -(1 -U n-k,n ) (k/n) 2 k/n 1/n s k/n f ′ s k/n ds + o p k -1/2 .
The derivative of function

f equals f ′ (x) = (λ -1)x λ-2 log x + x λ-2 , then Ω n1 = (λ -1) (n/k) 1 1/k t λ-2 (log t) [1 -V n (1 -kt/n) -kt/n] dt + (n/k) 1 1/k t λ-2 [1 -V n (1 -kt/n) -kt/n] dt + o p k -1/2 ,
and

Ω n2 = (λ -1) (n/k) [k/n -(1 -U n-k,n )] 1 1/k t λ-1 log tdt + (n/k) [k/n -(1 -U n-k,n )] 1 1/k t λ-1 dt + o p k -1/2 = λ -2 (n/k) [k/n -(1 -U n-k,n )] + o p k -1/2 .
Fix 0 < τ < 1/2, then using approximation [START_REF] Delaigle | Estimation of boundary and discontinuity points in deconvolution problems[END_REF], in Lemma 1, yields

√ kΩ n1 = (λ -1) n/k 1 1/k t λ-2 (log t) B n (1 -kt/n) dt + n/k 1 1/k t λ-2 B n (1 -kt/n) dt + √ k Ω n1 (τ ) + o p (1) , and
√ kΩ n2 = -λ -2 n/kB n (1 -k/n) + √ k Ω n2 (τ ) + o p (1) , where √ k Ω n1 (τ ) = (λ -1) O p n -τ (k/n) 1/2-τ (n/k) 1/2 1 0 t λ-2+(1/2-τ ) |log t| dt + O p n -τ (k/n) 1/2-τ n/k 1 0 t λ-2+(1/2-τ ) dt, and
√ k Ω n2 (τ ) = λ -2 O p n -τ n/k (k/n) 1/2-τ .
For λ > 1,

1 0 t λ-2+(1/2-τ ) |log t| dt = (λ -1/2 -τ ) -2 and 1 0 t λ-2+(1/2-τ ) dt = (λ -1/2 -τ ) -1
are finite integrals. Then both quantities √ k Ω n1 and √ k Ω n2 are equal to O p (k -τ ) for all large n, which tends in probability to 0 as n → ∞. Recall that up to now we have showed that

√ kΩ n1 = (λ -1) n/k 1 1/k t λ-2 (log t) B n (1 -kt/n) dt + n/k 1 1/k t λ-2 B n (1 -kt/n) dt + o p (1)
,

and √ kΩ n2 = λ -2 n/kB n (1 -k/n) + o p (1)
.

It remains to prove that

I n := (λ -1) n/k 1/k 0 t λ-2 (log t) B n (1 -kt/n) dt + n/k 1/k 0 t λ-2 B n (1 -kt/n) dt,
converges, in probability, to 0. Indeed, since

E |B n (1 -ks/n)| ≤ ks/n, then E |I n | ≤ (λ -1) 1/k 0 t λ-2+1/2 (|log t| + 1) dt. Since 1/k 0 t λ-2+1/2 (|log t| + 1) dt = 2 (2λ -1) 2 k -λ+ 1 2 (2λ -log k + 2λ log k + 1) ,
which tends to 0 as n → ∞, then I n converges to 0 in probability. This completes the proof of Theorem 1.

Proof of Theorem 2.

To establish the asymptotic normality of α, given in [START_REF] Daouia | A γ-moment approach to monotonic boundaries estimation: with applications in econometric and nuclear fields[END_REF] , we proceed by similar arguments as for µ n in the proof of Theorem 4.

Proof of Theorem 3

Let us divide the integral (2), in two parts, as follows:

µ = µ 1,n (k) + µ 2,n (k) ,
where

µ 1,n (k) := k/n 0 Q (1 -s) ds and µ 2,n (k) := 1 k/n Q (1 -s) ds.
Recall that, in Section 1 formula ( 16), we have defined estimator µ n of µ by

µ n = µ 1,n (k) + µ 2,n (k) ,
where

µ 1,n (k) := (k/n) (n c/k) 1/ α α α -1 + d c -β/ α (k/n) β/ α-1 β -1 and µ 2,n (k) := 1 n n i=k+1
X n-i+1,n .

To simplify notations, let us set

Z ni := √ n k/n (nc/k) 1/α { µ i,n (k) -µ i,n (k)} , i = 1, 2. (33)
First, we consider Z n1 . It is easy to verify that, as n → ∞

µ 1,n (k) = {1 + o p (1)} k n (nc/k) 1/α α α -1 ,
and, under the condition (6), we have

µ 1,n (k) = {1 + o p (1)} k n α α -1 (n c/k) 1/ α .
It follows that

µ 1,n (k) -µ 1,n (k) = {1 + o p (1)} k n α α -1 (n c/k) 1/ α - α α -1 (nc/k) 1/α .
Let us write Z n1 = T 1n + T 2n , where

T 1n := {1 + o p (1)} √ k α α -1 - α α -1 , and
T 2n := {1 + o p (1)} √ k (n c/k) 1/ α (nc/k) 1/α -1 .
We begin by showing that T 2n P → 0, as n → ∞. First observe that T 2n may be rewritten into

T 2n = {1 + o p (1)} √ k (nc/k) 1/ α-1/α ( c/c) 1/ α-1/α -1 .
Assumptions (i) and (ii) of Theorem 1 imply that k 1/2 / log (n/k) → ∞. Also, the asymptotic normality of α gives αα = O p k -1/2 . Therefore (1/ α -1/α) log (nc/k) P → 0, this implies that (nc/k) 1/ α-1/α P → 1, as n → ∞. On the other hand, we have

c/c -1 = α -1 ( α -α) log n k + o p k -1/2 log n k .
Since c is a consistent estimator of c, then Taylor's expansion gives

( c/c) 1/ α-1/α -1 = α -1 (1 + o p (1)) ( α -α) ( c/c -1) , as n → ∞.
It suffices now to show that √ k ( c/c) 1/ α-1/α -1 converges to 0 in probability. Indeed, again by using the fact that α -

α = O p k -1/2 , yields √ k ( c/c) 1/ α-1/α -1 = O p (1)   k -1/2 log n k + o p   log n k √ k     ,
which tends in probability to 0, because we already have √ k/ log (n/k) → ∞. Now, we consider the term T 1n . Since α is a consistent estimator of α, then it is easy to show that

T 1n = - 1 + o p (1) (α -1) 2 √ k ( α -α) , as n → ∞.
From 1, we infer that

T 1n = - (1 + o p (1)) α (α -1) 2 (η 1 W 1n + η 2 W 2n + η 3 W 3n ) , as n → ∞.

It follows that

Z n1 = - α (α -1) 2 {η 1 W 1n + η 2 W 2n + η 3 W 3n } + o p (1) . ( 34
)
Let us now consider the asymptotic distribution of Z n2 , in [START_REF] Gardes | A moving window approach for nonparametric estimation of the conditional tail index[END_REF]. It is shown that

Z n2 = - 1 k/n B n (1 -s) dQ (1 -s) k/nQ (1 -k/n) + o p (1)
.

On the other hand, from [START_REF] Byström | Managing extreme risks in tranquil and volatile markets using conditional extreme value theory[END_REF], we have

Q (1 -k/n) ∼ (nc/k) 1/α , as n → ∞, it follows that Z n2 = W 4n + o p (1) . (35)
Combining ( 34) and ( 35) achieves the proof of Theorem 3.

Proof of Theorem 4.

Now, we investigate the asymptotic normality of µ n given in [START_REF] Daouia | Robust nonparametric frontier estimators: qualitative robustness and influence function[END_REF]. Since W in , i = 1, ..., 4 are sequences of centred Gaussian r.v.'s, then

√ n k/n (nc/k) 1/α { µ n -µ} d → N 0, ΓΣΓ t , as n → ∞, where Γ := - α (α -1) 2 η 1 , - α (α -1) 2 η 2 , - α (α -1) 2 η 3 , 1 ,
Γ t is the transpose of Γ and Σ is the variance-covariance matrix of the vector (W 1n , ..., W 4n ) defined by

=       1 1 λ 2 -1 λ 2 λ 3 -1 λ 2 -1 1 λ 2 -1 λ 1 2λ-1 -1 λ 2 1 (2λ-1) 2 -1 λ 3 λ-1 λ 2 λ 3 -1 λ 2 1 (2λ-1) 2 -1 λ 3 2 (2λ-1) 3 -1 λ 4 -1 λ 2 -1 λ-1 λ -1 λ 2 2 2-α       .
Note that the elements of were obtained after tedious computations of the limits of the expectations E [W in W jn ] for i, j = 1, 4 (i ≤ j). Finally, a standard calculation of the product ΓΣΓ t yields

ΓΣΓ t = α 2 β 4 (α -1) 4 (α -β) 4 + 2 2 -α + 2αβ 2 (α -1) 2 (α -β) 2 ,
which is denoted by σ 2 (α, β) . This completes the proof of Theorem 4.

Conclusion

The main objective of this paper was to propose a bias-reduced estimator for the expectation of a heavy-tail distribution. This was achieved on the basis of the bias-reduction of the first and second order parameter estimators of regularly varying distributions and the corresponding high quantiles estimators. In addition, the newly introduced estimator is asymptotically normal, making confidence intervals easily constructible.

A Appendix: Auxiliary results

Lemma 2 Let k = k n be a sequence of integers satisfying (6) and f

(x) = x λ-1 log x, λ > 1. Then, uniformly on s ∈ [1/n, k/n] , we have f ′ ϕ n (s) 1 -U n-k,n -f ′ s 1 -U n-k,n = o p (1) s k/n λ-2 log s k/n , as n → ∞. Proof. We have n (1 -U n-k,n ) /k p → 0, as n → ∞, then f ′ ϕ n (s) 1 -U n-k,n -f ′ s 1 -U n-k,n = (1 + o p (1)) f ′ ϕ n (s) k/n -f ′ s k/n .
A straightforward calculation of the derivative of f yields

f ′ ϕ n (s) k/n -f ′ s k/n = s k/n λ-2 (λ -1) ϕ n (s) s λ-2 -1 log ϕ n (s) k/n + (λ -1) log ϕ n (s) s + ϕ n s λ-2 -1 . (36)
Observe now, that inequalities (31) imply

min 1 --V n (1 -s) s , 0 < ϕ n (s) s -1 < max 1 -s -V n (1 -s) s , 0 .
We have sup

On the other hand, we infer that sup

1/n≤s≤k/n s ϕ n (s) -1 p → 0 as n → ∞. (38)
By applying the mean value theorem to the functions x → log x and x → x λ-1 respectively, then by using [START_REF] Gardes | On the estimation of the functional Weibull tail-coefficient[END_REF] and [START_REF] Gardes | Functional nonparametric estimation of conditional extreme quantiles[END_REF], we show readily that, as n → ∞ Proof. We only show the first result. The second one is obtained by similar arguments. Recall that

Ω * n1 := -(n/k) k/n 1/n f ′ ϕ n (s) 1 -U n-k,n -f ′ s 1 -U n-k,n × 1 -V n (1 -s) 1 -U n-k,n - s 1 -U n-k,n
ds.

Using Lemma 2, we get, as n → ∞

Ω * n1 = o p (

  for convergence in probability, convergence in distribution and equality in distribution, while N a, b 2 denotes the normal distribution with expectation a and variance b 2 .

1 p→

 1 s -V n (1s)| s p → 0 as n → ∞, it follows that sup 1/n≤s≤k/n ϕ n (s) s -0 as n → ∞.

→ 0 ,

 0 as n → ∞.[START_REF] Ghorbel | Predictive performance of conditional extreme value theory in value-at-risk estimation[END_REF] By using equations (39) and (40) together, we show that, uniformly in s ∈ [1/n, k/n] , the right-hand side of equation[START_REF] Gardes | Nonparametric estimation of the conditional tail copula[END_REF] , is equal to o p (1as n → ∞.

2 1 0s λ- 2 1 0 1 0s λ- 3 / 2 -τ |log s| ds and 1 0

 21211321 |log s| |1 -V n (1ks/n) ks/n| ds.Making use of approximation[START_REF] Delaigle | Estimation of boundary and discontinuity points in deconvolution problems[END_REF] yields√ kΩ * n1 = o p (1) (n/k) s λ-2 |log s| |B n (1ks/n)| + (ks/n) 1/2-τ O p n -τ ds.In other words√ kΩ * n1 = o p (1) (n/k) s λ-2 |log s| |B n (1ks/n)| ds. + o p k -τ 1 1/k s λ-3/2-τ |log s| ds.The expectation of the first term of right-hand side of the previous equation is less than or equal too p (1) (n/k) 1/2 |log s| E [|B n (1ks/n)|] ds.Using the fact thatE |B n (1ks/n)| ≤ (ks/n)1/2 , we show that the previous quantity is less than or equal to o p (1) s λ-3/2 |log s| ds. Since both integrals s λ-3/2-τ |log s| ds are finite, then √ kΩ * n1 = o p (1) , as n → ∞.