Reducing turbulent transport in tokamaks by combining intrinsic rotation and the low momentum diffusivity regime

Abstract

Based on the analysis of a large number of high-fidelity nonlinear gyrokinetic simulations, we propose a novel strategy to improve confinement in tokamak plasmas by combining up-down asymmetric flux surface shaping with the Low Momentum Diffusivity (LMD) regime. We show that the intrinsic momentum flux driven by up-down asymmetry creates strong flow shear in the LMD regime that can significantly reduce energy transport, increasing the critical gradient by up to $25\%$. In contrast to traditional methods for generating flow shear, such as neutral beam injection, this approach requires no external momentum source and is expected to scale well to large fusion devices. The experimental applicability of this strategy in spherical tokamaks is addressed via simulations by considering actual equilibria from MAST and a preliminary equilibrium from SMART.

Introduction–It is well known that strong toroidal rotation (specifically $E\times B$ rotation shear) can reduce turbulent transport [1, 2, 3, 4, 4, 5, 6, 7, 6, 8, 9, 10, 11, 11, 12, 13, 14, 11, 15, 16, 17, 18, 19] and stabilize MagnetoHydroDynamic (MHD) instabilities [20, 21, 22, 23, 24, 21, 22], thereby improving the prospects for tokamak fusion power plants. Traditional methods of driving plasma rotation such as Neutral Beam Injection (NBI) [25, 26, 27] or Radio Frequency (RF) waves [28, 29, 30, 31, 32] are not expected to scale well to large devices [33]. Taking NBI as an example, we see that the ratio of the momentum $p_{\text{NBI}}$ over the energy $E_{\text{NBI}}$ injected by the beam scales as $p_{\text{NBI}}/E_{\text{NBI}}\sim 1/\sqrt{E_{\text{NBI}}}$ [34]. Thus, the more energetic beams needed to penetrate into larger and denser plasmas will be less effective at driving rotation.

An intriguing alternative is to utilize intrinsic rotation [35, 36, 37], spontaneously driven by the plasma turbulence itself under certain conditions. The symmetry properties of gyrokinetics [5, 38, 39, 40] show that the only way to create flow of the order of the sound speed in the tokamak core is to break the up-down symmetry of magnetic flux surfaces about the midplane [41, 42]. This approach is independent of external sources and is expected to scale well to large devices (as it persists in the limit $\rho^{*}=\rho_{i}/a\to 0$, where $\rho_{i}$ is ion gyroradius and $a$ is tokamak minor radius). One can calculate the radial profile of rotation in the steady-state operation of a tokamak by requiring the total toroidal angular momentum flux to be zero, which is well approximated by considering only the ion contribution ($\Pi_{i}=0$) due to their large mass compared to electrons. Since the momentum flux depends on both the toroidal ion rotation $\Omega_{i}$ and its radial derivative $d\Omega_{i}/dx$, it is common to Taylor expand $\Pi_{i}$ about $\Omega_{i}=0$ and $d\Omega_{i}/dx=0$ to find [41]

$$\Pi_{i}\left(\Omega_{i},\frac{d\Omega_{i}}{dx}\right)\simeq\Pi_{i,int}-n_{i}m_{i}R^{2}_{0}D_{\Pi i}\frac{d\Omega_{i}}{dx}-n_{i}m_{i}R^{2}_{0}P_{\Pi i}\Omega_{i}=0,$$ (1)

where $\Pi_{i,int}=\Pi_{i}(0,0)$ is the intrinsic momentum flux generated by up-down asymmetry and the second and third terms are the diffusive and pinch terms, respectively. Here, $n_{i}$ is ion density, $m_{i}$ is the ion mass, $R_{0}$ is the major radius of the tokamak, $x$ is the radial coordinate, $D_{\Pi i}$ is the angular momentum diffusion coefficient and $P_{\Pi i}$ is the angular momentum pinch coefficient. Since the pinch term is typically smaller than the diffusive term [43, 44, 45, 46], we neglect it for now, but show in Ref. [47] that this will only cause an underestimation of the flow shear created by up-down asymmetry. The strength of the ion momentum diffusivity is typically compared with the ion heat diffusivity $D_{Qi}$ through the ion Prandtl number

$$\text{Pr}_{i}=\frac{D_{\Pi i}}{D_{Qi}}.$$ (2)

This allows one to estimate the rotation shear being driven by a given turbulence level. A low Prandtl number means that a given momentum source (external or intrinsic) drives larger rotation shear at a given turbulence level. Without considering the pinch term, Eq. (1) indicates that maximizing $\Pi_{i,int}$ or minimizing ion momentum diffusivity will strengthen velocity shear. Maximizing $\Pi_{i,int}$ has been the subject of previous work [42]. The focus of this letter is to minimize the Prandtl number so as to help stabilize turbulence. In conventional aspect ratio tokamaks, $\text{Pr}_{i}\simeq 1$ has been experimentally observed [48, 49], motivating many theoretical works to assume this [50, 51, 41, 42]. Recently, however, it has been shown by gyrokinetic simulations that a significantly lower Prandtl number can be achieved [52, 53], especially at tight aspect ratio and low safety factor [54, 55, 56]. However, the Prandtl number calculations in many works were simplified, since either parallel momentum flux [54] or the toroidal angular momentum flux generated by flow shear parallel to the magnetic field [55, 56] was calculated. A more complete calculation must consider all components of toroidal angular momentum flux [38, 13] generated by toroidal flow shear as actually present in tokamaks [57]. A thorough scan of the dependence of the Prandtl number on the various geometric properties of tight aspect ratio tokamaks has not been performed, nor has a self-consistent study of the stabilization effect of intrinsic flow shear by up-down asymmetry on plasma turbulence. Despite recent experimental tokamak data showing that one can obtain a Prandtl number lower than one [7, 58, 59, 44, 60, 45, 46], these results are often hard to interpret given the challenge of separating the three terms in Eq. (1), especially for tight aspect ratio spherical tokamaks [61, 45]. All these theoretical and experimental issues call for a comprehensive theoretical study of toroidal angular momentum transport and how this can be combined with rotation drive mechanisms to stabilize turbulence.

We perform a thorough study of the Low Momentum Diffusivity (LMD) regime using a large number of high-fidelity nonlinear (NL) local gyrokinetic simulations. We first use a circular geometry to determine the parameter regime that achieves LMD. We then consider a tilted elliptical geometry in order to drive rotation from up-down asymmetry. By scanning the flow shear $d\Omega_{i}/dx$ to find the value at which momentum flux is zero, we show that turbulence can self-consistently drive an intrinsic flow shear that significantly reduces the heat flux. The experimental applicability of this approach is addressed via simulations of a MAST tokamak equilibrium and a preliminary SMART equilibrium. To our knowledge, this work not only provides the first comprehensive numerical study on toroidal angular momentum diffusivity in the LMD regime, but also demonstrates the first potentially practical method of driving strong rotation shear in future tokamak power plants.

Methods–We use the well-benchmarked code GENE [62, 63] to perform a large number of flux tube (local) gyrokinetic simulations and model core plasma turbulence. We consider the Miller representation [64] to parameterize the flux surfaces. To make our analysis as accurate as possible, we consider the actual toroidal flow shear (including its parallel and $E\times B$ components) to drive momentum flux as well as calculate the exact toroidal angular momentum flux (according to Eq. (A.2) in Ref. [65]).

Analysis–We start with circular concentric flux surfaces, an inverse aspect ratio $\epsilon=r_{0}/R_{0}=0.36$ ($r_{0}$ is the minor radial location of the flux tube), adiabatic electrons, equal electron and ion temperatures, and a density gradient of $R_{0}/L_{n}=2.22$. These parameters are the same as the Cyclone Base Case (CBC) [66] (except the tighter aspect ratio) and produces Ion Temperature Gradient (ITG) driven turbulence. We perform a parameter scan over safety factor $1.05\leq q\leq 4.55$, magnetic shear $0.1\leq\hat{s}\leq 1.6$ and ion temperature gradient $4.96\leq R_{0}/L_{Ti}\leq 12.96$, comprising a 3D data set with 576 NL gyrokinetic simulations, where $L_{n}$ and $L_{T}$ are the characteristic gradient scale lengths of density and temperature, respectively. To drive non-zero toroidal angular momentum flux and compute the momentum diffusivity, we impose toroidal flow shear $d\Omega_{i}/dx$ including a parallel and a perpendicular component of $\omega_{\perp}=-(r_{0}/q)d\Omega_{i}/dx=0.12c_{s}/R_{0}$, where $c_{s}=\sqrt{T_{e}/m_{i}}$ is the sound speed, and $T_{e}$ is electron temperature. We use a simulation box with widths $L_{x}\times L_{y}=200\rho_{i}\times 200\rho_{i}$ and number of grid points $(n_{k_{x}},n_{k_{y}},n_{z},n_{v_{||}},n_{\mu})=(192,64,64,32,9)$, where $\rho_{i}$ is the ion gyroradius, $(x,y,z,v_{||},\mu)$ are the radial, binormal, parallel, parallel velocity, and magnetic moment directions, respectively, $k_{x}$ and $k_{y}$ are the wavevectors in $x$ and $y$ directions. The heat flux and toroidal angular momentum flux, $\hat{Q}_{i}$ and $\hat{\Pi}_{i}$, are respectively normalized to gyroBohm units $c_{s}n_{i}T_{i}(\rho_{i}/R_{0})^{2}$ and $c^{2}_{s}m_{i}n_{i}R_{0}(\rho_{i}/R_{0})^{2}$, where $T_{i}$ is the ion temperature.

Figure 1 (a) shows the constant Prandtl number manifold $\text{Pr}_{i}=0.5$ that we construct by linearly interpolating the $\text{Pr}_{i}$ data from our 3D parameter scan. In simulations with up-down symmetric geometries, the Prandtl number is calculated by

$$\text{Pr}_{i}=\frac{\hat{\Pi}_{i}}{\hat{Q}_{i}}\frac{R_{0}}{L_{Ti}}\frac{\epsilon}{q\,\omega_{\perp}}\frac{c_{s}}{R_{0}}.$$ (3)

Below the manifold, corresponding to lower values of $q$ at fixed $\hat{s}$ and $R_{0}/L_{Ti}$, $\text{Pr}_{i}$ is smaller than $0.5$. We choose to define this region as the LMD regime. Only the nonlinearly unstable region of the parameter space is presented in the figure. We see that $\hat{s}\simeq 1$ and low $R_{0}/L_{Ti}$ favor a low $\text{Pr}_{i}$, reflected by a larger range in $q$ with low $\text{Pr}_{i}$. The heat flux shown by the color of the manifold changes primarily with $R_{0}/L_{Ti}$, indicating that $q$ and $\hat{s}$ do not significantly affect $\hat{Q}_{i}$ in the LMD regime. To show the parameter dependence of the Prandtl number more clearly, several contours of $\text{Pr}_{i}$ are shown in Fig. 1 (b)-(d) by holding either $q$, $\hat{s}$ or $R_{0}/L_{Ti}$ constant. It is clear that $\text{Pr}_{i}$ has a non-trivial dependence on these parameters. Figure 1 (c) and (d) show that $\text{Pr}_{i}$ depends strongly on $q$, and Fig. 1 (b) shows that at low $q$, $\text{Pr}_{i}$ becomes small over a large range in $\hat{s}$ and $R_{0}/L_{Ti}$. This feature, together with a small aspect ratio $1/\epsilon=2.78$ suggests that a more poloidally angled magnetic field reduces $\text{Pr}_{i}$. While the dependence of $\text{Pr}_{i}$ on $R_{0}/L_{Ti}$ is weaker, being closer to marginal stability reduces $\text{Pr}_{i}$. The Prandtl number depends strongly on magnetic shear, with a factor of $\sim 2$ decrease giong from $\hat{s}=0.1$ to $\hat{s}\simeq 1$, as shown in Fig. 1 (d). Though not shown here, only considering the parallel component of $\hat{\Pi}_{i}$ results in an overestimate of $\text{Pr}_{i}$ [47]. Adding kinetic electrons or changing to Trapped Electron Mode (TEM) turbulence retains the above features in $\text{Pr}_{i}$, while increasing the aspect ratio $1/\epsilon$ strongly increases $\text{Pr}_{i}$ [47].

Figure 1: Visualizations of the Prandtl number. (a) shows the surface with $\text{Pr}_{i}=0.5$, below which we define as the LMD regime. The color map on the surface denotes the heat flux $\hat{Q}_{i}$. The contours in the bottom plane represent lines of constant $q$ values on the surface. (b)-(d) provide contours of Prandtl number fixing one of our three parameters: (b) $q=1.55$, (c) $\hat{s}=1.0$, and (d) $R_{0}/L_{Ti}=6.96$. The red contour lines denote $\text{Pr}_{i}=0.5$.

Figure 2: Subplots (a) and (b) show $\hat{\Pi}_{i}$ and $\hat{Q}_{i}$ as a function of $\omega_{\perp}$ obtained with NL simulations for $\hat{s}=0.1$ (red), $\hat{s}=0.4$ (green), $\hat{s}=0.8$ (blue) and $\hat{s}=1.2$ (black). Subplots (c) and (d) show a comparison of $\hat{\Pi}_{i}/\hat{Q}_{i}$ as a function of $\omega_{\perp}$ between the QL model (black), NL simulations with four cases of $\omega_{\perp}=\bigl{\{}0,-0.1,-0.3,\omega^{{\Pi}_{i}=0}_{\perp}\bigl{\}}c_{s}/R_{0}$ (red), and the full NL scans shown in (a) and (b) for $\hat{s}=0.4$ (green) and $\hat{s}=0.8$ (blue), respectively. All results are for $q=2.05$, $R_{0}/L_{Ti}=10.96$. The error bars (calculated by a rolling average on the time traces) are shown for some example cases.

Next, we combine LMD with a source of intrinsic momentum flux to drive strong flow shear. To predict the self-consistent rotation gradient that could be achieved in experiments, we perform NL gyrokinetic simulations with adiabatic electrons for tilted elliptical magnetic flux surfaces with $\epsilon=0.36$, elongation $\kappa=1.5$, and elongation tilt angle $\theta_{\kappa}=\pi/8$. This tilt angle was identified to generate the strongest intrinsic momentum flux in a CBC-like equilibrium [42].

Without external momentum sources, steady state operation of a tokamak requires $\hat{\Pi}_{i}=0$. To enforce this for a given equilibrium in local gyrokinetic simulations, we need to scan over the flow shear $\omega_{\perp}$ to find the value at which $\hat{\Pi}_{i}=0$. This is illustrated in Fig. 2 (a), where we see that as we increase the magnitude of the (negatively valued) flow shear, $\hat{\Pi}_{i}$ gradually drops to zero and then changes sign. The flow shear for which $\hat{\Pi}_{i}=0$ is the expected steady-state value in experiment. We see that a higher $\hat{s}$ results in a higher magnitude of the steady-state flow shear $|\omega_{\perp}|$. This is both because higher $\hat{s}$ lowers the Prandtl number (reflected in Fig. 1 (d)) and because higher $\hat{s}$ increases the intrinsic momentum flux (reflected by values of $\hat{\Pi}_{i}$ at $\omega_{\perp}=0$ in Fig. 2 (a)) [40, 42]. As expected, the heat flux is reduced by a larger $|\omega_{\perp}|$ [54] as shown in Fig. 2 (b).

Making a 4D NL scan of $(q,\hat{s},R_{0}/L_{Ti},\omega_{\perp})$ is computationally intensive. To reduce the numerical costs, we use a QuasiLinear (QL) model developed specifically to estimate $\hat{\Pi}_{i}/\hat{Q}_{i}$ [65]. This allows us to scan $\omega_{\perp}$ and find the steady-state flow shear $\omega^{{\Pi}_{i}=0}_{\perp}$ such that $\hat{\Pi}_{i}(\omega^{\Pi_{i}=0}_{\perp})=0$. We then compare two NL simulations, one with $\omega_{\perp}=\omega^{\Pi_{i}=0}_{\perp}$ and one with $\omega_{\perp}=0$, to determine the heat flux reduction by intrinsic rotation. A benchmark of the flux ratio $\hat{\Pi}_{i}/\hat{Q}_{i}$ between NL simulations and our QL model is shown in Fig. 2 (c) and (d). We see that our QL model predicts $\omega^{{\Pi}_{i}=0}_{\perp}$ well. Note that we use a linear interpolation between neighboring data points. To estimate the Prandtl number for up-down asymmetric cases, we perform two NL simulations at $\omega_{\perp}=\bigl{\{}-0.1,-0.3\bigl{\}}c_{s}/R_{0}$. As explained in Ref. [47], we adjust Eq. (3) to be

$$\overline{\text{Pr}}_{i}=\Delta\left(\frac{\hat{\Pi}_{i}}{\hat{Q}_{i}}\right)\frac{R_{0}}{L_{Ti}}\frac{\epsilon}{q\Delta\omega_{\perp}}\frac{c_{s}}{R_{0}},$$ (4)

so that the Prandtl number is proportional to the slope in Fig. 2 (c) and (d), where $\Delta$ indicates the difference of quantities between their values at $\omega_{\perp}=-0.1c_{s}/R_{0}$ and $-0.3c_{s}/R_{0}$. Hence, four NL simulations are performed for each $(q,\hat{s},R_{0}/L_{Ti})$, i.e. $\omega_{\perp}=\bigl{\{}0,-0.1,-0.3,\omega^{{\Pi}_{i}=0}_{\perp}\bigl{\}}c_{s}/R_{0}$ (shown in Fig. 2 (c) and (d) by the red lines). Note that the QL model is only used to obtain $\omega^{{\Pi}_{i}=0}_{\perp}$.

Figure 3 summarizes the effectiveness of up-down asymmetry combined with the LMD regime to create flow shear that stabilizes turbulence. When going from subplot (a), (b), (c) to (d), one gradually moves into the LMD regime. As we can see, the intrinsic flow shear strengthens, together with the decrease in Prandtl number, enabling a reduction in the heat flux by the intrinsic flow shear. Outside the LMD regime ($\overline{\text{Pr}}_{i}\gtrsim 0.5$), illustrated by Figs. 3 (a) and (b), only weak flow shear ($|\omega_{\perp}|\lesssim 0.2c_{s}/R_{0}$) is produced and the heat flux is barely affected. In the LMD regime ($\overline{\text{Pr}}_{i}\lesssim 0.5$), shown in Fig. 3 (c) and (d), the Prandtl number is low, enabling a strong intrinsic flow shear ($|\omega_{\perp}|\gtrsim 0.3c_{s}/R_{0}$) that significantly reduces the heat flux and can increase the critical gradient by up to $25\%$. Therefore, a combination of LMD with up-down asymmetry can create strong stabilizing intrinsic flow shear in the tokamak core.

Figure 3: The heat flux $\hat{Q}_{i}$ for the self-consistent value of flow shear $\omega^{\Pi_{i}=0}_{\perp}$ (ensuring $\hat{\Pi}_{i}=0$) driven by up-down asymmetric flux surfaces (blue), compared to the same cases without any flow shear $\omega_{\perp}=0$ (red) as a function of $R_{0}/L_{Ti}$ for (a) $q=2.05,\hat{s}=0.1$, (b) $q=2.05,\hat{s}=0.4$, (c) $q=2.05,\hat{s}=0.8$, and (d) $q=1.05,\hat{s}=0.8$. The Prandtl numbers (black) are also shown for each equilibrium. The flow shear values $\omega^{\Pi_{i}=0}_{\perp}$ (unit $c_{s}/R_{0}$) required to achieve $\hat{\Pi}_{i}=0$ are indicated by the numbers neighboring each blue data point.

Figure 4: Subplots (a) and (b) show the radial profiles of ion temperature and safety factor measured for MAST shot 24600 at $t=0.28s$. The profiles are obtained from TRANSP interpretative analysis, using the data processing code Pyrokinetics [67]. The three vertical dashed lines indicate the three considered radial locations $\rho_{pol}=\left\{0.71,0.77,0.84\right\}$. The inset in (b) shows the MAST experimental flux surface shape at these three radial locations as well as the artificially tilted elliptical geometry. (c) shows the experimentally measured toroidal flow shear $|\omega_{\perp}|$ produced by NBI (blue line) and the GENE predictions for intrinsic flow shear (stars) that could be created intrinsically if MAST could create the tilted elliptical geometries. (d) shows the heat flux with flow shear using the tilted elliptical geometry, where the vertical dashed lines indicate the predicted steady-state intrinsic flow shear.

To check the experimental feasibility of this strategy, we take an experimental equilibrium from the Mega Ampere Spherical Tokamak (MAST). We consider shot 24600 at time $t=0.28s$ (with an on-axis NBI injection), which has the ion temperature profile shown in Fig. 4 (a) and the flux surfaces and $q$ profile shown in (b). The experimentally measured kinetic profiles are averaged over $0.27s\leq t\leq 0.29s$ to reduce the numerical uncertainties. We use $\rho_{pol}=\psi^{1/2}_{n}$ as the radial coordinate, where $\psi_{n}$ is the normalized poloidal flux. As we can see, the $q$ profile is low over a wide radial range, which, together with the tight aspect ratio of MAST, facilitates LMD. No significant MHD instabilities are observed at around $t=0.28s$ [68] and the system is in a quasi-steady state. We perform gyrokinetic simulations at three radial locations $\rho_{pol}=\left\{0.71,0.77,0.84\right\}$, which are dominated by ITG driven turbulence. We tested other radial locations, but they are dominated by either Micro-Tearing Modes (MTM) or Electron Temperature Gradient (ETG) modes, which are computationally much more challenging for NL gyrokinetic simulations. Kinetic electrons are included in these simulations. Our simulations indicate that these three radial locations are close to or in the LMD regime, with $\text{Pr}_{i}=\left\{0.58,0.40,0.41\right\}$.

To verify that one can expect the LMD regime to reduce turbulence, we artificially modified the MAST flux surfaces to be tilted ellipses as shown in Fig. 4 (b) without considering the external torque from NBI, and performed NL flow shear scans to find $\omega^{\Pi_{i}=0}_{\perp}$. These flow shear values are shown in Fig. 4 (c), alongside the experimental flow shear profile. This illustrates that experimentally significant flow shear levels driven by external NBI sources could actually be driven intrinsically if MAST could achieve tilted elliptical flux surfaces. Figure 4 (d) shows that this intrinsic flow shear almost entirely quenches the ITG turbulence and thus significantly reduces the heat flux.

While it is difficult to create strongly up-down asymmetric flux surfaces in MAST, a new spherical tokamak under construction, SMART [69, 70, 71, 72, 73], has stronger shaping capabilities. Figure 5 (a) shows a preliminary up-down asymmetric geometry that can potentially be achieved in SMART with the addition of its advanced coils. We consider one radial location at $\rho_{pol}=0.79$, dominated by ITG turbulence, and perform a similar NL flow shear scan as we did for MAST. The important parameters are $R_{0}/L_{Ti}=6.32$, $q=1.34$, $\hat{s}=1.25$, $\epsilon=0.39$. A strong intrinsic flow shear of $|\omega^{\Pi_{i}=0}_{\perp}|=0.23c_{s}/R_{0}$, together with one order of magnitude heat flux reduction (shown in Fig. 5 (b) and (c), respectively) are observed, further demonstrating the experimental significance of our method.

Figure 5: Subplot (a) shows an up-down asymmetric equilibrium that can be achieved by the SMART tokamak with its advanced coil set. Subplots (b) and (c) show the ratio $\hat{\Pi}_{i}/\hat{Q}_{i}$ and $\hat{Q}_{i}$ as functions of flow shear from NL GENE simulations at $\rho_{pol}=0.79$.

Conclusions–In this letter, we propose and analyze a new method to drive strong intrinsic rotation shear in tokamaks: combining the LMD regime with up-down asymmetry. Using a large number of NL and QL gyrokinetic simulations, we self-consistently calculate the flow shear and show that it can significantly stabilize ITG turbulence. The experimental feasibility of this method is addressed by considering a MAST experimental equilibrium and a planned SMART equilibrium. Moreover, our results actually present a lower bound for the impact of the intrinsic flow shear as we neglected the pinch term in our analysis [47] and any external momentum drive. Based on our preliminary study, SMART should be able to test the main predictions in this Letter. This novel approach for generating strong flow shear does not rely on external momentum injection and therefore could be an attractive way of creating strong rotation in future large spherical tokamaks like STEP [74].

Acknowledgement–The simulations in this work are performed on CSCS Daint, Cineca Marconi and CSCS Lumi clusters, with a total computational cost of about $3\times 10^{7}$ CPU hours. The authors thank Prof. Ben McMillan, Dr. Antoine Hoffmann, Dr. Francis Casson, Arnas Volcokas and Alessandro Balestri for the fruitful discussions. This work has been carried out within the framework of the EUROfusion Consortium, partially funded by the European Union via the Euratom Research and Training Programme (Grant Agreement No. 101052200 - EUROfusion). The Swiss contribution to this work has been funded by the Swiss State Secretariat for Education, Research and Innovation (SERI). Views and opinions expressed are however those of the author(s) only and do not necessarily reflect those of the European Union, the European Commission or SERI. Neither the European Union nor the European Commission nor SERI can be held responsible for them. This work was supported in part by the Swiss National Science Foundation and by the EPSRC Energy Programme (Grant Number EP/W006839/1).