Maya N Nair

Irene Palacio

Arlensiú Celis

Alberto Zobelli

Alexandre Gloter

Stefan Kubsky

Jean-Philippe Turmaud

Matthew Conrad

Claire Berger

Walter De Heer

Edward H Conrad

Amina Taleb-Ibrahimi

Antonio Tejeda

Band Gap Opening Induced by the Structural Periodicity in Epitaxial Graphene Buffer Layer

The epitaxial graphene buffer layer on the Si face of hexagonal SiC shows a promising band gap, of which the precise origin remains to be understood. In this work, we correlate the electronic to the atomic structure of the buffer layer by combining angle resolved photoemission spectroscopy (ARPES), scanning tunneling microscopy (STM), and high-resolution scanning transmission electron microscopy (HR-STEM). We show that the band structure in the buffer has an electronic periodicity related to the structural periodicity observed in STM images and published X-ray diffraction. Our HR-STEM measurements show the bonding of the buffer layer to the SiC at specific locations separated by 1.5 nm. This is consistent with the quasi 6 × 6 periodic corrugation observed in the STM images. The distance between buffer C and SiC is 1.9 Å in the bonded regions and up to 2.8 Å in the decoupled regions, corresponding to a 0.9 Å corrugation of the buffer layer. The decoupled regions are sp 2 hybridized. Density functional tight binding (DFTB) calculations demonstrate the presence of a gap at the Dirac point everywhere in the buffer layer, even in the decoupled regions where the buffer layer has an atomic structure close to that of graphene. The surface periodicity also promotes band in the superperiodic Brillouin zone edges as seen by photoemission and confirmed by our calculations.

Graphene is a promising candidate for nanoelectronics because of its exceptional electronic properties. [1][2][3][4][5] However, the absence of energy gap in this zero-gap semiconductor limits its applications in digital electronics. This problem has therefore motivated many different approaches to opening a gap ranging from chemical modification by doping or functionalization with atoms, molecules, or a substrate to patterning graphene into nanoribbons. 1,[6][7][8][9][10][11][12][13][14][15][16][17] In scalable nanoribbons, the precise control of the electronic properties requires edges controlled at an atomic level, which is still beyond the accuracy of current lithographic techniques. In fact, the observed transport gap for lithographic ribbons is due to localization effects caused by disordered ribbon edges rather than finite size effects. 16 A particular interesting variety of graphene ribbons, 17 grown on SiC facet templates, have their position defined by lithography prior to growth. Their ballistic transport 1 indicates smooth edges, and the well-defined mini-ribbons that border the facets show a band gap. 18,19 The problem of tailoring the band gap in ribbons by controlling their width to nanometer dimensions could be avoided if a bottom-up graphene structure could be grown with an intrinsic band gap. Such a variety indeed exists on silicon carbide and is called the "buffer" layer. When SiC(0001) is annealed at temperatures below the optimal temperature for growing graphene, an atomic layer of carbon atoms with honeycomb structure appears. 20,21 This buffer layer is semi-conducting with electronic properties different from those of graphene as shown in Figure 1a. Upon further annealing, another graphene layer forms below the existing buffer layer so that the buffer converts to a conducting graphene layer with a new buffer layer under it at the SiC interface. The top graphene layer exhibits the ideal dispersion of graphene [22][23][24] (Figure 1b) and is n-doped with the Dirac point -0.4 eV below the Fermi level (E F ). The position of the Dirac point has been proposed to be caused by pinning E F by interface states. [25][26][27] This situation shifts graphene bands in energy with respect the buffer layer. The buffer layer preserves the π band of graphene although it becomes gapped by more than 0.5 eV, 28 twice the previous value in epitaxial graphene samples. 29 The buffer layer has been extensively studied. The atomic structure has been studied by STM, 21,[35][36][37][38]40,43 LEED, 22,38,40,41 and other diffraction techniques 30,43 and the electronic structure with ARPES 22,23,28,29 and by STEM-EELS spectroscopy. 19 Many calculations have also been performed to understand its structure and its electronic properties. 20,26,27,[31][32][33][34] The buffer layer is observed as a (6 × 6) reconstruction in STM, [35][36][37][38]40 while LEED pattern indicates a (6√3 × 6√3)R30°2 2,35,[39][40][41] structure, but the apparent contradiction can be explained. 27 From the electronic point of view, the buffer layer exhibits a band gap, as observed by ARPES 22,28 or by STS. 42 This band gap could in principle appear due to a complete covalent bonding of the buffer layer with the substrate, suppressing the π band of graphene and rendering the system insulating. However, such a global covalent bonding is hard to reconcile with the fact that hydrogen intercalation lifts the buffer from the substrate, decoupling it and converting it in quasi-free-standing monolayer graphene. 24 The current explanation for buffer layer band gap, while preserving the honeycomb lattice, is to attribute the gap to the formation of covalent C-Si bonds between the buffer and the substrate on particular sites of the honeycomb lattice that are in near coincidence with Si atoms in the substrate surface. 21,22,34 The gap's origin has also been suggested to be the result of a pseudo nanoribbon network of regions decoupled from the SiC. 28 Given the potential technological interest of having a gapped material entirely made of graphene, understanding the origin of the band gap in the buffer layer is important if the gap is to be used and controlled in relevant applications. Here we show that the electronic structure of the buffer layer is strongly modified by the observed 6 × 6 periodicity. Structural analysis is performed to quantify the presence of defects on the surface, the bonding to the substrate, and the correlation between the buffer's electronic and structural periodicities. For this we combine STM, cross-sectional-scanning transmission electron microscopy (X-STEM), and angle resolved photoemission spectroscopy (ARPES) together with density functional based tight binding (DFTB) calculations to give important insights into the nature of the buffer's observed band gap.

A well-suited technique to explore the formation of band gaps is angle-resolved photoemission that allows to map the band structure over the whole reciprocal space. Figure 2a shows a schematic of reciprocal space indicating the graphene Brillouin zone and its high symmetry points as well as a set of the smaller Brillouin zones of a (6 × 6) SiC incommensurate period recently observed by X-ray diffraction. 43 The incom-mensurate density modulation has a fundamental wavelength of 6(1 + δ)a SiC (where a SiC is bulk SiC lattice constant and δ = 0.037). The band structure of the buffer layer measured along the ΓKM direction (K y ) of graphene (red arrow in Figure 2a) is shown in Figure 2b. Figure 2c shows an ARPES cut through the K point in a direction perpendicular to ΓK, K x (blue arrow in Figure 2a). The spectrum shows σ bands in the buffer layer as expected from the presence of in-plane σ bonds in the buffer honeycomb structure (Figure 2b). The spectra show no evidence of a metallic linear dispersing π band either in the ΓK direction or in the perpendicular direction through the K point. A close up of the bands at the K point shows the gapped spectral feature ε 1 (Figure 2d) in agreement with ref 28. Its spectral weight is distributed along a parabola, with a wide gap at the Dirac point of more than 0.4-0.5 eV in the occupied part of the band structure (estimated from the energy of the first electronic state to E F ), confirming the semiconducting properties of the buffer layer. The highest spectral weight in ε 1 is reached at the energies around -0.4 eV and -1.4 eV, which correspond roughly to the binding energies of the g 1 and Figure 1. Electronic band structure through the K point in a direction perpendicular to ΓK (see Figure 2a for the Brillouin zone scheme) for (a) buffer layer and (b) graphene layer above the buffer layer. A band gap opening at the K point is clearly present in the buffer layer as its characteristic electronic states do not reach the Fermi level (panel a). The overall band structure is shifted toward higher binding energies when metallic graphene is grown on top of the semiconducting buffer layer (panel b). The shift can be interpreted as being due to a Schottky barrier between a metal and a p-doped semiconductor. Equivalent states in the buffer and in graphene system are indicated by dashed lines to highlight the shift in the electronic states upon graphene growth on top of the buffer layer. The graphene dispersion is marked in panel b by dashed lines. g 2 states clearly observed in the earlier stages of growth and that have been previously observed in UHV-grown buffer layer. 22 The second derivative of Figure 2b-d is shown, respectively, in Figure 2e-g. This is the usual way to analyze photoemission data to enhance the electronic dispersions. The spectral weight around -0.4 eV of the ε 1 feature has a very small dispersion, but the spectral weight feature around -1.4 eV disperses significantly and seems to reach maxima and minima at reciprocal space positions close to a (6 × 6) SiC superperiodicity. This is demonstrated in Figure 2g where the vertical dashed lines, marking the (6 × 6) SiC zone boundaries perpendicular to the ΓK direction, line up with maxima in the ε 1 band intensities. The periodicity of ε 1 is further observed in the E(k) dispersions along the ΓM 6×6SiC directions shown in the Supporting Information (Figure S.2).

The (6 × 6) SiC periodicity of ε 1 is also seen in the constant energy cuts of the electronic structure. Figure 3 shows constant energy cuts at four different energies. The constant energy cuts in Figure 3b and c are near the graphene K point at energies, where ε 1 has stronger spectral weight. The spectral weight in these constant energy cuts is not constant in reciprocal space, indicating that the state observed there disperses. Figure 3b andc show that the spectral weight of the p z -character of ε 1 features is indeed periodic, with a spatial period much larger than the lattice parameter of graphene. This p z -derived state is sensitive to the structural superperiodicity of the (6 × 6) SiC reconstruction. Regardless of the precise spatial periodicity (as different commensurate and incommensurate lattice parameters may be defined), if a π-derived band was sensitive to a strong superperiodic potential, a gap could appear at the Brillouin zone edges of the superperiodicity due to Bragg reflection on them. Such a situation could explain the fact that the spectral weight of the ε 1 feature is discontinuous in energy. We first try to figure out the periodicity of this p z -derived state only visible close to the K point of graphene. ε 1 is different from previously observed g1 and g2 states of the buffer layer 22 as these states are visible in a larger region of the reciprocal space with a different symmetry (see Supporting Information (Figure S.1). ε 1 is compatible with the orientation of SiC Brillouin zones, with a size close to that of the (6 × 6) SiC ,as seen in Figure 3b andc. High spectral weight is not expected in every Brillouin zone as the photoemission intensity is not necessarily periodic, [44][45][46] as particularly evident in the paradigmatic case of the (7 × 7) reconstruction of Si(111), where the spectral weight of the surface state does not seem to follow the (7 × 7) periodicity. 44 This is due to the influence of the structure factor in photoemission, analogous to the situation in X-ray diffraction. The (6 × 6) SiC is thus a plausible periodicity according to the constant energy cut, whose atomic origin is to be identified. It is also important to note that the incommensurate modulation of the buffer electronic density has been shown to be modulated by reciprocal lattice vectors of SiC. 43 The periodicity of the buffer's σ band is different from the p z bands. Figure 3d ande show constant energy cuts at -2.3 eV and -3.45 eV near Γ (Figure 3a) where σ band exists. Figure 3d and e show that the highly dispersive σ band is periodic. Any true periodicity must be compatible with the constant energy cuts shown in Figure 3d and e, as they are intercepting the same electronic state, as observed in Figure 3a. It can be observed that the highly dispersive σ band follow bigger Brillouin zones than the π-derived band, oriented as those of SiC and close to the expected ones for (√3 × √3)R30°S iC . A similar periodicity has been also observed by other groups. 37,38 Indications on the atomic origin of the periodicities can be obtained by combining STM and X-ray diffraction. Figure 4a shows the characteristic STM topographic image of the buffer layer. It shows that a honeycomb graphene layer is modulated by an apparent (6 × 6) SiC reconstruction (about 1.8 nm periodicity) that was previously reported for the buffer layer 38,39,[47][48][49][50] and compatible with the incommensurate modulation with a period close to the (6 × 6) SiC cell reconstruction observed by X-ray diffraction. 43 The structural interpretation of the periodicity felt by the π-derived states is found in Figure 4a, where the (6 × 6) SiC unit cell indicates a pseudoperiodicity in STM images. STM images also allow to discard a possible band gap opening in the buffer layer due to defects. Figure 4b shows a very regular honeycomb lattice with no evidence for hexagon-pentagon-heptagon defects that were suggested to cause a band gap. 51 We can also rule out that the buffer layer is uniformly strongly coupled to the SiC substrate, which would completely suppress the π-band. Information on the coupling can indeed be provided by cross-sectional STEM images, which are sensitive essentially to the atomic structure, and not convoluted with electronic structure effects like STM. Figure 4c shows a STEM cross sectional view of the buffer-SiC interface along the crystallo-graphic zone axis [110] SiC . The image was obtained with a 60 keV electron beam to avoid sample damage. Figure 4c corresponds to a high angle annular dark field (HAADF) collection mode that is sensitive to the heaviest elements. Bright spots are therefore Si atoms on the substrate and C atoms in the buffer layer. The crosssectional image conclusively shows that the buffer layer is pinned to the substrate with an apparent periodicity of about 1.5 nm, which corresponds to the projection along the [1-10] SiC of the (6 × 6) SiC lattice parameter. The 1.8 nm period is in perfect agreement with of the (6 × 6) SiC lattice parameter found in STM.

In the place where the buffer carbon is pinned to the SiC, the buffer-substrate distance is reduced to ∼1.9 Å (Figure 4d). The pinning distance is very close to the Si-C bond length in SiC (1.89 Å) implying a strong sp 3 covalent bonding of the buffer, pinning the carbon to Si in the SiC interface as suggested by Emtsev et al. 22 In the decoupled regions, the distance of the buffer layer to the substrate is much larger, ∼2.8 Å. This distance is only 15% less than the distance between van der Waals bonded layers in graphite. The result is that the buffer layer has a total vertical corrugation close to 1 Å. Additional structural information on the buffer-SiC interface can be obtained using a different STEM collection mode, that is, annular bright field (ABF), which is sensitive to light elements like carbon. Figure 4e shows the result of an ABF measurement. Both Si and C atoms are clearly resolved on the substrate, and since the atomic resolution is very high, it is possible to extract information on the subsurface structural modification when the buffer layer is grown on SiC. To quantify deviations with respect to the bulk structure, we have extracted the intensity profile from Figure 4e and plotted it in Figure 4f (the procedure is described in the Methods section). In this way, we observe a clear separation of the Si and C atoms in the SiC bilayer located at several nanometer from the interface. The distances between the Si and C layers are measured at about 60 pm that is in fair agreement with the expected 63 pm along [001] SiC for a bulk structure. Deviations from the bulk structure become visible for the last three SiC bilayers and are substantial for the interface SiC bilayer just below the buffer layer, where the peak separation is not any more visible.

In addition to the structural observations, numerical simulations support both the electronic and structural data described above. We have performed density functional tight binding approach (DFTB) calculations, which combines high accuracy from DFT based parameters with the reduced computational cost of standard tight binding (see the Methods section for more details). This technique is thus advantageous for describing large unit cells with complex bonding reconfiguration such as the (6√3 × 6√3)R30°S iC buffersubstrate interface (that corresponds to a 13 × 13 graphene supercell). The resulting optimized atomic positions from the calculations are shown on the inset of Figure 4c. A histogram of the calculated z atomic coordinates of the atoms is shown in Figure 4g (the z-coordinates of all the atoms are summed over the unit cell). These calculations show that C and Si layers exhibit a certain disorder close to the interface, especially at the last SiC bilayer, explaining the change of contrast in the experimental profile (Figure 4f). The calculated buffer layer presents a corrugation, with a periodicity and an intensity comparable to the STEM-images of Figure 4c, e, andf. Besides this structural information, the atomic configuration of atoms in the buffer layer indicates that atoms in the decoupled regions are sp 2 hybridized, indicating that the buffer layer has the π-related band of graphene.

The electronic structure derived from this extended supercell model results in highly folded bands. To understand how the electronic structure of the buffer layer is modified with respect to that of graphene, we employed an unfolding postprocessing which provides a primitive cell effective band structure (see the Methods section for more details). Band unfolding for epitaxial graphene over SiC has been already presented in the work of Kim and co-workers using a DFT plane waves approach. 27 Here we use a localized basis set, suitable for larger clusters and providing an easier interpretation of the states by projecting on a subset of orbitals or atoms. We considered the spectral 2a). Color scale represents the spectral weight. The graphene band structure has been superimposed for reference by aligning the deepest σ band. Vertical dotted lines in a represent the center of several Brillouin zone of the (6 × 6) SiC superperiodicity. An energy broadening of 10 meV has been applied to the bands. Left energy scales refer to the Fermi level as obtained in the supercell model, the right scales to the Fermi level in monolayer graphene. A close-up of the electronic structure near the Fermi level is shown in c for the ΓKM direction and in d for the perpendicular direction through K. The density of states for the buffer layer and for graphene is shown in panels e and f, corresponding to the same energy windows used in panels b and d, respectively. contribution for sole carbon atoms at the buffer layer, which represents a good approximation for the limited penetration of photons in ARPES. Figure 5 shows the band structure of the (6√3 × 6√3)R30°S iC reconstruction unfolded over the graphene Brillouin zone as well as the density of states. The perfect graphene band structure has been superimposed by aligning the deepest σ band (white continuous line). The left energy scale refers to the Fermi level as obtained in the supercell model for the buffer layer, the right scale to the Fermi level in perfect graphene. Note that SiC has an intrinsic polarization, which, in the absence of dipole corrections, gets equilibrated in supercell calculations by a surface charge overestimation. This is due to the periodic boundary condition that is imposed to the potential in the supercell. This effect leads to a spurious upper shift of the Fermi level at the buffer layer. To compare with experimental data it is therefore more pertinent to refer all calculated energies from the Dirac point energy of the perfect graphene.

In agreement with previous results, we observe that σ bands in the buffer layer and in graphene disperse similarly. The translation symmetry breaking in the buffer supercell is responsible for a residual spectral weight at the unfolded σ bands, which is further modulated by the apparent (6 × 6) SiC reconstruction. This effect is clearly visible in Figure 5a, where the center of several Brillouin zones of the (6 × 6) SiC reconstruction has been indicated with vertical dotted lines.

The π band of the buffer layer can also be identified, whose dispersion is significantly different to that of graphene. The π band splits at the Γ point due to rehybridization with SiC bands (for SiC allotropes band structures and surface states see for instance ref 52). For higher wave numbers, we observe a maximum spectral weight at K where a "X"-shaped feature arises about 1.2 eV below the graphene Dirac cone. A similar remarkably intense "X"-shaped feature is also experimentally observed both in the buffer and in the graphene monolayer (see Figure 1).

Above the Dirac cone position, the system presents a calculated 0.4 eV band gap, clearly visible as a dark line in the band structure but also in the density of states. Note that DFTB based calculations usually underestimate the band gap values and 0.4 eV might be seen as a lower bound. As mentioned previously, the position of the Fermi level in the calculation is overestimated (it is located 1 eV above the Dirac point energy of the graphene, it is in the conduction band of the buffer layer). Nevertheless, the real position of the Fermi level can be estimated with the help of experimental observations. It has been reported that graphene on top of the buffer layer is slightly n doped leading to a Fermi level several hundreds of meV above the graphene Dirac cone. This provides a constraint on the buffer layer Fermi level position which therefore lies within the computed band gap indicating the semiconducting character of the buffer layer.

The weakly dispersive bands observed experimentally between the intense "X"-shaped feature and the band gap also appear in these calculations. In the ARPES measurement, these states have maxima roughly distributed along a parabola (see Figure 2d,g) and have a lower velocity than the dispersion of isolated graphene would predict. The calculated bands (Figure 5c,d and Figure S.4) have more complicated and sharper details, notably due to lack of disorder, incommensurability, and other broadening source that are not included in the calculations. Nevertheless, the weakly dispersive character and the appearance of discontinuous contrast in the spectral weight are also observed in the calculation. The comparison of the DOS calculated for the graphene and the buffer layer is also of interest. In the occupied domain, both DOS match rather well, except in the range from 0 to -1 eV. In this energy domain, narrow peaks are only present for the buffer layer. Their energies are very comparable with the experimentally reported maxima energies of -0.4 and -1.4 eV (Figure 2d), notably keeping in mind that the Fermi level in the calculated DOS lie several hundred of meV above the zero. Such electronic states may be related to local orders in the (6√3 × 6√3)R30°supercell and possibly with the apparent 6 × 6 reconstruction. Such structural periodicities are known to be present as measured by several different structural techniques STM, TEM, LEED, or XRD (Figure 4), corresponding to a strong interaction between the buffer layer and the underlying substrate, evident from the formation of out-of-plane covalent bonds. We thus conclude that not only the σ bands of graphene are affected by the structural periodicities of SiC but also the π bands, which are at the origin the band gap opening.

In conclusion, we have analyzed the band gap opening in the graphene buffer layer by employing a complementary set of experimental and theoretical approaches that include STM, STEM, and ARPES with DFTB based calculations. The combined results from this approach allow us to conclude that the gap around the Dirac point does not originate from a significant density of hexagon-pentagon-heptagon defects or the complete pinning of the buffer layer to the substrate. As suggested by STM images, cross-sectional STEM demonstrates that the buffer layer is only pinned at some precise locations. In between these pinning regions, separated by 1.5 nm, the buffer shows areas of decoupled graphene, with sp 2 hybridization. Nonetheless, these decoupled areas exhibit a band gap, as the whole buffer layer system. The main result is that the buffer layer π band shows band gaps aligned with the substrate 6 × 6 nanometric superperiodicity. The induced band gap in this periodic graphene system resembles the gap opening in graphene nanoribbons. 19 However, the difference between these two systems is that in graphene nanoribbons, the periodicity plays no role in the band gap opening. In the buffer layer, the situation is completely different. The band gap opens in the π band due to a nanoscale superperiodicity. These results indicate new ways of controlling graphene's band structure simply by transferring or growing graphene on top of self-organized or nanostructured systems.

METHODS. Experimental Details. Graphene samples were prepared from a n-doped 4H-SiC(0001) substrate in a controlled silicon sublimation furnace. 53 The semiconducting buffer layer graphene in this letter was prepared using the same procedure as in Nevius et al. 28 Graphitized samples were transported in air and thermally annealed at 500 K in UHV prior to measurements. ARPES measurements were carried out at the CASSIOPEE beamline at SOLEIL synchrotron radiation source. The beamline is well equipped with a Peterson PGM monochromator with an energy resolution E/ΔE of 25 000 for energies below 100 eV. The detector is a Scienta R4000 detector with a base resolution of ΔE < 1 meV. Binding energies are referenced to the Fermi level of a metal in electrical contact with the sample. The photon energy used for the constant energy cuts in Figure 3 was 70 eV. The photon energy for Figure 1a and b was 70 and 150 eV, respectively. The STM experiments were performed with a variable temperature (VTSTM) from Omicron. Bias voltages were applied to the grounded PtIr tip. STEM-EELS measurements were carried out using a STEM NION microscope and an in-house modified GATAN spectrometer. To limit the electron beam damage, most of the experiments were done with 60 keV electron beam, unless otherwise specified. Atomic resolution was obtained in STEM mode by correcting the 3/5 order spherical aberration of the objective lens. Figure 4c is a high angle annular dark field (HAADF) image made by collecting scattered electrons with angles between 80 and 180 mrads. Figure 4e is an STEM-ABF (annular bright field) measurement performed at 200 keV with a resolution close to 70 pm. The intensity profile from this figure has been obtained by integrating over 5 nm in a direction parallel to the surface after applying a deconvolution filter to sharpen the atomic columns (such integration increases the signal-to-noise ratio and it is enough to average over several reconstruction periods). To compare with the STEM-HAADF profile (Figure 4d), the STEM-ABF profile contrast has been reversed (Figure 4f). Sample preparation for STEM is described elsewhere. 19 Computational Details. We consider a commensurate (6√3 × 6√3)R30°SiC supercell where the silicon terminated upper SiC surface matches the carbon buffer layer's area (corresponding to a graphene 13 × 13 supercell). To avoid spurious interactions from the bottom surface, we use a 12 atomic layer model of the SiC substrate (one SiC-6H unit cell along the c direction) with dangling bonds at the lower carbon surface saturated by hydrogen atoms. The total size of the system was 1742 atoms. All calculations have been performed using the DFTB scheme in its self-consistent extension, 54 as implemented in the dftb+ code. 55 We employed the pbc-0-3 parameter set that has been demonstrated to provide a good representation of structural and electronic properties for SiC in respect to DFT-LDA 56 (band dispersions are only slightly underestimated). Atomic positions have been relaxed until forces were below 10 -4 eV/Å on each atom. To compare the buffer layer electronic band structure with respect to the perfect graphene bands, we employed the unfolding postprocessing method developed by Lee and coworkers. 57,58 This technique, originally derived for Wannier and linear combination of atomic orbitals basis functions, can be directly employed in the framework of DFTB. 59 In the representation of Figure 5,an energy broadening of 10 meV has been applied to the bands.

Figure 2 .

 2 Figure 2. Electronic band structure of the buffer layer, with the binding energy relative to E F = 0. (a) Scheme of the Brillouin zones of the system. The large black hexagon corresponds to graphene Brillouin zone. The small green hexagons correspond to a pseudo (6 × 6) SiC . The red(blue) arrow corresponds to the ARPES scan directions in b(c). (b) ARPES measurement along ΓK direction showing the band gap at the K point and the presence of the σ band characteristic of in-plane bonds in graphene-like structures. (c) ARPES measurement through the K point in a direction perpendicular to ΓKM. (d) Zoom of c (from the dotted rectangle in c) showing the detail of the electronic structure of the buffer layer. A faint Dirac cone at a domain with graphene is also present, as indicated by the dashed lines crossing themselves at the Dirac point. A dotted parabola marks the spectral weight of the ε1 spectral feature of the buffer layer. It has intensity maxima at approximately -0.4 and -1.4 eV. (e-g) Second derivatives of b-d, respectively. A low intensity linear dispersion is also observed due to the presence of small domains of graphene. The high symmetry points of the pseudo (6 × 6) SiC periodicity are indicated on the top of the panel g.

Figure 3 .

 3 Figure 3. Constant energy cuts through the electronic structure the buffer layer. (a) E(k) dispersion relation along the ΓK direction of graphene. Dashed lines show the energies for the constant energy cuts in b-e: (b) -1.1 eV, (c) -1.7 eV, (d) -2.3 eV, and (e) -3.45 eV. The hexagons in the constant energy cuts in panels b and c show (6 × 6) SiC unit cells. It can be appreciated that the spectral weight follows a periodicity of (6 × 6) SiC . Panels d and e show (√3 × √3)R30°S iC Brillouin zones corresponding to the symmetry of the spectral weight in this energy range. The Supporting Information includes some line profiles along these figures to make the observed periodicities more explicit.

Figure 4 .

 4 Figure 4. Top and side views of the buffer layer structure. (a) STM image (I = 0.25 nA, V = -1.5 V) of the buffer layer, showing the quasi-(6 × 6) SiC periodicity (black diamond). (b) Close-up on the image shows the graphene honeycomb structure without hexagon-pentagon-heptagon defects. (c) STEM-HAADF image of a cut through the buffer layer graphene-SiC interface. The dotted line is a guide to the eye above the buffer layer to show the bonding points of the buffer layer to the substrate. The overall periodicity between bonding points is ∼1.6 nm, compatible with a (6 × 6) SiC period. The atomic model in the inset below shows the DFTB relaxed model. (d) Vertical profiles of the interface in c have been obtained by averaging the STEM-HAADF data over several regions where the buffer is attached to the substrate ("contact") and where it is decoupled ("bump"). The red arrows in panel c indicate where the "contact" profiles were measured. (e) STEM-ABF image resolving Si and C atoms, showing the buffer layer on the right (z = 0). (f) Profile obtained from e, showing the contrast distribution for the last six SiC bilayers and the buffer layer. (g) Theoretical number of atoms versus the z coordinate (normal to the surface). The last SiC bilayer and the buffer layer have been indicated.

Figure 5 .

 5 Figure 5. Band structure of the (6√3 × 6√3)R30°S iC buffer layer reconstruction unfolded over (a) ΓKM and (b) along the perpendicular direction through K in the extended primitive graphene Brillouin zone (see the reciprocal space scheme in Figure 2a). Color scale represents the spectral weight. The graphene band structure has been superimposed for reference by aligning the deepest σ band. Vertical dotted lines in a represent the center of several Brillouin zone of the (6 × 6) SiC superperiodicity. An energy broadening of 10 meV has been applied to the bands. Left energy scales refer to the Fermi level as obtained in the supercell model, the right scales to the Fermi level in monolayer graphene. A close-up of the electronic structure near the Fermi level is shown in c for the ΓKM direction and in d for the perpendicular direction through K. The density of states for the buffer layer and for graphene is shown in panels e and f, corresponding to the same energy windows used in panels b and d, respectively.

ACKNOWLEDGMENTS This work was supported by the French Agence Nationale de la Recherche (ANR), project CoRiGraph, ref ANR-12-BS04-0017 and by the French American Cultural Exchange council through a Partner University Fund project. The authors acknowledge support from CASSIOPEE beamline staff and F. Nicolas from SOLEIL synchrotron. C.B. acknowledges partial financial support from the EU graphene Flagship program. Financial support was provided to W.d.H. and C.B. by the AFSOR under Grant No. FA9550-13-0217 and NSF No. 1506006. E.H.C. would like to acknowledge support from the NSF under Grant No. DMR-1401193.


Francisco Bernal
email: francisco.bernal@polytechnique.edu

Emmanuel Gobet
email: emmanuel.gobet@polytechnique.edu

Jacques Printems

Volatility uncertainty

Keywords: chaos expansion, uncertainty quantification, stochastic control, stochastic programming, Swing options, Monte Carlo simulations. MSC 2010: 93Exx, 62L20, 41A10, 90C15, 49L20

Introduction

Statement of the problem

Swing or Take-Or-Pay contracts are very common in energy markets. When the underlying commodity is gas, they give the possibility to the buyer of such a contract to buy/exchange the gas on the spot market (or the forward market) at a given price, with a flexibility on the cumulative volume on a given period, say [0, T ]. The valuation problem can be split into two steps: identifying the model from data, and computing the swing option price given the model. The second step is usually quite demanding since it requires to solve a non-linear Hamilton-Jacobi-Bellman (HJB) equation: for some related numerical methods, see [CSK01, JRT04, Kep04, BEBD + 06, BBP09]. Here we address the problem of swing valuation when the volatility of the model is unperfectly identified.

As a model, we consider the Schwartz-Smith model [START_REF] Schwartz | Short-term variations and long-term dynamics in commodity prices[END_REF] simplified to a single factor (see also [START_REF] Brennan | The price of convenience and the valuation of commodity contingent claims[END_REF]- [START_REF] Schwartz | The stochastic behavior of commodity prices: implications for valuation and hedging[END_REF]): namely, we assume that the forward gas price F (s, t) (fixed at time s for receiving and paying one unit of gas at time t > s) obeys the one-factor dynamics (for any maturity t) dF (s, t) F (s, t)

= σe -α(t-s) dW s , 0 ≤ s < t, (1.1)

given a initial forward curve (F (0, t) : t ≥ 0). The parameter σ > 0 is the volatility and α > 0 is the mean-reversion coefficient that describes the rate at which the short-term deviations are expected to disappear (see (2.2) later, for the link with mean-reverting Ornstein-Uhlenbeck process). The statistical inference of (α, σ) is usually handled with Maximum Likelihood estimation or Kalman filtering [Sch97, SS00, MT02]; standard statistical results [START_REF] Harvey | Forecasting, Structural Time Series Analysis, and the Kalman Filter[END_REF][START_REF] Rao | Statistical inference for diffusion type processes[END_REF] show that with a large number of data, the estimators are asymptotically normal (for the Local Asymptotic Normality property, see [START_REF] Gobet | LAN property for ergodic diffusion with discrete observations[END_REF]). Therefore, naturally, one may consider that the estimation errors on (α, σ) are approximately Gaussian, with small covariance/variances. Here to simplify the analysis, we consider that α is known and only σ is subjected to estimation error. In other words, we have

σ d ≈ N (µ, ν 2 ), (1.2)
with small variance ν > 0 and some µ > 0. This comes from a Central Limit Theorem result where typically ν 2 is proportional to the inverse of the number of available data.

The model-certain price v σ (0) at time 0 of the swing contract depends naturally on σ, and we wish now to account for the variability of price with respect to σ; in other words, to quantify the uncertainty in σ. Typically our aim is to efficiently estimate the distribution of the random variable v σ (0). Because v σ (0) solves a non-linear equation (for each given σ), it raises new issues that have not yet been addressed before.

From the point of view of financial economics, evaluating the mean excess between the uncertain swing option price v σ (σ) and the deterministic price v σ (0) is helpful to account for the risk of model uncertainty. This is similar to Credit Value Adjustment in finance, and other related adjustments, collectively called xVA (see [START_REF]Review of the credit valuation adjustment risk framework[END_REF]), corresponding (in the case of CVA) to accounting for the possibility of a counter-party's default within a bilateral contract. In our case, E [(v σ (σ) -z) + ] with z ≥ v σ (0) is referred to as UVA (Uncertainty Value Adjustment). The level z could be taken as v σ (0) = E [v σ (σ)], or as a given quantile of v σ (σ). This measure may serve to compute a deposit in order to face a quite unfavorable pricing of the swing option. Some illustrations are presented in Section 6.

Background results in uncertainty quantification

Many numerical strategies for uncertainty quantification are possible, of which we single out the most usual ones. For definiteness, we assume that a set of κ uncertain parameters σ is being propagated by a PDE, which are usually modeled as independent random distributions, and we say the random space is κ-dimensional. In the case of an uncertain function, modeled by a random field, dimension reduction is typically performed first via Karhunen-Loève expansions (KLE) truncated to a finite number κ of random parameters. Here are the different numerical strategies.

• Simulation-based methods, consisting in sampling M different σ according to (1.2) and for each of them, solving (v σ (m) (0) : 1 ≤ m ≤ M ). Sample statistics can then be performed on the ensemble of simulations. These methods are straightforward to implement provided that a deterministic solver for the PDE with fixed parameters is available (hence they are called non-intrusive). Moreover, they can handle arbitrary uncertainty distributions and high-dimensional uncertainty spaces with no particular modification. On the other hand, the convergence rate is very slow (typically as M -1/2 ), a problem which is compounded if each deterministic sample is costly to compute. Strategies to reduce the statistical error (and thus to speed up simulations) include Latin hypercube [START_REF] Loh | On latin hypercube sampling[END_REF], quasi Monte Carlo [START_REF] Niederreiter | Random number generation and quasi-Monte-Carlo methods[END_REF] and the Multilevel method [START_REF] Teckentrup | Further analysis of multilevel Monte Carlo methods for elliptic PDEs with random coefficients[END_REF].

• Local expansion-based methods, consisting in deriving Taylor/perturbation methods, taking advantage (or assuming) that the uncertain parameter σ fluctuates little. See [START_REF] Kleiber | The Stochastic Finite Element Method[END_REF]. It boils down to a sensitivity analysis which appears to be quite delicate in the case of HJB equation, because of the non-linear property of the equation (see (2.4) later).

• Generalized Polynomial Chaos (gPC) methods [START_REF] Xiu | Fast numerical methods for stochastic computations: a review[END_REF][LK10], which involve a functional expansion (the gPC expansion) of the solution in a basis in the random space-ideally orthogonal with respect to the joint distribution of parameters. There are two main approaches to gPC, differing in how the residual to the PDE is treated. In the first one-stochastic collocation-the residual of the gPC expansion is forced to vanish on a predefined set of nodes in the random space. This involves the discretization of the random space, which can be costly in high dimensions (κ 4)-although sparse grids can be sometimes used-but has the advantage of being non-intrusive. The second approach is called stochastic Galerkin, where the residual is forced to be orthogonal to each of the basis functions in the random space retained in the gPC expansion. This method is usually considered more accurate than stochastic collocation and does not involve grids in the random space; but the projection of the residual results in a system of κ PDEs which are different from the original problem (hence it is an intrusive method).

In this work we follow the last approach using polynomial chaos expansion (PCE) based on Gaussian noise (due to (1.2)). To the best of our knowledge, the above techniques have been investigated in the case of linear problems (like [START_REF] Mikulevicius | Linear parabolic stochastic PDE and Wiener chaos[END_REF]) and scarcely in the case of nonlinear PDEs (in fluid mechanics see [START_REF] Mikulevicius | Global L2-solutions of stochastic Navier-Stokes equations[END_REF], in stochastic optimization and HJB equations see [START_REF] Briand | Simulation of BSDEs by Wiener Chaos Expansion[END_REF][START_REF] Huschto | Solving stochastic optimal control problems by a Wiener chaos approach[END_REF], both cases are not related to Uncertainty Quantification strictly speaking).

Model and polynomial chaos expansion

Modeling the energy contract: HJB equation with uncertainty

To better represent (1.1), define

X s = σe -αs s 0 e αu dW u , (2.1)
which is an Ornstein-Uhlembeck process, starting from 0, and with zero long-term average:

dX s = -αX s ds + σdW s , X 0 = 0.
The process X is Gaussian and its limiting distribution as time goes to infinity is N (0, σ 2 2α ). Then, by the Itô formula we have

F (s, t) = F (0, t) exp e -α(t-s) X s - σ 2 4α e -2α(t-s) -e -2αt .
(2.2)

The gas spot price is

S t = F (t, t) = F (0, t) exp X t - σ 2 4α 1 -e -2αt .
(2.

3)

The swing contract has an expiration equal T > 0, and the fixed price to buy gas is given by the strike K > 0. At any date, the swing owner has the right to buy q t units of gas at the price K, so that his instantaneous profit is

f (t, X t , q t ) = q t (S t -K).
Here the volume flexibility is allowed in the range q t ∈ [q m , q M ] := Q ⊂ [0, +∞) at any time t ∈ [0, T ]. The cumulative gas withdrawn is

Y t = t 0 q s ds.
The value function is defined as

v σ (t, x, y) = sup (qs) t≤s≤T ,qs∈Q E P (X T , Y T ) + T t f (s, X s , q s )ds | X t = x, Y t = y (2.4)
where the expectation is taken under the risk-neutral measure defining the forward contract/spot price. The term P plays the role of a penalty forcing the cumulative withdraw Y T to finish in a given range (depending on the swing contract), typically we can take P with polynomial growth in its variables. Under smoothness assumptions on v, we know that v solves the following Hamilton-Jacobi-Bellman PDE:

   ∂ t v σ + sup q∈Q L σ,q v σ + f = 0, (t, x, y) ∈ [0, T ] × R × R, v σ (T, x, y) = P (x, y),
where

L σ,q Ψ(t, x, y) = σ 2 2 ∂ 2
x Ψ(t, x, y) -αx∂ x Ψ(t, x, y) + q∂ y Ψ(t, x, y).

By changing t T -t, we get another PDE (v σ (t, .) = u σ (T -t, .)) with initial condition:

   ∂ t u σ = σ 2 2 ∂ 2 xx u σ -αx∂ x u σ + sup q∈Q q ∂ y u σ + s(t, σ)e x -K , u σ (0, x, y) = P (x, y), (2.5)
where

s(t, σ) = F (0, T -t)e -σ 2 4α (1-e -2α(T -t) ) (see Equation (2.3))
. This is a degenerate diffusion along y compounded by advection in that direction (in the sup term).

Boundary conditions. Equation (2.5) is defined on the halfplane y ≥ 0. However, the solution is only needed in a bounded region I around the origin where contract scenarios with significant probability occur. For computational purposes, a finite computational domain D = [x m , x M ] × [0, y M ] with x m < 0 < x M must be defined which contains I. The boundary conditions (BCs) on ∂D are to be chosen in such a way that they do not affect the solution inside I; and they are compatible with the initial condition P .

Along the y direction, (2.5) is purely hyperbolic and, since q ≥ q m ≥ 0, information travels downward. Consequently, no BC is required at the downwind side of the boundary (y = 0), while the BC on the upwind side cannot affect I as long as it is farther away than q M T , i.e. if y M > sup y∈I +q M T .

Along the x direction, (2.5) is diffusive and information travels instantaneously; the goal is to place the BC far enough so that its effect inside I is negligible. In order to get a rough estimate of x L := max (-x m , x M ), let us heuristically split (2.5) into three parts: a diffusion term σ 2 2 ∂ 2 xx u σ ; an advection term -αx∂ x u σ , which advects information away from I; and a nonlinear term, which we neglect. Based on the fundamental solution of the diffusion part only, the effect b > 0 of a Dirichlet BC u(t , x L , y) on a point (x I , y I ) ∈ I at time t > t ("boundary pollution") is assumed bounded as

b |u(t , x L , y)| 2πσ 2 (t -t ) exp - |x I -x L | 2 2σ 2 (t -t ) . (2.6)
In this way, given a prescribed tolerance b, a rough estimate of x L can be obtained. Let us define Z σ (t, x, y) := ∂ y u σ (t, x, y) + s(t, σ)e x -K, (2.7) which may be denoted simply by Z σ when there is no ambiguity. Then

sup q∈Q qZ = q M (Z) + -q m (Z) -= Z q M 1 Z≥0 + q m 1 Z<0
and the HJB equation is

∂ t u σ = σ 2 2 ∂ 2 xx u σ -αx∂ x u σ + Z σ q M 1 Z σ ≥0 + q m 1 Z σ <0 , u σ (0, x, y) = g(y).
(2.8)

Hermite Polynomial chaos expansion

From (1.2), we assume an equality and we may write

σ 2 (ξ) = (µ + νξ) 2 , (Q)
leading to the approximation (1.2) as ν → 0. Here ξ is a standard Gaussian random variable N (0, 1), with density

p(ξ) = 1 √ 2π e -ξ 2 /2 .
The uncertain parameter being related to Gaussian distribution, we derive a PCE using Hermite polynomials. This is the purpose of the subsequent presentation.

Hermite polynomials

Consider the orthogonal basis of Hermite polynomials [LK10, Appendix B.1.2 p.502]: they are given by the Rodrigues formula

φ n (ξ) = (-1) n e ξ 2 2 ∂ n ξ [e -ξ 2 /2 ], n ≥ 0.
Therefore,

       φ 0 (ξ) = 1, φ 1 (ξ) = ξ, φ 2 (ξ) = ξ 2 -1, φ 3 (ξ) = ξ 3 -3ξ, φ 4 (ξ) = ξ 4 -6ξ 2 + 3, φ 5 (ξ) = ξ 5 -10ξ 3 + 15ξ, φ 6 (ξ) = ξ 6 -15ξ 4 + 45ξ 2 -15,
and

E [φ n (ξ)φ m (ξ)] = +∞ -∞ φ n (ξ)φ m (ξ)p(ξ)dξ = n!δ nm .
(2.9)

In particular, we set

V k := E φ 2 k (ξ) = k!.
(2.10)

These polynomials satisfy a three-term relation, as any family of orthogonal polynomials; here it takes the form:

ξφ n (ξ) = φ n+1 (ξ) + nφ n-1 (ξ), n ≥ 1, φ 1 (ξ) = ξφ 0 (ξ).
(2.11)

From this, we deduce

ξ 2 φ n (ξ) = φ n+2 (ξ) + (2n + 1)φ n (ξ) + n(n -1)φ n-2 (ξ), n ≥ 2, φ 3 (ξ) = (ξ 2 -3)φ 1 (ξ),
(2.12) that we will use subsequently to derive the system of non-linear PDEs to solve.

Chaos expansion, Monte-Carlo evaluation of the distribution of u

σ (•) Since u σ(ξ) : ξ ∈ R → L 2 (L 2 ([0, T ] × D L ; R), p(ξ)dξ)
, it can be decomposed onto the Hermite basis with coefficients in L 2 ([0, T ] × D L ; R) (see [START_REF] Le Maître | Spectral methods for uncertainty quantification: with applications to computational fluid dynamics[END_REF]):

u σ(ξ) (.) = +∞ k=0 ũk (.)φ k (ξ), (2.13)
with (for any t, x, y)

ũk (t, x, y) = E u µ+νξ (t, x, y)φ k (ξ) /k!, E |u σ(ξ) (t, x, y)| 2 = +∞ k=0 ũ2 k (t, x, y)k!.
(2.14) Thus, in the decomposition (2.13), (ũ 2 k (t, x, y)k!) k≥0 must form a convergent series for almost every (t, x, y) ∈ [0, T ] × D L . Following this PCE approach, our goal is now to identify the sequence of function coefficients ũk (•), k ≥ 0.

Once the coefficients identified (see Subsection 2.2.3 below), note that it is very easy (and efficient) to sample the distribution of a functional of u σ (•). For instance, if one is interested by the distribution of a single value u σ (t, x, y), it is enough to sample i.i.d. standard Gaussian ξ 1 , . . . , ξ M and to evaluate the approximation n k=0 ũk (t, x, y)φ k (ξ) as those sample points {ξ = ξ m , m = 1, . . . , M } and for some large order n: it will provide an empirical measure of u σ (t, x, y) under the uncertainty measure for σ. If instead of a single value one needs a functional of the solution as a function of (t, x, y) ∈ [0, T ] × D L , the methodology is unchanged. This approach is significantly faster than the naive one consisting to solve the whole PDE for each sampled ξ = ξ m , since only the φ k (ξ) have to be computed and the {ũ k : 0 ≤ k ≤ n} do not depend anymore on the noise.

If only the expectation of u σ (t, x, y) is required, it is enough to consider ũ0 (t, x, y) in view of (2.14), so that no extra simulation is needed. The variance is also obtained simply, without extra effort, by using the relation

Var(u σ(ξ) (t, x, y)) = +∞ k=1 ũ2 k (t, x, y)k! ≈ n k=1 ũ2 k (t, x, y)k!.
Higher order moments can be computable explicitly too, see [LK10, Appendix C]. This makes this PCE approach potentially much cheaper than the crude Monte-Carlo method.

Representation of the chaos expansion from the HJB equation

We truncate (Galerkin projection) the expansion (2.13) at the truncation level n ≥ 0:

u σ(ξ) (t, x, y) ≈ u ξ n (t, x, y) := n k=0 ũk (t, x, y)φ k (ξ).
The above residual at order n is orthogonal to {φ i (ξ)} n i=0 . From (2.8) and (2.14), we obtain a PC system (k = 0, 1, . . . , n)

E φ k (ξ)∂ t u ξ n (t, x, y) = E φ k (ξ) σ 2 (ξ) 2 ∂ 2 xx u ξ n (t, x, y) -E φ k (ξ)αx∂ x u ξ n (t, x, y) (2.15) + E φ k (ξ)Z ξ n (t, x, y) q M 1 Z ξ n (t,x,y)≥0 + q m 1 Z ξ n (t,x,y)<0 , (2.16)
where we set (similarly to (2.7))

Z ξ n (t, x, y) = Z t, x, y, u ξ n (t, x, y), ξ = ∂ y u ξ n (t, x, y) + s(t, σ(ξ))e x -K.
We may add a similar treatment of the initial and boundary conditions. The following transformations of (2.15)-(2.16) are obtained using extensively the orthogonality relation (2.9) and the definition of u n , assuming all the necessary smoothness on u (and thus on ũk ).

Term with time-derivative. The time-derivative and the first-order derivative in y (available in y) can be transformed into semi-Lagrangian derivative. In order to get rid of advection along a degenerate direction of the diffusion (y), let us define the semi-Lagrangian derivative along a characteristic Γ(t) : y(t) = y(0) + qt, t ≥ 0, as

D (q) t u(t, x, y) := D (q) t u(t, x, y ∈ Γ(t)) = ∂ t u(t, x, y) + q∂ y u(t, x, y).
It can be discretized for instance as

D (q) t u(t, x, y) ≈ u(t, x, y) -u(t -∆t, x, y -q∆t) ∆t . Set W + k (t, x, y) = E φ k (ξ)[s(t, σ(ξ))e x -K]1 Z ξ n (t,x,y)≥0 , W - k (t, x, y) = E φ k (ξ)[s(t, σ(ξ))e x -K]1 Z ξ n (t,x,y)<0 , H + jk (t, x, y) = E φ j (ξ)φ k (ξ)1 Z ξ n (t,x,y)≥0 . Since Z ξ n (t, x, y) = ∂ y u ξ n (t, x, y) + s(t, σ(ξ))e x -K, we derive E φ k (ξ)∂ t u ξ n -E φ k (ξ)Z ξ n (t, x, y) q M 1 Z ξ n (t,x,y)≥0 + q m 1 Z ξ n (t,x,y)<0 = E φ k (ξ)(∂ t u ξ n -q M ∂ y u ξ n (t, x, y))1 Z ξ n (t,x,y)≥0 + E φ k (ξ)(∂ t u ξ n -q m ∂ y u ξ n (t, x, y))1 Z ξ n (t,x,y)<0 -q M E φ k (ξ)[s(t, σ(ξ))e x -K]1 Z ξ n (t,x,y)≥0 -q m E φ k (ξ)[s(t, σ(ξ))e x -K]1 Z ξ n (t,x,y)<0 = E φ k (ξ)D (-q M ) t u ξ n (t, x, y)1 Z ξ n (t,x,y)≥0 + E φ k (ξ)D (-qm) t u ξ n (t, x, y)1 Z ξ n (t,x,y)<0 -q M W + k (t, x, y) -q m W - k (t, x, y) = E φ k (ξ)D (-qm) t u ξ n (t, x, y) + E φ k (ξ)(D (-q M ) t u ξ n (t, x, y) -D (-qm) t u ξ n (t, x, y))1 Z ξ n (t,x,y)≥0 -q M W + k (t, x, y) -q m W - k (t, x, y) = V k D (-qm) t ũk (t, x, y) + n j=0 H + jk (t, x, y)(D (-q M ) t ũj (t, x, y) -D (-qm) t ũj (t, x, y)) -q M W + k (t, x, y) -q m W - k (t, x, y).
(2.17)

We could also write the equation in terms of D (-q M ) t ũk (t, x, y). For this, set

H - jk (t, x, y) = E φ j (ξ)φ k (ξ)1 Z ξ n (t,x,y)<0 .
Then, similar computations give

E φ k (ξ)∂ t u ξ n -E φ k (ξ)Z ξ n (t, x, y) q M 1 Z ξ n (t,x,y)≥0 + q m 1 Z ξ n (t,x,y)<0 = V k D (-q M ) t ũk (t, x, y) + n j=0 H - jk (t, x, y)(D (-qm) t ũj (t, x, y) -D (-q M ) t ũj (t, x, y)) -q M W + k (t, x, y) -q m W - k (t, x, y). (2.18)
Interpolating the formulas (2.17) and (2.18) with a parameter θ ∈ R yields

E φ k (ξ)∂ t u ξ n -E φ k (ξ)Z ξ n (t, x, y) q M 1 Z ξ n (t,x,y)≥0 + q m 1 Z ξ n (t,x,y)<0 = V k D (-θqm-(1-θ)q M ) t ũk (t, x, y) + n j=0 [θH + jk (t, x, y) -(1 -θ)H - jk (t, x, y)](D (-q M ) t ũj (t, x, y) -D (-qm) t ũj (t, x, y)) -q M W + k (t, x, y) -q m W - k (t, x, y). (2.19)
Term with first space derivative. We simply have

E φ k (ξ)αx∂ x u ξ n (t, x, y) = αxV k ∂ x ũk (t, x, y).
(2.20)

Term with first second derivative. The transformation of E φ k (ξ) σ 2 (ξ) 2 ∂ 2 xx u ξ n (t,
x, y) depends on the model of uncertainty on ξ → σ(ξ). In the case (Q) where σ 2 (ξ) = (µ + νξ) 2 , in view of the three and four-terms relation (2.11)-(2.12) and using

V k = E [φ 2 k (ξ)] = k! (see (2.10)), we have (for 0 ≤ k ≤ n) E φ k (ξ) σ 2 (ξ) 2 ∂ 2 xx u ξ n = µ 2 2 V k ∂ 2 xx ũk + ν 2 2 n j=0 E ξ 2 φ j (ξ)φ k (ξ) ∂ 2 xx ũj + µν n j=0 E [ξφ j (ξ)φ k (ξ)]∂ 2 xx ũj , = µ 2 2 V k ∂ 2 xx ũk + ν 2 2 (k + 2)(k + 1)V k ∂ 2 xx ũk+2 + (2k + 1)V k ∂ 2 xx ũk + V k ∂ 2 xx ũk-2 + µν (k + 1)V k ∂ 2 xx ũk+1 + V k ∂ 2 xx ũk-1 = V k ν 2 2 ∂ 2 xx ũk-2 + µν∂ 2 xx ũk-1 + 1 2 µ 2 + (2k + 1)ν 2 ∂ 2 xx ũk + µν(k + 1)∂ 2 xx ũk+1 + (k + 2)(k + 1) 2 ν 2 ∂ 2 xx ũk+2 .
Summary. Now, plug the above into (2.16)-(2.19)-(2.20): after dividing by V k = k! > 0, the PC system in the Hermite basis with the semi-Lagrangian scheme approximated by forward finite differences is

(k = 0, 1, . . . , n) - ν 2 2 ∂ 2 xx ũk-2 -µν∂ 2 xx ũk-1 - µ 2 + (2k + 1)ν 2 2 ∂ 2 xx ũk -(k + 1)µν∂ 2 xx ũk+1 - (k + 2)(k + 1)ν 2 2 ∂ 2 xx ũk+2 + αx∂ x ũk + 1 ∆t ũk (t, x, y) - 1 ∆t ũk (t -∆t, x, y + [θq m + (1 -θ)q M ]∆t) + 1 k!∆t n j=0 [θH + jk (t, x, y) -(1 -θ)H - jk (t, x, y)] × ũj (t -∆t, x, y + q m ∆t) -ũj (t -∆t, x, y + q M ∆t) - q m k! W - k (t, x, y) - q M k! W + k (t, x, y) = 0.

Implementation

Experimentally, it turns out that the interpolation parameter θ in (2.19) does not noticeably affect the results. Thus, we fix θ = 1 henceforth. The PC system reads:

- ν 2 2 ∂ 2 xx ũk-2 -µν∂ 2 xx ũk-1 - µ 2 + (2k + 1)ν 2 2 ∂ 2 xx ũk -(k + 1)µν∂ 2 xx ũk+1 - (k + 2)(k + 1)ν 2 2 ∂ 2 xx ũk+2 + αx∂ x ũk + 1 ∆t ũk (t, x, y) - 1 ∆t ũk (t -∆t, x, y + q m ∆t) + 1 k!∆t n j=0 H + jk (t, x, y) ũj (t -∆t, x, y + q m ∆t) -ũj (t -∆t, x, y + q M ∆t) - q m k! W - k (t, x, y) - q M k! W + k (t, x, y) = B k (t, x, y), (3.1)
where for convenience we have introduced B k (t, x, y) = 0 (it will be nonzero with the benchmark problem to be introduced in Section 4). The above system poses the challenge that the nonlinear terms {H + ij , W ± k } couple all of the components of the solution, ũ0 , . . . , ũn . This means that the discretized version of (3.1) with finite differences (FD) will in principle have full bandwidth. Indeed, a well known drawback of PC is that the method may give rise to taxing algebraic equations [START_REF] Xiu | Fast numerical methods for stochastic computations: a review[END_REF], which is only compounded by the nonlinear nature of (3.1).

In order to (partially) circumvent this issue, we put forward an FD semi-implicit time-stepping scheme, whereby the linear terms are evaluated implicitly, and the nonlinear ones, explicitly (i.e. based on the solution on the previous timestep). Formally:

- ν 2 2 ∂ 2 xx ũk-2 m+1 F D -µν ∂ 2 xx ũk-1 m+1 F D - µ 2 + (2k + 1)ν 2 2 ∂ 2 xx ũk m+1 F D -(k + 1)µν ∂ 2 xx ũk+1 m+1 F D - (k + 2)(k + 1)ν 2 2 ∂ 2 xx ũk+2 m+1 F D + αx∂ x ũk m+1 F D + 1 ∆t ũk (t, x, y) m+1 F D - 1 ∆t ũk (t -∆t, x, y + q m ∆t) m+1 F D + 1 k!∆t n j=0 H + jk (t, x, y) m F D × ũj (t -∆t, x, y + q m ∆t) m F D -ũj (t -∆t, x, y + q M ∆t) m F D - q m k! W - k (t, x, y) m F D - q M k! W + k (t, x, y) m F D = B m+1 k . (3.2)
In (3.2), [Lv] m+1 F D and [Lv] m F D are proper FD discretizations of some operator L evaluated at current time t m+1 = (m + 1)∆t and at previous time t m = m∆t, respectively. In this way, the discretized system (3.2) has a pentadiagonal, rather than full, coupling. Let the spatial domain [x m = -x L , x M = x L ] × [0, y M ] be discretized into a grid of (1 + n x )(1 + n y ) nodes (x i , y j ) = (-x L + i∆x, j∆y), where ∆x = 2x L /n x and ∆y = y M /n y ; and let [U k ] m+1 ij be the FD approximation to ũk (t m+1 , x i , y j ). Furthermore, the random dimension is approximated as Ω ≈ [-L w , L w ] and discretized into n w + 1 equispaced nodes. We have employed the following FD approximations:

• central second-order-i.e. O(∆ 2 x)-FDs for the implicit second derivative,

∂ 2 xx ũk (t m+1 , x i , y j ) m+1 F D = [U k ] m+1 i+1,j + [U k ] m+1 i-1,j -2[U k ] m+1 i,j (
∆x) 2 ;

• second-order upwind FDs for the implicit first derivative along x, c∂ x ũk (t m+1 , x i , y j )

m+1 F D = max (c, 0) 3[U k ] m+1 ij -4[U k ] m+1 i-1,j + [U k ] m+1 i-2,j 2∆x 13 + min (c, 0) -[U k ] m+1 i+2,j + 4[U k ] m+1 i+1,j -3[U k ] m+1 ij 2∆x ;
• second-order central FDs for the explicit first derivative along y (notice that there is no need for an upwind scheme since it is computed explicitly),

Z(t m , x i , y j , ξ) m F D = n k=0 [U k ] m i,j+1 -[U k ] m i,j-1 φ k (ξ) 2∆y + s(t m , σ(ξ))e x i -K.
The nonlinear terms are approximated by numerical integration. Since they involve the indicator function, straightforward trapezoidal quadrature is preferable to higherorder methods, which are more prone to oscillations when dealing with discontinuous integrands (due to the indicator). For example,

H + sk ((t m+1 , x i , y j )) m F D = L w n w φ s (-w L )φ k (-w L )p(-w L )1 [Z(t m ,x i ,y j ,-w L )] m +φ s (w L )φ k (w L )p(w L )1 [Z(t m ,x i ,y j ,w L )] m + 2 nw r=1 φ s (ξ r )φ k (ξ r )p(ξ r )1 [Z(t m ,x i ,y j ,ξr)] m .
When choosing the timestep, we mind ∆t < q M ∆y-i.e. prevent the interpolation of semilagrangian derivatives from "jumping" over cells. (A CFL upper bound ∆t < 2∆x/|αx L | is not necessary since the derivatives in x are evaluated implicitly.)

The interpolation scheme is a critical aspect of semilagrangian differentiationspecially so if the former is to be monotone, see [START_REF] Warin | Some non monotone schemes for time dependent Hamilton-Jacobi-Bellman equations in stochastic control[END_REF]. In particular, the condition

1 ∆t ũk (t -∆t, x, y + q m ∆t) m+1 F D → 0 as ∆t → 0
must hold for consistency. (In our Matlab implementation, this was achieved with cubic-spline interpolation-but not necessarily so with the other available methods.) The resulting linear system to be solved at each time step has a block pentadiagonal structure (there are linear algebra routines tailored to such systems, see [START_REF] Benkert | An efficient implementation of the Thomas algorithm for block pentadiagonal systems on vector computers[END_REF]).

Hamilton-Jacobi PDEs such as (2.5) are well known to often lack a solution in the classical sense [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF]. For that reason, numerical schemes for solving scalar equations such as (2.5) must be designed to pick up the unique, sought-for viscosity solution. With FDs in particular, this is enforced with monotone discrete operators [START_REF] Barles | Convergence of approximation schemes for fully nonlinear second order equations[END_REF][START_REF] Tourin | An introduction to finite difference methods for PDEs in finance[END_REF] In contrast to the scalar case-and to the best of our knowledge-there is not a systematic theory of viscosity solutions for systems of equations such as (3.1) [START_REF] Liu | Viscosity solutions of fully nonlinear parabolic systems[END_REF]. Therefore, we will always assume that either a classical solution does exist, or that our numerical method is able to pick up the viscosity solution, otherwise. With this proviso, we have dropped the requirement of monotonicity of the finite difference operators. For instance, (3.3) and cubic-spline interpolation are not monotone, but they are more accurate than monotone versions, which typically are just O(∆x).

The code has been written in vectorized style in Matlab R2017 and run on a laptop.

Numerical validation

In order to assess empirically the accuracy and stability of the numerical scheme, we consider the PDE in (2.5) with an additional source term B:

∂ t u - σ 2 (ξ) 2 ∂ 2 xx u + αx∂ x u -Z q M 1 Z>0 + q m 1 Z<0 = B(t, x, y, ξ). (4.1)
Analogously to Section 2.2.3, we derive the PC equations for (4.1) up to truncation level n. For convenience, let us rewrite the resulting PC system (k = 0, . . . , n) without the semilagrangian approximation as

∂ t ũk - ν 2 2 ∂ 2 xx ũk-2 -µν∂ 2 xx ũk-1 - µ 2 + (2k + 1)ν 2 2 ∂ 2 xx ũk -(k + 1)µν∂ 2 xx ũk+1 - (k + 2)(k + 1)ν 2 2 ∂ 2 xx ũk+2 + αx∂ x ũk - 1 k! E φ k (ξ)Z t, x, y, u n (ξ) q M 1 Z>0 + q m 1 Z<0 = 1 k! E[φ k (ξ)B(t, x, y, ξ)]. (4.2)
The difference between (4.2) and (2.16) is the presence of the vector of source terms

E[φ k (ξ)B(t, x, y, ξ)]/k! (k = 0, . . . , n).
Let us choose as exact solution of the benchmark problem (or simply, the benchmark)

u ex (t, x, y, ξ) = nex k=0 ṽk (t, x, y)φ k (ξ), (4.3)
with

ṽk (t, x, y) = 1 + t (k + 1) √ k! sin 2πx 3 sin πy 2 ,
which complies with (2.14). Moreover, we let u ex (t, x, y, ξ) serve as initial condition and Dirichlet BCs for the benchmark problem (4.1). (Note that it can no longer be associated with the fair price of a gas swing option). We highlight from the onset that results turn out to be remarkably insensitive to discretization in the random dimension (provided it is adequate). Therefore, we set L w = 13, n w = 200 in all the ensuing experiments. Regarding the spatial domain, we set x m = -3, x M = 3 y M = 4, and take Dirichlet BCs along {|x| = 3, 0 ≤ y ≤ y M } and {x m ≤ x ≤ x M , y = y M }. The spatial discretization is discussed in terms of either n x and n y , or of ∆x = 6/n x and ∆y = 4/n y , equally. The remaining parameters are: F (0, t) = .15 1 + (sin 2πt)/3 , K = .1, q m = .1, q M = .7, µ = 1, α = .3 and (unless stated otherwise) ν = .3. The variables n and T remain to be set, as well as the value n ex for the u ex in (4.3). Given N ev = 500 fixed, uniformly scattered points across the domain, the error of the PC expansion is characterized by In a first set of experiments, we enforce the truncation error of the PC expansion to be zero by always setting n ex = n, and focus on the discretization error with respect to ∆t, ∆x and ∆y. This is shown both with n = n ex = 0 (left side of Figure 1), and with n = n ex = 16 (right side). A few remarks are in place:

RM S(t) = 1 N ev Nev j=1 n k=0 ṽk (t, x j , y j ) -ũk (t, x j , y j ) 2 .
• Using cubic-spline interpolation, consistency with respect to ∆t is clear (also for smaller values than shown on Figure 1). Furthermore, stability is apparent.

• While the accuracy noticeably improves with decreasing ∆t and ∆x, it is barely sensitive to ∆y.

• The pentadiagonal coupling notwithstanding, accuracy is scarcely affected by the number n + 1 of PC equations in the system. However, when n is very large, the evaluation of high-order Hermite polynomials may overflow (not shown). This is a well known issue which can be ameliorated by ad-hoc procedures, for instance by resorting to the Cauchy-Schwartz bound

|H ± jk (t, x, y)| ≤ E [|φ j (ξ)||φ k (ξ)|] ≤ j! √ k!.
Then, if the numerical evaluation of |H ± jk (t, x, y)| overshoots the above bound, H ± jk (t, x, y) must be replaced by sign H ± jk (t, x, y)

√ j! √ k!.
• With n = 0, the PC system is effectively an HJB equation. While the FD operators being used are not monotone, there is no need for them in the first place, since the benchmark solution is guaranteed to be smooth.

In Figure 1, the RMS error decreases roughly linearly w.r.t. ∆t as n x ∞ and n y ∞. In order to better assess the convergence rate of the numerical scheme, we produce an ensemble of simulations with discretization values randomly picked from uniform distributions in the ranges 0 ≤ n = n ex ≤ 20, 50 ≤ ∆x ≤ 200, 50 ≤ ∆y ≤ 200, and log (10 -3 ) ≤ log (∆t) ≤ log (.025). Then, we fit the resulting RM S(t = 1) errors to the nonlinear model

RM S(t = 1) = 1 + c 4 n c 1 ∆t + c 3 (∆x) c 2 ,
yielding as best-fit constants c 1 ≈ .26, c 2 ≈ 2.26, c 3 ≈ .55 and c 4 ≈ .0025, roughly in agreement with the expected second order convergence rate along x. This fit is plotted in Figure 2, and it turns out to be distinctly the best (in a least-squares sense) over several plausible nonlinear convergence models, including floating powers of ∆t and ∆y (not shown). (The outliers in Figure 2 are simulations with small ∆t and large ∆x.)

In the second part of this numerical study we focus on the effect of the truncation level n. Henceforth, we set n ex = 16 and allow n to differ from it.

In Figure 3, the longer-term behaviour of the RMS error is plotted. Observe that: • Asymptotically, all the instances of the RM S error tend to grow at a constant, sublinear rate (check the dashed lines of slope one on Figure 3). This means that relative errors are actually decreasing with increasing t (recall that the benchmark solution (4.3) grows linearly with t). Moreover, it strongly suggests stability of the numerical PC code for longer terminal times.

0 1 2 3 4 5 6 7 8 9 z = ( 1 + c 4 n )( c 1 t + c 3 x c 2 )
• When the PC solution lives in a smaller random space than the benchmark (i.e. when n < n ex = 16), errors are almost always larger when ν is larger. On the other hand, when n = n ex = 16, the RM S error is virtually insensitive to ν. Moreover, truncating beyond n ex does not improve accuracy (as expected).

A well known issue of PC is that the moments of the PC solution (specially from the variance onwards) tend to deteriorate rather rapidly as t ∞, unless the truncation level is also constantly raised or the random basis adaptively changed [START_REF] Gerristma | Timedependent generalized polynomial chaos[END_REF]. Figure 4 shows the propagation of the error of the mean and that of the variance (with respect to the mean and variance of the benchmark solution) at (x 0 , y 0 ) = ( √ 2, √ 3). The relative error of the variance grows faster than that of the expectation and-unless n is very close to n ex -becomes O(1) by time t 2.

Finally, in Figure 5 we study the dependence of the PC numerical solution on the random variable ξ ∼ N (0, 1). We focus on the interval -3 ≤ ξ ≤ 3 which has an aggregate probability of over 99%. When n n ex = 16, the error (with respect to the benchmark) is small and increases with T and with (∆t, ∆x, ∆y). On the other hand, when n is too low, the truncation error overrides that of discretization and can be unacceptably large (O(1) and bigger). Summing up, the numerical scheme (3.2) has been validated with a smooth solution.

Uncertainty quantification on gas swing option

We consider the gas swing option priced under one-factor dynamics in [START_REF] Bardou | Optimal quantization for the pricing of swing options[END_REF]. The parameters are: T = 1 (in years) , µ = .7, α = 4, F 0t = 20 (a constant), q m = 0, q M = 2190, Q max = 1900, Q min = 1300 and K = {5, 10, 15, 20}. In a slight departure from our own framework, the contract described there does not allow to adjust the gas drawing rate continuously, but only once a day over one year (i.e. it depends on 365 discrete variables rather than on a continuous function y(t)). There are two cases:

• Unconstrained contract: in this scenario, there are no upper or lower limits to the total amount of gas that can be withdrawn over the entire contract period (the swing option contract is in fact equivalent to a string of 365 call options). Therefore, P (x, y) = 0 in (2.5) (so that the initial condition is smooth).

• Constrained contract: this is the most usual contractual scenario, when there are strict limits on the total amount of gas that can be withdrawn, i.e.

Q min ≤ Y T ≤ Q max .
(5.1) Strictly speaking, (5.1) amounts to a penalty function which is zero for Q min ≤ y ≤ Q max and infinity outside that interval. For numerical purposes, it is approximated by the double ramp function where A is a user-defined positive constant. We take A = 10000 as recommended by [START_REF] Bardou | Optimal quantization for the pricing of swing options[END_REF]. Note that A = 0 recovers the unconstrained case.

P (x, y) = A s(0, σ)e x min (0, y -Q min ) -max (0, y -Q max ) , 20 -3 -2 -
The solution is to be evaluated at y T = 0 (no initial gas) and

x T = µ 2 4α (1 -e -2αT ) + log F 0t (T ) F 0t (0) ≈ .03061.
The initial conditions for the PC equations now read:

ũk (0, x, y) = Ae x k! min (0, y -Q min ) -max (0, y -Q max ) E s σ(ξ), 0 φ k (ξ) , (5.2)
which is again approximated by trapezoidal quadrature for simplicity. Numerical experiments show that taking L w = 10 and n w = 100 is sufficient in all the cases. Due to the nonsmooth initial condition, in this problem linear interpolation performs better than cubic splines, so that we use the former instead. x M 15µ/ √ 2α for the solution to be unaffected by the finite computational domain.

Computational domain and BCs. As discussed in Section 2.1, the computational domain D = [x m , x M ]×[0, y M > q M T ] must be large enough that the numerical solution at point (x T , y T ) at time T is not affected by the BCs enforced up to time T . We set y M = 1.3q M T = 2847, and choose as BCs: u µ+νξ n (t, x, y) = P (x, y), if (x, y) ∈ ∂D and t = 0, 0, if (x, y) ∈ ∂D and 0 < t ≤ T .

In this way, only the effect of the initial condition along ∂D must be reckoned with. Owing to the e x factor in (5.2), u µ+νξ n (T, x T , y T ) is rather insensitive to x m < 0 (we set x m = -1.24), but x M > 0 must be handled with more care. The reason is that the influence of the initial condition on the boundary decays as O(e -x 2 M /(2µ 2 ) ), but the BC itself grows as O(e x M ). The value of x M is estimated as follows: we wish to keep the "boundary pollution" (let us call it b) smaller than of the solution at (T, x T y T ). According to Table 1

(more on it later), b = 30 > .001 × E[u µ+νξ n (T, x T , y T )] ≥ .001 × |ũ k (T, x T y T )| (0 ≤ k ≤ n) is a conservative bound.
Combining the estimate (2.6) and the BC (at t = 0) given by (5.2) yields x M ≈ 15µ/ √ 2α ≈ 3.72 (i.e. 15 standard deviations of the OU process (2.1)). Figure 6 shows the effect of the "boundary pollution" for several x M . unconstrained (A = 0) constrained (A = 10000) K u µ 0 (T, x T , y Using the PC code in "scalar mode". The PC algorithm can, in principle, also be used for scalar HJB equations such as (2.5)-and hence for pricing swing options without uncertainty, too. This is possible by letting n = 0, ν = 0, and µ = σ in (3.1). In the sequel, u σ 0 (•) is to be understood as the numerical solution of an HJB equation with volatility σ obtained by so using the PC code (in "scalar mode").

Since some of the FD operators in Section 3 are not monotone, u σ 0 (•) is not guaranteed to pick up viscosity solution. (This is particularly relevant in the constrained case, where the initial conditions are non-differentiable.) On the other hand, comparison with the results of a separate viscosity-capturing HJB code (based on Howard's algorithm [START_REF] Tourin | An introduction to finite difference methods for PDEs in finance[END_REF]) shows no difference with u σ 0 (•)-while the latter enjoys a higher convergence rate w.r.t. ∆x. Therefore, we stick to the PC code henceforth.

In order to get a feeling about the several discretization parameters involved, we calculate the swing option prices for K = {5, 10, 15, 20} with no uncertainty (i.e. u σ 0 (T, x T , y T )), and compare it with the values given in [START_REF] Bardou | Optimal quantization for the pricing of swing options[END_REF]. In Table 1, u EX has been calculated with the Black-Scholes formula, and u RQ with a quantization method (see [START_REF] Bardou | Optimal quantization for the pricing of swing options[END_REF] for details). In the remaining of the paper, we will use either the

• "coarser discretization": ∆t = 1/365, n x = 160 (and n y = 200), or the • "finer discretization": ∆t = 1/(2 × 365), n x = 320 (and n y = 200 as well).

Reference solution. In order to assess the quality of the PC solution, we will fix ν first, and then compare u µ+νξ n (T, x T , y T ) inside a given interval of ξ with the HJB solution u σ 0 (T, x T , y T ), where σ = µ + νξ. Specifically, for each K we precompute a set of {u σ i 0 (T, x T , y T )} 90 i=1 HJB solutions with the "finer" discretization at 90 equispaced volatilities σ 1 = 0.1, . . . , σ 90 ≈ 1.75. For a given ν, the comparable interval of ξ is thus (σ 1 -µ)/ν ≤ ξ ≤ (σ 90 -µ)/ν. 2. (ν = .5; "finer" discretization)

We start the numerical study by investigating the interplay between the discretization and the truncation level, n. In order to better highlight the effect of n, we first pick the relatively large value ν = .5. The dashed curves in Figure 7 are based on a subset of 70 of those HJB simulations. They represent u σ 0 (T, x T , y T ) in the interval 0.1 ≤ σ ≤ 1.3 (i.e. centered around µ). It can be seen that the delta of the swing option is nearly constant.

The red curves in Figure 7 are the PC approximation u µ+νξ n (T, x T , y T ) for n = {2, 4, 8}. The fit to the reference HJB curves improves as n grows. In order to systematically check this, we track the RMS error between each pair of curves, defined as (we include a superindex K to specify the pair of curves):

K n = 2 n = 4 n = 8 n =
R K = 1 70 70 i=1 u K,σ i 0 (T, x T , y T ) -u K,σ i n (T, x T , y T ) 2 .
(5.4)

The values of R K for the "finer" and "coarser" discretization are listed on Table 2. As expected, R K usually improves (i.e. diminishes) as the discretization is refined and n grows. The exception are the entries with K = 5, for which neither discretization is fine enough at ν = .50. (Since Hermite polynomials oscillate with ever larger amplitude as |ξ| increases, numerical errors accumulate at the tails of the higher-order components of u µ+νξ n (•)-check also Figure 5.) Figure 8 zooms in on the fit between u µ+νξ n (T, x T , y T ) and HJB reference curve with K = 15, for n = {2, 4, 8, 16} and the "finer" discretization. It is apparent that as n grows, the errors accumulate on the higher end of ξ. Therefore, the truncation level can be increased only so much for a given discretization. For this reason, when computing the cdf of the swing prices under uncertainty, the integration w.r.t. ξ must be cut short above. Furthermore, since σ = µ + νξ must be positive, swing prices for ξ < -µ/ν are meaningless. In sum, we introduce the following truncated cdf:

CDFT[u µ+νξ n , z] := 1 √ 2π ξ M max (ξm,-µ/ν) e -ξ 2 /2 1 u µ+νξ n (T,x T ,y T )≤z dξ,
where ξ m and ξ M are respectively negative and positive enough that the distribution is well represented. (Note that this may not be possible if the cut-off |µ/ν| is too small). By differentiating CDFT[u µ+νξ n

, z] with respect to z, one obtains the approximate pdf of the Swing option price under uncertainty. The error tends to accumulate at the tails, thus having a larger adverse effect on the higher moments of the distribution. This is illustrated in Figure 9, where the more realistic uncertainty ν = .1 has been In Figure 9, and in all the remaining simulations in the paper, the following values have been used: K = 15, ξ m = -3.5, ξ M = 3.5. (We point out that no smoothing has been applied either to the truncated cdfs or to the pdfs.)

Uncertainty Value Adjustment

Consider now the following practical situation. The swing option is sold at the price v σ (0) corresponding to a medium estimation of the volatility: this may be an equilibrium price in the market. However, the seller of the contract wants to assess the risk for volatility misspecification: he could consider a worst-case analysis by computing sup σ v σ (σ) but this may be too conservative. More reasonably, he incorporates the distribution on σ (coming from his calibration/estimation procedure) and aims at computing E [(v σ (σ) -z) + ] with z ≥ v σ (0). The above quantity can be seen as an Uncertainty Value Adjustment (UVA), whose purpose is to give a financial evaluation of the possibility of misspecification of σ in the swing problem, beyond a risk threshold z. In finance (like for CVaR), z is a given by a quantile at 95% or 99%.

We proceed to explore this idea using the gas swing option in the previous Section. Let us consider three quantiles q = {99%, 95%, 90%}. For each of them, and given a PC expansion u µ+νξ n (T, x T , y T ), the risk threshold z q is defined by CDFT[u µ+νξ n , z q ] = q/100.

After solving for z q with Newton's method, the UVA is calculated as

UVA[u µ+νξ n , q] = 1 √ 2π ξ M
max (ξm,-µ/ν) min 0, u µ+νξ n -z q e -ξ/2 dξ. , q]/E[u µ+νξ n (T, x T , y T )] %. For the Monte Carlo (MC) references, the entries are 100×UVA[u σ 0 , q]/E[u σ 0 (T, x T , y T )] %.

swing option. They are compared with the "Monte Carlo" reference based on stochastic collocation of the HJB equation. A relatively coarse discretization and low n are enough to get adequate acccuracy. With µ = .7 and ν = .25 the cut-off value |µ/ν| = 2.8 is already small (less than 3 standard deviations of the distribution of ξ). In order to investigate UVA at larger uncertainties than 10%, we set µ = 1.3 (and keep K = 15, ξ m = -3.5, ξ M = 3.5, and the other parameters as before). See Table 4 for the UVA percentages in that case. ν q = 99% q = 95% q = 90% 10% . 0062 Table 4: UVA (in terms of a percentage over the expected price of swing option) for µ = 1.3, K = 15, and the other parameters as in Section 5. (The "coarse" discretization was employed.)

We conclude that settling for n = 2 and the "coarser" discretization is adequate for an accurate estimation of UVA. On a laptop and in Matlab, it takes about 3-4 minutes to calculate u µ+νξ

Figure 1 :

 1 Figure 1: Convergence of RM S(t = 1) w.r.t. ∆t, ∆x and ∆y (separately), with n = n ex .

Figure 2 :

 2 Figure 2: Best fit of RM S(t = 1) error to a nonlinear convergence rate model. Colours: n = 0 (white), 1 ≤ n ≤ 6 (green), 7 ≤ n ≤ 13 (blue), 14 ≤ n ≤ 20 (red); and n ex = n.

Figure 3 :

 3 Figure 3: Propagation of the RM S error. Auxiliary dashed lines have slope 1.

Figure 4 :

 4 Figure 4: Error of the expectation (left) and variance (right) of the PC expansion for various truncation levels (∆t = .005, n x = 200, n y = 100). Dashed lines have slope 1.

Figure 5 :

 5 Figure 5: Relative error in the random dimension with respect to the benchmark (n ex = 16) at (x 0 , y 0 ) = ( √ 2, √ 3) and: t = 1 (left), t = 5 (middle) , and t = 10 (right). Dashed curves represent a rougher discretization (∆t = .02, n x = 90, n y = 60). Solid curves are a finer discretization (∆t = .005, n x = 150, n y = 100). For n = 4, both sets are nearly indistinguishable to the naked eye.

Figure 6 :

 6 Figure 6: Illustration of the "boundary pollution" effect with K = 5 and n = 8:x M 15µ/ √ 2α for the solution to be unaffected by the finite computational domain.

Figure 7 :

 7 Figure7: Comparison between the PC expansion with increasing n, and the HJB reference curves, at (T, x T , y T ). With n > 4, the improvement becomes harder to tell with the naked eye; use the zoom and/or check Table2. (ν = .5; "finer" discretization)

Figure 8 :

 8 Figure 8: Comparison of the PC expansions ("finer" discretization) with increasing n to the reference HJB curve for K = 15 and ν = .5. For the sake of clarity, only a few individual values of the PC curves are displayed instead of the continuous curves.

Figure 9 :

 9 Figure 9: Pdf of the swing option price (with K = 15) under 10% uncertainty on the volatility. The error of discretization overrides that of truncation. The black curve is the reference "Monte Carlo" pdf, based on stochastic collocation of the HJB equation.

Table 1 :

 1 Using the PC algorithm (3.2) with n = 0 and ν = 0 as an HJB solver for pricing gas swing options without uncertainty. (∆t = 1/365, n x = 320, n y = 200, x M = 3.72). References u EX and u RQ are the solution of the slightly different problem in[START_REF] Bardou | Optimal quantization for the pricing of swing options[END_REF].

		T ) u EX	u µ 0 (T, x T , y T ) u RQ
	5 33274	32760	29361	29342
	10 22298	21844	19987	19866
	15 11760	11381	10867	10698
	20 4181	3966	2580	2680

Table 2 :

 2 Improvement of the PC solution as n grows and the discretization is refined, as measured by R K (5.4). For K = 5 and n = 16, both discretizations are inadequate.

	16

Table 3 :

 3 Table 3 compiles UVA approximations for K = 15 and ν = 10% using several PC expansions. The entries are the percentage of UVA relative to the mean price of the ν = .10 and K = 15 q MC n = 2 n = 4 n = 8 n = 16 99% .0067 .0072 .0070 .0071 .0069 .0071 .0069 .0071 .0069 95% .0453 .0476 .0463 .0470 .0457 .0470 .0457 .0470 .0457 90% .1031 .1075 .1045 .1065 .1036 .1066 .1035 .1066 .1036 disc. coarser finer coarser finer coarser finer coarser finer Relative (percentual) UVA: 100×UVA[u µ+νξ n

  .0062 .0424 .0424 .0973 .0974 20% .0126 .0128 .0857 .0858 .1964 .1960 30% .0197 .0204 .1322 .1323 .3014 .2990 40% .0267 .0276 .1832 .1797 .4156 .4042

	n	2	4	2	4	2	4

(•) and the UVA. While the risk premium against misspecification of the volatility is quite small when the calibration is precise, it grows roughly linearly as the calibration uncertainty increases.

* This research is part of the Chair Financial Risks of the Risk Foundation, the Finance for Energy Market Research Centre (FiME) and the ANR project CAESARS (ANR-15-CE05-0024).