# Bipolar Thermoelectricity in Bilayer-Graphene–Superconductor Tunnel
Junctions

L. Bernazzani Corresponding author<EMAIL_ADDRESS>Dipartimento di Fisica ”E. Fermi”, Università di Pisa, I-56127 Pisa, Italy
Fachbereich Physik, Universität Konstanz, D-78457 Konstanz, Germany  G.
Marchegiani  Quantum Research Center, Technology Innovation Institute, Abu
Dhabi, P.O. Box 9639, United Arab Emirates  F. Giazotto  NEST, Istituto di
Nanoscienze CNR & Scuola Normale Superiore, I-56127 Pisa, Italy  S. Roddaro
Dipartimento di Fisica ”E. Fermi”, Università di Pisa, I-56127 Pisa, Italy
NEST, Istituto di Nanoscienze CNR & Scuola Normale Superiore, I-56127 Pisa,
Italy  A. Braggio https://sites.google.com/site/alessandrobraggio/Home NEST,
Istituto di Nanoscienze CNR & Scuola Normale Superiore, I-56127 Pisa, Italy.

###### Abstract

We investigate the thermoelectric properties of a hybrid nanodevice composed
by a 2D carbon-based material and a superconductor. This system presents
nonlinear bipolar thermoelectricity as induced by the spontaneous breaking of
the Particle-Hole (PH) symmetry in a tunnel junction between a BiLayer
Graphene (BLG) and a Bardeen-Cooper-Schrieffer (BCS) superconductor. In this
scheme, the nonlinear thermoelectric effect, predicted and observed in SIS′
junctions, is not affected by the competitive effect of the Josephson
coupling. From a fundamental perspective, the most intriguing feature of this
effect is its bipolarity. The capability to open and control the BLG gap
guarantees improved thermoelectric performances, that reach up to $1\,$mV/K
regarding the Seebeck coefficient and a power density of $1~{}{\rm nW}/\mu{\rm
m}^{2}$ for temperature gradients of tens of Kelvins. Furthermore, the
externally controlled gating can also dope the BLG, which is otherwise
intrinsically PH symmetric, giving us the opportunity to investigate the
bipolar thermoelectricity even in presence of a controlled suppression of the
PH symmetry. The predicted robustness of this system could foster further
experimental investigations and applications in the near future, thanks to the
available techniques of nano-fabrication.

##  I Introduction

Thermoelectricity is conventionally investigated in materials where there is a
well-defined dominant carrier [1, 2, 3]. Indeed, the Boltzmann theory of
transport demonstrates that in the linear response regime the Seebeck
coefficient has the same sign of the charge of the dominant carrier [1, 4]. In
other words, thermoelectricity seems to be present only in systems where the
PH symmetry is broken to a certain extent [5]. This becomes paradigmatic, in
superconducting systems, where it should be expected that PH symmetry and the
dissipationless flow of Cooper pairs would implicitly limit any thermoelectric
effect. Notwithstanding, a few works [6, 7] anticipated that even
superconducting systems could exhibit a thermoactive behaviour. Along this
route, theoretical [8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18] and experimental
[19, 20, 21] research has shown that nonlocality or phase coherence may still
trigger intriguing thermoelectric phenomena in superconducting and hybrid
systems. Furthermore, some authors even reported the existence of Absolute
Negative Conductance (ANC) in superconducting tunnel junctions in out-of-
equilibrium conditions or under microwave irradiation [22, 23, 24, 25].

More recently, instead, in a number of theoretical papers [26, 27, 28, 29],
the spontaneous breaking of the PH symmetry was predicted in Superconductor-
Insulator-Superconductor’ (SIS′) tunnel junctions subject to a finite thermal
gradient. In these works, a strong bipolar thermoelectric effect was clearly
identified in the nonlinear regime of the thermal bias. These latter
predictions were confirmed experimentally [30, 31] providing also with methods
to make superconducting non-volatile memories [26, 32] or other spin-active
devices [33].

However, there are several drawbacks concerning the physical realisation of
the cited effect in asymmetric SIS′ junctions. The most serious issue is given
by the Josephson coupling, which needs to be substantially suppressed [30], as
theoretically anticipated in Ref. 28. Indeed the Josephson effect may
completely kill the thermoelectric generation of power by shunting the
junction at low biases, i.e. $|eV|\ll(\Delta_{R}+\Delta_{L})$, here
$\Delta_{R/L}$ are the values of the bulk superconducting gap function on the
two sides of the SIS′ junction, therefore restoring the PH symmetry and
suppressing the thermoelectric effect. This limitation may be even more
serious in out-of-equilibrium conditions since the Josephson coupling presents
anomalous jumps when thermal gradients are applied to asymmetric SIS′
superconducting junction [34, 35]. Hence, the possibility of making a junction
using instead a semiconductor coupled to superconductors would be beneficial
in this respect. Furthermore, this nonlinear thermoelectric effect manifests
strongly at the so-called ”matching peak” conditions
$V_{peak}=\pm|(\Delta_{R}-\Delta_{L})/e|$, where the singularities in the BCS
DoS are aligned. So the possibility to engineer in detail the semiconducting
gaps would be beneficial for the performance of the device.

Nevertheless, the aforementioned bipolar thermoelectric effect requires that
the semiconductors forming the junction have a gap comparable with those of
the superconductors which they are coupled to, i.e. they should have gaps of
the order of few meV [36, 37]. On the other hand, elemental or binary compound
semiconductors have gaps larger than hundreds of meV [36]. This fact pushed us
to consider a completely different platform, i.e. carbon-based low dimensional
materials, that can exhibit smaller bandgaps, whose value can even be tuned in
some cases by electrostatic gating [38].

Hereafter we demonstrate the bipolar thermoelectricity in hybrid junctions
composed of BiLayer Graphene (BLG) coupled to BCS superconductors (S) through
an insulating barrier (I). Notably, the bandgap of BLG is controllable by
electrical gating [39], making possible to adapt the bandgap to the best
operating conditions, which are connected with the superconducting gap and
with the temperature of the hot lead (BLG), by changing just the electric
gates voltages. This peculiar phenomenon increases the flexibility and the
performances of the device. For instance, it allows one to switch on and off
thermoelectricity in the junction. Moreover, the replacement of BLG in place
of one superconducting electrode relaxes also the constraints on the hot
electrode temperature, leading to an even wider applicability of the
thermoactive elements, which in principle could operate up to hundred kelvins.

##  II Device description

(a)

(b)

Figure 1:  a) Schematics of the proposed device. The BiLayer Graphene (BLG) is
placed on top of a Silica (SiO2) dielectric ($\varepsilon_{b}=3.9$) 100-nm-
thick, beneath which there is the bottom gate (doped Si). BLG is encapsulated
in a hBN monolayer but is Ohmic coupled with a Au electrode on the left side.
The Niobium (Nb) superconducting electrode is deposited partially on top of
the BLG. hBN constitutes the tunnelling barrier with the superconducting
electrode. Finally, an oxide ($\varepsilon_{t}=7.5$)[40], approximately
190-nm-thick, is deposited above the electrodes and a top-gate made from Au is
placed on top of this dielectric. b) Band scheme of the BLG [41, 42, 43]. A
gap has been opened using the top/bottom gates. With the dashed lines we also
show how the BLG Fermi energy $E_{F}^{BLG}$ moves along the bands, due to an
asymmetric gating effect. We assumed $v_{F}\approx 10^{6}\,\rm m/s$ as the
Fermi velocity in graphene.

Low dimensional systems have attracted a lot of interest recently [44] and
intriguingly have been already discussed as promising materials for improving
thermoelectric performance in the pioneering work of Ref. 45. BLG, in
particular, is a two-dimensional material formed by two planes of graphene
stacked onto each other. It naturally inherits the intrinsic PH symmetry of
graphene if it is not extrinsecally doped. However its bands deviate from the
conventional linear dispersion relations of graphene [46, 47]. In fact, by
tight-binding computations [48, 49] one can show that they are characterised
by a nonlinear dispersion relation [50]. Similarly to monolayer graphene, BLG
is a semimetal, and, in the pristine form, valence and conduction bands touch
each other at the high symmetry K, K’ points. Notably, by applying an electric
field perpendicular to the BLG plane, a gapped phase can be induced,
effectively behaving as an electrically controllable semiconductor. In Fig. 1a
we display a sketch of the setup discussed below, where the BLG is coupled to
a top and a bottom gate that may be tuned externally. The gate-controlled
resultant bandgap is highlighted in Fig. 1b. In the absence of gating, the BLG
is gapless and the dispersion of the lowest-energy band is approximately
quadratic $\epsilon(k)=\pm(\sqrt{1+4(v_{F}\hbar
k)^{2}}-1)/2\approx\pm(v_{F}\hbar k)^{2}$ (orange solid curves) [41], with
$v_{F}$ the Fermi velocity. BLG’s gap can be opened and controlled by applying
a differential gating between the backgate $V_{g,B}$ and the topgate $V_{g,T}$
[41, 42] (solid red curves). However, gating may also induce charge doping at
the BLG planes, spoiling the perfect PH symmetry and paving the way to the
intriguing possibility of investigating a bipolar thermoelectric device in the
presence of a controlled breaking of the PH symmetry as well. As reported in
[42] this tunable bandgap can be computed using an approximate formula,
according to which it is given by [41, 42],

$2E_{g}=\frac{|U|t_{\perp}}{\sqrt{U^{2}+t_{\perp}^{2}}}\,,$ (1)

where $U$ is the potential energy difference between the two graphene planes,
also called the gating asymmetry parameter, and $t_{\perp}\approx 0.4\,$eV is
the inter-layer hopping energy term in the tight binding approximation of
BLG’s electronic bandstructure. The approximate energy dispersion in the
gapped state [41] $\epsilon(k)\approx\pm\sqrt{U^{2}/4+\pm(v_{F}\hbar
k)^{4}/t_{\perp}^{2}}$ is shown in Fig. 1b (solid red curve).

This paper deals with the theoretical study of a BLG-I-S tunnel junction,
whose schematic representation is given in Fig. 1a. The BLG is on the left
side of the junction and is deposited on the Silica substrate. We suppose to
be able to heat up the BLG electronic temperature in order to induce an
effective thermal gradient across the junction by means of irradiation by a
photon source or by the Joule effect due to normal currents flowing between
electrodes [51, 52] (this can be usually accomplished using other normal metal
electrodes, e.g. yellow contact in Fig. 1a). The BLG is in clean contact with
a gold (Au) electrode (see Fig. 3a), so that the voltage drop $V$ applied
across the structure localizes at the BLG-I-S junction. In other words, the Au
electrode and the BLG are equipotentials. The electrostatic gating is applied
with the aid of a top-gate $V_{g,T}$ capacitively-coupled to BLG through an
oxide and a $Si^{+}$ layer, with the role of the bottom gate, beneath the
Silica layer, as shown in Fig. 1a [40]111The top/bottom gates are referred to
the BLG Ohmic contacts in order to induce a proper gating on the BLG which is
not affected by the bias applied to the junction.. The right electrode
composing the junction is a Bardeen-Cooper-Schrieffer (BCS) superconductor
(grey in Fig. 1a). In order to have a superconducting gap commensurate with
those of BLG and, at the same time, a good range of operating temperatures, we
chose Nb for the BCS superconductor. This metal has a critical temperature of
$T_{c}\approx 9\,\rm K$. It is deposited partially on the encapsulated BLG in
order to provide a tunnel contact with it. The encapsulation of BLG is made
from hexagonal Boron Nitride (hBN) monolayer. This material constitutes proper
tunnel barriers since it is a good insulator with an almost perfect
crystallographic order [54].

##  III Bipolar thermoelectricity with PH symmetry

In the tunnelling framework, we can compute the charge and heat currents
according to the following formula [55, 37, 56, 26, 27]

$\displaystyle\begin{bmatrix}I\\\
\dot{Q}_{L}\end{bmatrix}=\frac{4\pi|\mathcal{T}|^{2}}{\hbar}\int_{-\infty}^{+\infty}d\varepsilon$
$\displaystyle\begin{bmatrix}-e\\\
\varepsilon\end{bmatrix}N_{BLG}(\varepsilon)N_{S}(\varepsilon+eV)\cdot$
$\displaystyle\cdot[f_{L}(\varepsilon)-f_{R}(\varepsilon+eV)]$ (2)

where $|\mathcal{T}|$ is the transmittivity of the tunnelling barrier, which
in our case is given by the hBN layer between the BLG and the superconducting
electrode. In writing Eq. (III), we further assumed that the barrier
transmissivity is independent of the gate voltages $V_{g,T}$ and $V_{g,B}$
applied to the BLG. This approximation is reasonable at the energy scales
considered in this work, where the induced BLG bandgaps and the applied
voltage bias are the order of few meV, but may fail for a very strong gating.
The electrons’ gas in the semi- and superconductive leads follows the Fermi-
Dirac distributions
$f_{i}(\epsilon,T_{i})=[1+\exp{(\epsilon/k_{B}T_{i})}]^{-1}$, where we assume
that carriers are in quasi-equilibrium, being well described by a standard
quasiparticle distribution with an electronic temperature $T_{i}$ [57]. We
introduce also the tunnelling barrier bias voltage $V$.

In Eq. III the Density of States (DoS) of the two electrodes are
$N_{BLG/S}(E)$. In our case the right electrode is a BCS superconductor,
therefore

$N_{S}(\varepsilon)=N_{F}\,\Biggl{|}\Re{\Biggl{[}\frac{\varepsilon+i\Gamma_{D}}{\sqrt{(\varepsilon+i\Gamma_{D})^{2}-\Delta_{R}(T_{R})^{2}}}\Biggr{]}\Biggr{|}}$
(3)

with $N_{F}$ the DoS of the metal in the normal state at the Fermi energy,
$\Delta_{R}(T_{R})$ is the superconducting gap and $\Gamma_{D}$ is the Dynes’
phenomenological correction to BCS DoS [58]. On the other hand, the left
electrode is composed by the BLG whose DoS for $E\ll t$ can be approximated
with the following expression [41, 59]:

$N_{BLG}(\epsilon)\approx\frac{t_{\perp}}{\sqrt{3}\,\pi\,t^{2}}\cdot\frac{|E|\,\Theta(|E|-E_{g})}{\sqrt{{E}^{2}-E_{g}^{2}}}$
(4)

where $\Theta(x)$ is the Heaviside step-function and $E=\epsilon+E_{F}^{BLG}$
with $E_{F}^{BLG}$ being the BLG’s Fermi energy. The tight-binding energies
for the in-plane hopping integral of the graphene is $t\approx 3\,$eV and the
out-of-plane hopping term is $t_{\perp}\approx 0.4\,$eV. We assume the
insulating barrier between BLG and superconducting electrodes being the hBN
222The hBN is commonly used to encapsulate graphene since its crystalline
configuration matches well with that of graphene diminishing the surface
disorder at the interfaces.. We also assume to have a tunnelling conductance
$G_{T}\approx(10\,{\rm k}\Omega)^{-1}$ 333The conductance
$G_{T}=4N_{F}e^{2}|\mathcal{T}|^{2}t_{\perp}/(\hbar\sqrt{3}t^{2})$ would
represent the differential conductance of the tunnel junction when the
superconductor is in the normal state and the BLG is gapless. for the hBN
separating the BLG electrode, above which the superconducting electrodes will
be deposited forming a junction of about $0.1\,\mu\rm{m}^{2}$, this choice is
consistent with values reported in Ref. 62. Despite a larger active surface
would lead to a higher conductance that would be beneficial for the
performance of the proposed junction, it is nonetheless important that the
barrier between the leads is a good insulating barrier, since we are going to
establish a sizeable temperature gradient across the junction. Furthermore,
enlargement of the active surface could also lead to unwanted effects such as
pinhole points [63].

Clearly, either asymmetric gating or charge traps at the oxide surface may
also induce a deviation in the midgap position of the BLG’s Fermi energy
$E_{F}^{BLG}$ partially breaking the PH symmetry. However, with two
independent top and bottom gates it is always possible to restore the PH
symmetry in the BLG. Indeed, since the bandgap $2E_{g}\approx U$ (for $U\ll
t_{\perp}$, see Eq. (1)) is induced mainly by the differential gating mode,
the common mode can be exploited to modify $E_{F}^{BLG}$ changing the BLG
doping [41, 42, 64, 65, 48, 40]. More specifically, we can find a line of
points in the $(V_{g,B},V_{g,T})$ plane for which the BLG turns out to be
fully PH symmetric and gapped. This is a straight line with a slope of
$-C_{b}/C_{t}$ where $C_{i}=\varepsilon_{0}\varepsilon_{i}/d_{i}$ are the
capacitances with top ($i=t$) or bottom ($i=b$) gates determined by the
thicknesses $d_{b}$($d_{t}$) oxide layers and their dielectric constant
$\varepsilon_{t}$ ($\varepsilon_{b}$) [40]. We remark that the charge
accumulated in the BLG through gating modifies the Fermi energy and so the
total junction differential conductance (see discussion in Sec. III). However,
these effects cannot be interpreted as a gating dependence of the
transmissivity of the tunnelling barrier.Therefore, in the first part of this
paper, we investigate the device in the gating conditions where the BLG
electrode is fully PH symmetric, i.e. $E_{F}^{BLG}=0$ according to our
notation. In this configuration, due to the PH symmetric nature of the
electrodes, one should expect that the junction could never display a linear
thermoelectric behaviour [5, 66]. Notwithstanding, when BLG’s electronic
temperature and the thermal gradients exceed finite thresholds, we predict the
possibility to induce spontaneously a thermoelectric regime in the junction
just by tuning BLG’s gap at values higher than Nb’s superconducting gap.
Furthermore, it is easy to check that the PH symmetry of the two leads, i.e.
$N_{i}(E)=N_{i}(-E)$, would result in the anti-symmetric nature of the current
with respect to the electric bias, i.e. $I(V)=-I(-V)$. This phenomenology is
similar to that reported in [26, 27, 30]. However, it is worth noting that the
possibility to tune the electrode gap through gating, is quite unique to this
system. In fact in the SIS′ case one should change the material composing the
electrode [30] or recurring to complex phase-controlled gaps in proximity
junctions [35].

The discussed bipolar thermoelectric effect is nonlinear since a finite
temperature difference is required to activate the junction. When
$E_{g}>\Delta_{R}(T_{R})$ the junction can be in the thermoactive regime,
i.e., displaying ANC $IV<0$, for an appropriate temperature difference between
the two terminals. In Fig. 2 we show the IV characteristics in the subgap
regime when the junction is thermoactive. The electronic temperature of the
BLG is assumed to be kept at a high temperature (e.g. $T_{L}=50\,$K in Fig.
2a) and the superconductor will be on the cold side with $T_{R}=T_{bath}=1\,$K
444In principle inverting the ratio between $E_{g}/\Delta_{0,R}$ also the
opposite temperature gradient may be assumed. But this is less practical since
the temperature of the superconducting lead should not exceed its critical
temperature in any case.. Similarly to the SIS′ case [27], we can identify two
different thermoelectric regimes that depend exclusively on the absolute
module $|V|$. These are the linear-in-bias regime, that is revealed by the
negative differential conductance at $V\approx 0$, and the nonlinear-in-bias
effect which manifests strongly with the biggest ANC at the matching peak
condition $|eV_{p}|=E_{g}(U)-\Delta_{R}(T_{bath})$.

(a)

(b)

Figure 2:  IV characteristics of the proposed junction in the thermoelectric
regime $E_{g}>\Delta_{R}(T_{R})$. The thermoelectricity is easily spotted
because $IV<0$ and this implies ANC. a) IV characteristics plotted by varying
the BLG’s gap (that is obtained by varying the gating voltages while keeping
the BLG in the PH symmetry line). Here the BLG temperature is kept at
$T_{L}=50\,$K b) IV characteristics plotted varying the electronic temperature
of the BLG while keeping the BLG gap fixed at $2E_{g}=8\,$meV. The height of
the peak, and consequently the strength of the thermoelectricity, grows by
increasing the hot lead temperature. Quite peculiarly, due to the nonlinearity
of the effect, it appears that the Seebeck voltage $V_{S}$, defined by
$I(V_{S})=0$, depends weakly on the BLG temperature being more connected to
the position of the matching peak $V_{p}$ instead, that in this system is
unaffected by the $T_{L}$. In the computation we assumed the superconductor to
reside at the bath temperature of $T_{R}=T_{bath}=1\,$K, a Dynes parameter
$\Gamma_{D}/\Delta=10^{-4}$, and tunneling conductance $G_{T}=(10~{}{\rm
k}\Omega)^{-1}$.

By looking at the peaks in Fig. 2a,b, we observe how their position is changed
following the evolution of the BLG gap $E_{g}$ that is modified by external
gating. Notably, we see that the height of the ANC peaks is weakly dependent
on the BLG gap (see 2a) but is more affected by the temperature of BLG (see
2b). Finally, we note that in this configuration where the BLG is hot and the
superconductor is cold, the hot temperature does not affect the position of
the peak since the $E_{g}$ depends only on the electrostatic potential $U$.
This is a substantial difference with respect to the SIS′ case [26, 27, 30].
As a consequence, it is also important that the thermoelectric power $P=-IV$
is a monotonously increasing function of the hot temperature in our system in
contrast to the SIS′ concept where the closing of the gap with temperature
limits the thermoelectricity to a well-defined range, spanning at most an
interval of few Kelvins.

Since we are not dealing with a linear thermoelectric effect, we cannot rely
on standard figures of merit typical for thermoelectric elements, such as ZT
[5, 68]. Hereafter we will focus on alternative figures of merit in order to
quantify the performances of the BLG-I-S thermoelectric junction in comparison
with other recent proposals [30]. In particular, it is interesting to show
them as in Fig. 3 by using contour plots. Therefore investigating how they
depend on both the BLG gap $E_{g}$ and temperature $T_{L}$. In Fig. 3a we show
the maximum thermocurrent $I_{M\\!A\\!X}=|I(\pm V_{p})|$, which is obtained at
the matching peak voltage $V_{p}$ and, correspondingly, in Fig. 3b the
thermoelectric power delivered at the same point $P(V_{p})=-I(V_{p})\cdot
V_{p}$ by the junction 555Even if the absolute value of those quantity depends
also on the Dynes’ parameter $\Gamma_{D}$ the general behaviour would be quite
independent on it.. The grey areas in the plots denote the regions with
$E_{g}<\Delta_{R}(T_{R})$, where the system is purely dissipative, thus
producing no net thermoelectric current. Both of the quantities are
monotonically increasing functions of the BLG temperature $T_{L}$ and quickly
saturate for $T_{L}\gtrsim E_{g}/k_{B}$ (not shown). These figures of merit
behave not monotonically with respect to the BLG gap $E_{g}$, displaying a
maximum for a temperature-dependent optimal value. We notice that the values
of the thermoelectric charge current are quite high for a junction of
$0.1\,\mu\rm m^{2}\,$. At the same time the output power, shown in figure 3b,
is notably high, since it reaches values above $750\,$pW for high-temperature
gradients and intermediate gap values which correspond to a power density of
roughly $1\,\rm{nW}/\mu\rm{m}^{2}$. Since it is computed at the matching peak,
it is natural to ask how this quantity would be affected by the quality of the
superconducting contact, leading to a higher Dynes’ parameter. We report,
without showing it, that the output power of the junctions is only reduced by
a logarithmic correction of about a factor of 2 when changing the Dynes’
parameter by two orders of magnitudes (from the $10^{-4}\Delta_{0}$ used here
to $10^{-2}\Delta_{0}$).

(a)

(b)

(c)

(d)

Figure 3:  Thermoelectric figures of merits for the proposed junction at the
PH symmetric point as a function of the BLG bandgap $2E_{g}$ and the BLG
electronic temperature $T_{L}$. a) Displays the maximum current in nA obtained
at the matching peak condition $I(V_{p})$. b) Shows the thermoelectric power
at the matching peak $-I(V_{p})V_{p}$ in units of pW. c) The thermodynamic
absolute efficiency computed at the matching peak $V_{p}$, where the output
power is roughly maximum. d) Displays the nonlinear Seebeck coefficient in
units of $\mu\rm V/K$. The remaining parameters are the same of Fig. 2.

Furthermore, it is interesting to discuss, looking at Fig. 3c, the absolute
efficiency of the thermoelectric engine at the maximum power. An useful
approximation for such quantity is the value it assumes at $V=V_{p}$, i.e.
$\eta_{p}\approx-I(V_{p})V_{p}/\dot{Q}_{L}(V_{p})$, since the junction
delivers maximum power approximately at the matching peak point. At the same
time, the heat losses in the junction is probably dominated by the
quasiparticle contribution $\dot{Q}_{L}(V_{p})$ computed using Eq. III
666Obviously the heat may potentially be transferred also by other
contributions such as phonon or even radiative contributions that are more
difficult to be estimated. Our computation of efficiency represents a sort of
optimal electronic efficiency of the device.. In this case the best
efficiencies are obtained at the lowest BLG temperatures $T_{L}$, where the
decrease in the heat losses is much bigger than the decrease of the
thermopower. Notably, the absolute efficiency reaches values up to $40\%$ and,
when the hot temperature is fixed, there is an optimal value for the gap. It
is worth noting that the efficiency here is the absolute thermodynamic
efficiency of the thermodynamic machine.

Finally, in Fig. 3d the thermoelectric performance, is shown in terms of
Nonlinear Seebeck Coefficient (NSC) which is computed as

$S=\frac{V_{S}}{T_{L}-T_{R}}$ (5)

[26], where the $V_{S}$ is computed numerically by finding the point at which
the IV characteristic crosses the zero current axis, i.e. $I(V_{S})=0$. This
is a natural extension, in a nonlinear picture, of the Seebeck coefficient
defined as $S=\lim_{\Delta T\to 0}\frac{\Delta V}{\Delta T}\bigl{|}_{I=0}$. By
the way, that is the solution of the linear thermoelectricity equation
$I=G\,\Delta V-G\,S\,\Delta T=0$ [5]. This quantity can be computed only in
the region where the junction is thermoelectric ($IV<0$), while the
dissipative regions are coloured in grey. The values of the NSC demonstrate
the strength of the thermoelectricity generation in this junction, that
reaches values higher than $1000\,\mu\rm V/K$ at the lowest temperature
gradient explored of $\approx 9\,$K. This graph reveals also a quite general
behaviour of this quantity which increases with the BLG bandgap. Looking at
panel 3d we can say that the highest values of the NSC can be reached at low
values of the temperature gradient where the thermopower is quite low (see
Fig. 3b) confirming once again, even in this nonlinear situation, a sort of
trade-off between efficiency and power production.

##  IV Partial breaking of PH symmetry

(a)

(b)  (c)

Figure 4:  a)IV characteristics for the PH asymmetric state of BLG (orange and
violet lines) for $T_{L}=50\,$K, and $V_{g,T}=-3\,$V. The green line is the PH
symmetric case. We notice that these characteristics show a thermoelectric
behaviour since in the neighbourhoods of the peaks $P=-IV>0$. b) Current at
zero bias (conventional thermocurrent) for $T_{L}=50\,$K. Purple areas denote
regions where $I(0)\approx 0$, a necessary condition for PH symmetry. c) NSC
$S^{+}$ for $T_{L}=20\,$K. Within purple regions the junction is not
thermoelectric and the NSC cannot be defined. Parameters: $T_{R}=1~{}$K,
$\varepsilon_{b}=3.9$, $\varepsilon_{t}=7.5$, $d_{b}=100$ nm, $d_{t}=192$ nm,
$d=0.35$ nm.

Having treated the PH symmetric case in Sec. III, now we consider the more
general case where the BLG is not perfectly gated at the PH symmetric line,
i.e., the aforementioned perfect charge neutrality line of the gates. The loss
of PH symmetry on the BLG side, going on with the stated conventions, can be
seen as the case where the Fermi energy $E_{F}^{BLG}\neq 0$. Indeed, according
to the self-consistent model of screening, adapted from Refs. 41, 42, 71, the
electrostatic asymmetry parameter $U$ [and so the bandgap, see Eq. (1)] and
the BLG Fermi energy should be computed self-consistently from the value of
the electrostatic gating voltages. We state here just the final formula that
can be derived from the more general treatments of Refs. 42, 71, 72.

A recent experiment [72] shows the important role of localized states, in
particular for designs with silicon dioxide substrate. These localized states
contribute to the charge density of the BLG but not to the screening of
electrostatic asymmetry parameter $U$ 777The presence of localized charge
states induced by impurities also determines a modification of the gating
conditions anyway these contributions can be simply accounted by including
appropriate offsets to the gate voltages, so we neglected them [72].. The
fundamental parameters that will be needed in the IV characteristics are the
electrostatic asymmetry parameter $U$ and the Fermi energy $E_{F}^{BLG}$ that
are given by 888The bandgap self-consistent expression of Eq. (6) is derived
using Eqs. (65) and (74) of Ref. [42] in the limit with $\varepsilon_{r}\sim
1$, $\Lambda\sim 1$ and assuming no mobile charge density when the Fermi
energy is in the gap. Correspondingly the Fermi-energy expression Eq. (7)
coincides with the chemical potential of Eq. (S13) of the supplementary
information of Ref. [72] where have been neglected terms involving the inverse
of BLG capacitance $C_{BLG}^{-1}$ being $C_{BLG}\gg C_{b},C_{t}$.

$\displaystyle
U\approx\frac{e\,\bigl{(}C_{b}V_{g,B}-C_{t}V_{g,T}\bigr{)}}{2C_{BLG}\bigl{[}1-\ln(|U|/4t_{\perp})/2\bigr{]}}\,,$
(6) $\displaystyle E_{F}^{BLG}\approx
e\frac{\,C_{b}V_{g,B}+C_{t}V_{g,T}}{C_{b}+C_{t}}\,,$ (7)

where $C_{BLG}=\varepsilon_{0}\varepsilon_{r}/d$ is the BLG interlayer
capacitance with $d\approx 0.35\,$nm the distance of the two graphene layers
and $\varepsilon_{r}\sim 1$ 999This simplified model may not be fully self-
consistent when the system is gated well outside the gap, but corrections are
expected to slightly renormalize the gap value not affecting physics
discussed.. The solution of the first equation has to be found numerically for
every gate configuration $(V_{g,B},V_{g,T})$. The experimental realisation of
the device determines the gate capacitances $C_{i}$ with $i=t,b$ . In order to
simplify the gating around the PH symmetric point, it is convenient to
consider experimental setups where $C_{b}\approx C_{t}$. In principle, by
plugging the these quantities into Eq. 1 and Eq. 4 one can easily determine
the BLG’s gap and the BLG DoS.

In this situation of partial breaking of the PH symmetry, a conventional
thermoelectric effect adds to the nonlinear contribution that we have outlined
in Sec. III. However, as we discuss below, this conventional contribution does
not completely spoil the bipolar thermoelectricity phenomenology. This result
shows that the bipolar thermoelectric effect is also stable against slight
perturbations of the PH symmetry of the electrodes. The stability is evident
by looking at Fig. 4a, where we fix $V_{g,T}=-3\,V$ and slightly change the
bottom voltage $V_{g,B}$ moving around the PH symmetry point (green line). In
this graph, in fact, we see that the current at the matching peak of the IV
characteristics is still thermoelectric ($IV<0$) even if there is a little
departure from the PH symmetry as signalled by a finite linear thermocurrent
contribution, i.e. $I(V=0)\neq 0$. The IV characteristics, in this case, are
no longer reciprocal, i.e., $I(V)\neq-I(-V)$, and they do not cross the origin
of the axes. Nevertheless, we notice that the system still displays a
behaviour connected to the spontaneous breaking of the PH symmetry. Indeed, at
$V=0$, the differential conductance $dI/dV$ is negative for both the orange
(violet) line cases. This negative conductance determines an electrical
instability. When the system is thermoelectric $IV<0$ this device effectively
operates as a battery; when connecting resistive load in parallel, and the
junction is able to support two opposite signs for the thermocurrent through
the load for the _same_ thermal gradient. In fact, if the load resistance is
sufficiently high we can find two different solutions for the current
circulating in the load (just imposing current conservation in the circuit as
in [27]) with opposite sign. This feature demonstrates that the thermoelectric
effect is still bipolar, i.e., two opposite signs of the thermocurrent for a
given thermal gradient are possible [30]. The intrinsic instability at
$V\approx 0$ and the bipolar nature of the effect are consequences of the
spontaneous breaking of the PH symmetry and would be important in applications
such as current-controlled thermoelectrical memories [30, 32] or even sensors
[76].

It is interesting to investigate how the zero bias thermoelectric current
$I(0)$ behaves, fixing a finite temperature gradient, as a function of the
voltage of the top and bottom gates, i.e., in the $(V_{g,B},V_{g,T})$ plane.
In Fig. 4b we see that this quantity is zero at the PH symmetry line (purple
region). The bright yellow lines parallel to that are the lines at which the
BLG gap equal the superconductor gap, and so the matching peak is more or less
at $V\approx 0$. This interpretation explains why $I(0)$ at these yellow lines
is so enhanced. Indeed we notice that the $I(0)$ current can be quite high and
even become comparable to the maximal thermocurrent generated by the bipolar
thermoelectric effect discussed before. This thermocurrent can be interpreted
as that from a conventional thermoelectric effect but presenting also a strong
nonlinear enhancement. In other words, we are discussing a thermoelectric
behaviour which is linear in the bias $V$ but nonlinear in the thermal
gradient.

Clearly, When the PH symmetry is broken the IV characteristics is no more
anti-symmetric and the bias values $V_{S}^{\pm}$, for which
$I(V_{S}^{\pm})=0$, are different depending on the sign of $V$, i.e.
$|V_{S}^{+}|\neq|V_{S}^{-}|$. As a consequence, we need to define a NSC which
is different depending on the positive/negative voltage branch such that
$S^{\pm}=V_{S}^{\pm}/(T_{L}-T_{R})$. Considering the NSC $S^{+}$, for
instance, we see that it is zero except for the green-yellow band in Fig. 4c.
The displayed behavior can be understood as follows. In the region of interest
$V_{g,B}\approx-V_{g,T}\approx 3$V, the bandgap is approximately constant,
reading $E_{g}\approx U/2\approx 6~{}$meV, as determined by the self-
consistent solution of Eq. (6). Hence, the overall behaviour of $S^{+}$ is
mainly determined by the evolution of the Fermi energy [see Eq. (7))], which
varies significantly once away from the PH symmetry line (black dotted in Fig.
4). This feature makes the green-yellow band asymmetric with respect to the PH
symmetry line. When considering $S^{+}$, we compute the zero-current state
with $V>0$. For $|V_{g,B}|<|V_{g,T}|$, the current-voltage characteristic
shifts towards higher current values (cf. 4a), so making the thermoelectric
peak not intersecting the $I=0$ axes above a threshold value of $E_{F}^{BLG}$.
Differently, for $|V_{g,B}|>|V_{g,T}|$, the characteristic is shifted in the
opposite direction, and the intersection persists until the matching peak
condition is at $V=0$. Since the value of the NSC mainly depends on $E_{g}$,
as discussed above, we notice that within the green and yellow region, $S_{+}$
depends only weakly on the gates (evolution determined by variation of the
Fermi energy). The maximum NSC reaches $300\,\mu\rm{V/K}$ for $T_{L}=20\,$K,
in agreement with Fig. 3d. Finally, we observe that these figures show that
small deviations from the PH symmetry would not completely spoil the bipolar
nonlinear thermoelectric effect, although one clearly loses the full
antisymmetry of the IV characteristics. This fact confirms the robustness of
the discussed bipolar thermoelectric effect, since it results to be stable
against small deviations from ideal gating conditions.

##  V Conclusions

In this article, we analysed the nonlinear thermoelectricity of a BLG-
Superconductor tunnelling junction. We demonstrated that by external gating it
is possible to tune the BLG gap and the thermoelectric effect by only
electrical means, achieving very promising values for the figures of merit,
e.g., NSC up to $1\,$mV/K at temperature of the BLG of $T_{L}=10\,$K. For such
temperature gradients, we predicted a thermoelectric generation of power of
hundreds of pW for $0.1\,\mu{\rm m}^{2}$ of active surface, which corresponds
to a surface power density of roughly $1\,{\rm nW}/\mu{\rm m}^{2}$ with a
maximal absolute thermodynamic efficiency up to $40\%$ for the electronic
degrees of freedom. Finally, we investigated how the different figures of
merit depend on the BLG electronic temperature and on the BLG gap in a
realistic range of variations of external parameters. The stated bipolar
thermoelectric phenomenology could survive even to a partial breaking of the
PH symmetry, as might be potentially induced by a non-ideal external
electrostatic gating of the BLG electrode. We discussed, for this more general
case as well, the thermoelectric properties of the device, showing indeed that
there is an interplay between the nonlinear bipolar thermoelectricity induced
by spontaneous breaking of the PH symmetry and the conventional unipolar
thermoelectricity as determined by a small PH breaking term. However, we
confirmed that this latter contribution can coexist with the former without
being detrimental for the bipolar thermoelectricity. We state one more time
that the phenomenology just outlined could lead to a wide versatility of this
sort of devices due to the tunability provided by the external gating
mechanism.

Finally, we hope that the sizable thermoelectric performances of this hybrid
device architecture, combined with its electrical flexibility, will trigger
the future experimental exploration of bipolar thermoelectricity at even
higher temperatures. Indeed, the tunability of the gap makes the device
adaptable to respond to unwanted changes in the temperature of the hot lead
and, at the same time, can be exploited to switch on and off its
thermoelectric behaviour. We have seen that hybrid superconducting-graphene
devices could play an important role in calorimetry and fast bolometry [77,
78, 79, 51]. We could envision similar applications for our device [76].
Furthermore, one could also exploit our system in thermoelectric nanodevices
similar to those designed in different platforms [80, 81, 82, 83] carrying the
advantage of higher operating temperatures. To conclude, we hope that the
peculiar features of the bipolar thermoelectric effect here reported will spur
further applications in the field of quantum technologies [32].

###### Acknowledgements.

The authors wish to thank for useful discussion Dr. A. Tomadin, Dr. T.
Novotny, Dr. K. Michaeli, and Prof. F. Strocchi. F.G. and A.B. acknowledge the
EU’s Horizon 2020 research and innovation program under Grant Agreement No.
800923 (SUPERTED) and No. 964398 (SUPERGATE) for partial financial support.
S.R. acknowledges the financial support of the Italian Ministry of University
and Research (PRIN project QUANTUM2D). A.B. acknowledges the Royal Society
through the International Exchanges between the UK and Italy (Grants No. IEC
R2 192166 and IEC R2 212041)

## References

  * Tritt [2004] T. M. Tritt,  _Thermal Conductivity: theory, properties, and applications_ (Kluwer Academic/Plenum Publishers, New York, 2004).
  * Tritt and Subramanian [2006] T. M. Tritt and M. A. Subramanian, Thermoelectric materials, phenomena, and applications: a bird’s eye view, MRS Bulletin 31, 188 (2006).
  * Macià-Barber [2015] E. Macià-Barber,  _Thermoelectric Materials: Advances and Applications_ (Pan Stanford Publishing, Taylor and Francis Group, 2015).
  * Grosso and Parravicini [2014] G. Grosso and G. P. Parravicini,  _Solid State Physics_ , 2nd ed. (Academic Press, Elsevier Ltd., 2014).
  * Benenti  _et al._ [2017] G. Benenti, G. Casati, K. Saito, and R. Whitney, Fundamental aspects of steady-state conversion of heat to work at the nanoscale, Phys. Rep. 694, 1 (2017).
  * Ginzburg [1991] V. Ginzburg, Thermoelectric effects in the superconducting state, Sov. Phys. Usp. 34, 101 (1991).
  * Smith  _et al._ [1980] A. D. Smith, M. Tinkham, and W. J. Skocpol, New thermoelectric effect in tunnel junctions, Phys. Rev. B 22, 4346 (1980).
  * Machon  _et al._ [2013] P. Machon, M. Eschrig, and W. Belzig, Nonlocal thermoelectric effects and nonlocal onsager relations in a three-terminal proximity-coupled superconductor-ferromagnet device, Phys. Rev. Lett. 110, 047002 (2013).
  * Ozaeta  _et al._ [2014] A. Ozaeta, P. Virtanen, F. S. Bergeret, and T. T. Heikkilä, Predicted very large thermoelectric effect in ferromagnet-superconductor junctions in the presence of a spin-splitting magnetic field, Phys. Rev. Lett. 112, 057001 (2014).
  * Sánchez  _et al._ [2018] R. Sánchez, P. Burset, and A. L. Yeyati, Cooling by cooper pair splitting, Phys. Rev. B 98, 241414 (2018).
  * Tabatabaei  _et al._ [2022] S. M. Tabatabaei, D. Sánchez, A. L. Yeyati, and R. Sánchez, Nonlocal quantum heat engines made of hybrid superconducting devices, Phys. Rev. B 106, 115419 (2022).
  * Pershoguba and Glazman [2019] S. S. Pershoguba and L. I. Glazman, Thermopower and thermal conductance of a superconducting quantum point contact, Phys. Rev. B 99, 134514 (2019).
  * Virtanen and Heikkilä [2004] P. Virtanen and T. Heikkilä, Thermopower in andreev interferometers, J. Low Temp. Phys. 136, 625 (2004).
  * Jacquod and Whitney [2010] P. Jacquod and R. Whitney, Coherent thermoelectric effects in mesoscopic andreev interferometers, Europhys. Lett. 91, 67009 (2010).
  * Kalenkov and Zaikin [2017] M. S. Kalenkov and A. D. Zaikin, Large thermoelectric effect in ballistic andreev interferometers, Phys. Rev. B 95, 024518 (2017).
  * Hussein  _et al._ [2019] R. Hussein, M. Governale, S. Kohler, W. Belzig, F. Giazotto, and A. Braggio, Nonlocal thermoelectricity in a cooper-pair splitter, Phys. Rev. B 99, 075429 (2019).
  * Kirsanov  _et al._ [2019] N. S. Kirsanov, Z. B. Tan, D. S. Golubev, P. J. Hakonen, and G. B. Lesovik, Heat switch and thermoelectric effects based on cooper-pair splitting and elastic cotunneling, Phys. Rev. B 99, 115127 (2019).
  * Ratnakar and Das [2021] A. Ratnakar and S. Das, Enhancement in tunneling density of states in a luttinger liquid: role of nonlocal interaction, Phys. Rev. B 104, 045402 (2021).
  * Eom  _et al._ [1998] J. Eom, C.-J. Chien, and V. Chandrasekhar, Phase dependent thermopower in andreev interferometers, Phys. Rev. Lett. 81, 437 (1998).
  * Jiang and Chandrasekhar [2005] Z. Jiang and V. Chandrasekhar, Quantitative measurements of the thermal resistance of andreev interferometers, Phys. Rev. B 72, 020502 (2005).
  * Tan  _et al._ [2021] Z. Tan, A. L. N. Kirsanov, A. Galda, V. Vinokur, M. Haque, A. Savin, D. Golubev, G. Lesovik, and P. Hakonen, Thermoelectric current in a graphene cooper pair splitter., Nat. Commun. 12, 138 (2021).
  * Aronov and Spivak [1975] A. G. Aronov and B. Z. Spivak, Photoeffect in a josephson junction, JETP Lett. 22, 101 (1975).
  * Gershenzon and Falei [1986] M. E. Gershenzon and M. I. Falei, Absolute negative resistance of a tunnel contact between superconductors with a nonequilibrium quasiparticle distribution function, JETP Lett. 44, 682 (1986).
  * Gershenzon and Falei [1988] M. E. Gershenzon and M. I. Falei, Absolute negative resistance in tunnel junctions of nonequilibrium superconductors, Sov. Phys. JETP 67, 389 (1988).
  * Gijsbertsen and Flokstra [1996] J. G. Gijsbertsen and J. Flokstra, Quasiparticle injection‐detection experiments in niobium, J. Appl. Phys. 80, 3923 (1996).
  * Marchegiani  _et al._ [2020a] G. Marchegiani, A. Braggio, and F. Giazotto, Nonlinear thermoelectricity with electron-hole symmetric systems, Phys. Rev. Lett. 124, 106801 (2020a).
  * Marchegiani  _et al._ [2020b] G. Marchegiani, A. Braggio, and F. Giazotto, Superconducting nonlinear thermoelectric heat engine, Phys. Rev. B 101, 214509 (2020b).
  * Marchegiani  _et al._ [2020c] G. Marchegiani, A. Braggio, and F. Giazotto, Phase-tunable thermoelectricity in a josephson junction, Phys. Rev. Research 2, 043091 (2020c).
  * Marchegiani  _et al._ [2020d] G. Marchegiani, A. Braggio, and F. Giazotto, Noise effects in the nonlinear thermoelectricity of a josephson junction, Appl. Phys. Lett. 117, 212601 (2020d).
  * Germanese  _et al._ [2022] G. Germanese, F. Paolucci, G. Marchegiani, A. Braggio, and F. Giazotto, Bipolar thermoelectric josephson engine, Nat. Nanotechnol. 17, 1084 (2022).
  * Germanese  _et al._ [2023] G. Germanese, F. Paolucci, G. Marchegiani, A. Braggio, and F. Giazotto, Phase control of bipolar thermoelectricity in josephson tunnel junctions, Phys. Rev. Appl. 19, 014074 (2023).
  * Giazotto  _et al._ [2021] F. Giazotto, F. Paolucci, A. Braggio, G. Marchegiani, and G. Germanese, Superconducting bipolar thermoelectric memory and method for writing a superconducting bipolar thermoelectric memory, Patent (21/12/2021), filing number: 102021000032042.
  * Germanese  _et al._ [2021] G. Germanese, F. Paolucci, G. Marchegiani, A. Braggio, and F. Giazotto, Spontaneous symmetry breaking induced thermospin effect in superconducting tunnel junctions, Phys. Rev. B 104, 184502 (2021).
  * Guarcello  _et al._ [2019] C. Guarcello, A. Braggio, P. Solinas, and F. Giazotto, Nonlinear critical-current thermal response of an asymmetric josephson tunnel junction, Phys. Rev. Applied 11, 024002 (2019).
  * Guarcello  _et al._ [2022] C. Guarcello, R. Citro, F. Giazotto, and A. Braggio, Temperature-biased double-loop josephson flux transducer, Phys. Rev. Applied 18, 014037 (2022).
  * Kittel [1986] C. Kittel,  _Introduction to Solid State Physics_ , 6th ed. (John Wiley & sons, Inc. (New York), 1986).
  * Tinkham [1996] M. Tinkham,  _Introduction to superconductivity_ , 2nd ed. (McGraw-Hill, 1996).
  * Bernazzani [2021] L. Bernazzani,  _Nonlinear Thermoelectricity in Nano-Devices_ , Master’s thesis, University of Pisa (2021).
  * Oostinga  _et al._ [2008] J. B. Oostinga, H. B. Heersche, X. Liu, A. F. Morpurgo, and L. M. K. Vandersypen, Gate-induced insulating state in bilayer graphene devices, Nature Mater. 7, 151 (2008).
  * Zhang  _et al._ [2009] Y. Zhang, T. Tang, C. Girit, Z. Hao, M. Martin, A. Zettl, M. Crommie, Y. R. Shen, and F. Wang, Direct observation of a widely tunable bandgap in bilayer graphene, Nature 459, 820 (2009).
  * Castro  _et al._ [2010] E. V. Castro, K. S. Novoselov, S. V. Morozov, N. M. R. Peres, J. M. B. L. dos Santos, J. Nilsson, F. Guinea, A. K. Geim, and A. H. Castro Neto, Electronic properties of a biased graphene bilayer, J. Phys. Condens. Matter 22, 175503 (2010).
  * McCann and Koshino [2013] E. McCann and M. Koshino, The electronic properties of bilayer graphene, Rep. Prog. Phys. 76, 056503 (2013).
  * [43] The lowest energy bands shown in the plot are obtained using Eq. 43 of Ref. [42] and relating the potential $U$ to the bandgap 2$E_{g}$ with Eq. (1).
  * Balandin [2011] A. Balandin, Thermal properties of graphene and nanostructured carbon materials, Nature Mater. 10, 569 (2011).
  * Dresselhaus  _et al._ [1999] M. S. Dresselhaus, G. Dresselhaus, X. Sun, Z. Zhang, S. B. Cronin, and T. Koga, Low-dimensional thermoelectric materials, Phys. Solid State 41, 679 (1999).
  * Novoselov  _et al._ [2005] K. Novoselov, A. Geim, S. Morozov, D. Jiang, M. Katsnelson, I. V. Grigorieva, S. V. Dubonos, and A. A. Firsov, Two-dimensional gas of massless dirac fermions in graphene, Nature 438, 197 (2005).
  * Sprinkle  _et al._ [2009] M. Sprinkle, D. Siegel, Y. Hu, J. Hicks, A. Tejeda, A. Taleb-Ibrahimi, P. Le Fèvre, F. Bertran, S. Vizzini, H. Enriquez, S. Chiang, P. Soukiassian, C. Berger, W. A. de Heer, A. Lanzara, and E. H. Conrad, First direct observation of a nearly ideal graphene band structure, Phys. Rev. Lett. 103, 226803 (2009).
  * McCann [2006] E. McCann, Asymmetry gap in the electronic band structure of bilayer graphene, Phys. Rev. B 74, 161403 (2006).
  * Sławińska  _et al._ [2010] J. Sławińska, I. Zasada, and Z. Klusek, Energy gap tuning in graphene on hexagonal boron nitride bilayer system, Phys. Rev. B 81, 155433 (2010).
  * Mak  _et al._ [2009] K. F. Mak, C. H. Lui, J. Shan, and T. F. Heinz, Observation of an electric-field-induced band gap in bilayer graphene by infrared spectroscopy, Phys. Rev. Lett. 102, 256405 (2009).
  * Vischi  _et al._ [2020] F. Vischi, M. Carrega, A. Braggio, F. Paolucci, F. Bianco, S. Roddaro, and F. Giazotto, Electron cooling with graphene-insulator-superconductor tunnel junctions for applications in fast bolometry, Phys. Rev. Applied 13, 054006 (2020).
  * Viljas and Heikkilä [2010] J. K. Viljas and T. T. Heikkilä, Electron-phonon heat transfer in monolayer and bilayer graphene, Phys. Rev. B 81, 245404 (2010).
  * Note [1] The top/bottom gates are referred to the BLG Ohmic contacts in order to induce a proper gating on the BLG which is not affected by the bias applied to the junction.
  * Maestre  _et al._ [2021] C. Maestre, B. Toury, P. Steyer, V. Garnier, and C. Journet, Hexagonal boron nitride: a review on selfstanding crystals synthesis towards 2d nanosheets, J. Phys. Mater. 4, 044018 (2021).
  * Tinkham [1972] M. Tinkham, Tunneling generation, relaxation, and tunneling detection of hole-electron imbalance in superconductors, Phys. Rev. B 6, 1747 (1972).
  * Ketterson and Song [1999] J. B. Ketterson and S. N. Song,  _Superconductivity_ (Cambridge University Press, 1999).
  * Giazotto  _et al._ [2006] F. Giazotto, T. T. Heikkilä, A. Luukanen, A. Savin, and J. P. Pekola, Opportunities for mesoscopics in thermometry and refrigeration: Physics and applications, Rev. Mod. Phys. 78, 217 (2006).
  * Dynes  _et al._ [1984] R. C. Dynes, J. P. Garno, G. B. Hertel, and T. P. Orlando, Tunneling study of superconductivity near the metal-insulator transition, Phys. Rev. Lett. 53, 2437 (1984).
  * Suárez Morell and Foa Torres [2012] E. Suárez Morell and L. E. F. Foa Torres, Radiation effects on the electronic properties of bilayer graphene, Phys. Rev. B 86, 125449 (2012).
  * Note [2] The hBN is commonly used to encapsulate graphene since its crystalline configuration matches well with that of graphene diminishing the surface disorder at the interfaces.
  * Note [3] The conductance $G_{T}=4N_{F}e^{2}|\mathcal{T}|^{2}t_{\perp}/(\hbar\sqrt{3}t^{2})$ would represent the differential conductance of the tunnel junction when the superconductor is in the normal state and the BLG is gapless.
  * Britnell  _et al._ [2012] L. Britnell, R. V. Gorbachev, R. Jalil, B. D. Belle, F. Schedin, M. I. Katsnelson, L. Eaves, S. V. Morozov, A. S. Mayorov, N. M. R. Peres, A. H. Castro Neto, J. Leist, A. K. Geim, L. A. Ponomarenko, and K. S. Novoselov, Electron tunneling through ultrathin boron nitride crystalline barriers, Nano Lett. 12, 1707 (2012).
  * Ramezani  _et al._ [2021] M. Ramezani, I. Sampaio, K. Watanabe, T. Taniguchi, C. Schönenberger, and A. Baumgartner, Superconducting contacts to a monolayer semiconductor, Nano Lett. 21, 5614 (2021).
  * Padilha  _et al._ [2012] J. E. Padilha, R. B. Pontes, and A. Fazzio, Bilayer graphene on h-bn substrate: investigating the breakdown voltage and tuning the bandgap by electric field, J. Phys. Condens. Matter 24, 075301 (2012).
  * Castro Neto  _et al._ [2009] A. H. Castro Neto, F. Guinea, N. M. R. Peres, K. S. Novoselov, and A. K. Geim, The electronic properties of graphene, Rev. Mod. Phys. 81, 109 (2009).
  * Dubi and Di Ventra [2011] Y. Dubi and M. Di Ventra, Colloquium: Heat flow and thermoelectricity in atomic and molecular junctions, Rev. Mod. Phys. 83, 131 (2011).
  * Note [4] In principle inverting the ratio between $E_{g}/\Delta_{0,R}$ also the opposite temperature gradient may be assumed. But this is less practical since the temperature of the superconducting lead should not exceed its critical temperature in any case.
  * Kheradsoud  _et al._ [2019] S. Kheradsoud, N. Dashti, M. Misiorny, P. P. Potts, J. Splettstoesser, and P. Samuelsson, Power, efficiency and fluctuations in a quantum point contact as steady-state thermoelectric heat engine, Entropy 21, 777 (2019).
  * Note [5] Even if the absolute value of those quantity depends also on the Dynes’ parameter $\Gamma_{D}$ the general behaviour would be quite independent on it.
  * Note [6] Obviously the heat may potentially be transferred also by other contributions such as phonon or even radiative contributions that are more difficult to be estimated. Our computation of efficiency represents a sort of optimal electronic efficiency of the device.
  * Alymov  _et al._ [2016] G. Alymov, V. Vyurkov, V. Ryzhii, and D. Svintsov, Abrupt current switching in graphene bilayer tunnel transistors enabled by van hove singularities, Sci. Rep. 6, 24654 (2016).
  * Icking  _et al._ [2022] E. Icking, L. Banszerus, F. Wörtche, F. Volmer, P. Schmidt, C. Steiner, S. Engels, J. Hesselmann, M. Goldsche, K. Watanabe, T. Taniguchi, C. Volk, B. Beschoten, and C. Stampfer, Transport spectroscopy of ultraclean tunable band gaps in bilayer graphene, Adv. Electron. Mater. 8, 2200510 (2022).
  * Note [7] The presence of localized charge states induced by impurities also determines a modification of the gating conditions anyway these contributions can be simply accounted by including appropriate offsets to the gate voltages, so we neglected them [72].
  * Note [8] The bandgap self-consistent expression of Eq. (6) is derived using Eqs. (65) and (74) of Ref. [42] in the limit with $\varepsilon_{r}\sim 1$, $\Lambda\sim 1$ and assuming no mobile charge density when the Fermi energy is in the gap. Correspondingly the Fermi-energy expression Eq. (7) coincides with the chemical potential of Eq. (S13) of the supplementary information of Ref. [72] where have been neglected terms involving the inverse of BLG capacitance $C_{BLG}^{-1}$ being $C_{BLG}\gg C_{b},C_{t}$.
  * Note [9] This simplified model may not be fully self-consistent when the system is gated well outside the gap, but corrections are expected to slightly renormalize the gap value not affecting physics discussed.
  * Paolucci  _et al._ [2023] F. Paolucci, G. Germanese, A. Braggio, and F. Giazotto, A highly-sensitive broadband superconducting thermoelectric single-photon detector, arXiv preprint arXiv:2302.02933 (2023).
  * Xia  _et al._ [2009] F. Xia, T. Mueller, Y. Lin, A. Valdes-Garcia, and P. Avouris, Ultrafast graphene photodetector, Nature Nanotech. 4, 839 (2009).
  * Mueller  _et al._ [2010] T. Mueller, F. Xia, and P. Avouris, Graphene photodetectors for high-speed optical communications, Nature Photon. 4, 297 (2010).
  * Lee  _et al._ [2020] G. H. Lee, D. K. Efetov, W. Jung, L. Ranzani, E. D. Walsh, T. A. Okhi, T. Taniguchi, K. Watanabe, P. K. D. Englund, and K. C. Fong, Graphene-based josephson junction microwave bolometer, Nature 586, 42 (2020).
  * Varpula  _et al._ [2017] A. Varpula, A. V. Timofeev, A. Shchepetov, K. Grigoras, J. Hassel, J. Ahopelto, M. Ylilammi, and M. Prunnila, Thermoelectric thermal detectors based on ultra-thin heavily doped single-crystal silicon membranes, Appl. Phys. Lett. 110, 262101 (2017).
  * Heikkilä  _et al._ [2018] T. T. Heikkilä, R. Ojajärvi, I. J. Maasilta, E. Strambini, F. Giazotto, and F. S. Bergeret, Thermoelectric radiation detector based on superconductor-ferromagnet systems, Phys. Rev. Applied 10, 034053 (2018).
  * Geng  _et al._ [2020] Z. Geng, A. P. Helenius, T. T. Heikkilä, and I. J. Maasilta, Superconductor-ferromagnet tunnel junction thermoelectric bolometer and calorimeter with a squid readout, J. Low Temp. Phys. 199, 585 (2020).
  * Varpula  _et al._ [2021] A. Varpula, K. Tappura, J. Tiira, K. Grigoras, O.-P. Kilpi, K. Sovanto, J. Ahopelto, and M. Prunnila, Nano-thermoelectric infrared bolometers, APL Photonics 6, 036111 (2021).


$\displaystyle=\dot{r}_{s}=2r\dot{r}=2\mathscr{g}_{r}.$

Finally, the time derivatives of $w_{h}$, $\xi_{1}$, $\xi_{2}$,
$\mathscr{s}_{L}$ and $\mathscr{c}_{L}$ are in Equation (38).

$\displaystyle w_{h}$ $\displaystyle=I\hat{z},$ $\displaystyle\dot{w}_{h}$
$\displaystyle=\dot{I}\hat{z}+I\dot{\hat{z}},$ (38) $\displaystyle\xi_{1}$
$\displaystyle=\frac{X}{a}+2p_{2},$ $\displaystyle\dot{\xi}_{1}$
$\displaystyle=\frac{\dot{X}}{a}+2\dot{p}_{2},$ $\displaystyle\xi_{2}$
$\displaystyle=\frac{Y}{a}+2p_{1},$ $\displaystyle\dot{\xi}_{2}$
$\displaystyle=\frac{\dot{Y}}{a}+2\dot{p}_{1},$ $\displaystyle\mathscr{s}_{L}$
$\displaystyle=\frac{Y}{r}=Y\mathscr{f}_{r},$
$\displaystyle\dot{\mathscr{s}}_{L}$
$\displaystyle=\dot{Y}\mathscr{f}_{r}+Y\dot{\mathscr{f}}_{r},$
$\displaystyle\mathscr{c}_{L}$ $\displaystyle=\frac{X}{r}=X\mathscr{f}_{r},$
$\displaystyle\dot{\mathscr{c}}_{L}$
$\displaystyle=\dot{X}\mathscr{f}_{r}+X\dot{\mathscr{f}}_{r}.$

In this section the computation of every $\mathscr{f}_{s}$ time derivative has
been explicitly noted. In the following this comments are omitted for the sake
of brevity, and it will be assumed all instances can be computed with the
corresponding derivatives of $s$ (to the needed order).

####  A.4.2 Third order derivatives

The third order time derivatives of the GEqOE can be computed with Equation
(39).

$\begin{split}\dddot{\nu}&=0,\\\
\dddot{q}_{1}&=-\ddot{I}\mathscr{s}_{L}-2\dot{I}\dot{\mathscr{s}}_{L}-I\ddot{\mathscr{s}}_{L},\\\
\dddot{q}_{2}&=-\ddot{I}\mathscr{c}_{L}-2\dot{I}\dot{\mathscr{c}}_{L}-I\ddot{\mathscr{c}}_{L},\\\
\dddot{p}_{1}&=\ddot{p}_{2}\left(\mathscr{d}-w_{h}\right)+2\dot{p}_{2}\left(\dot{\mathscr{d}}-\dot{w}_{h}\right)+p_{2}\left(\ddot{\mathscr{d}}-\ddot{w}_{h}\right)-\left(\ddot{\mathscr{f}}_{c}\xi_{1}+2\dot{\mathscr{f}}_{c}\dot{\xi}_{1}+\mathscr{f}_{c}\ddot{\xi}_{1}\right)\mathscr{U}-2\left(\dot{\mathscr{f}}_{c}\xi_{1}+\mathscr{f}_{c}\dot{\xi}_{1}\right)\dot{\mathscr{U}}-\mathscr{f}_{c}\xi_{1}\ddot{\mathscr{U}},\\\
\dddot{p}_{2}&=\ddot{p}_{1}\left(w_{h}-\mathscr{d}\right)+2\dot{p}_{1}\left(\dot{w}_{h}-\dot{\mathscr{d}}\right)+p_{1}\left(\ddot{w}_{h}-\ddot{\mathscr{d}}\right)+\left(\ddot{\mathscr{f}}_{c}\xi_{2}+2\dot{\mathscr{f}}_{c}\dot{\xi}_{2}+\mathscr{f}_{c}\ddot{\xi}_{2}\right)\mathscr{U}+2\left(\dot{\mathscr{f}}_{c}\xi_{2}+\mathscr{f}_{c}\dot{\xi}_{2}\right)\dot{\mathscr{U}}+\mathscr{f}_{c}\xi_{2}\ddot{\mathscr{U}},\\\
\dddot{\mathcal{L}}&=\ddot{\mathscr{d}}-\ddot{w}_{h}-\left(\ddot{\mathscr{f}}_{c}\Gamma+2\dot{\mathscr{f}}_{c}\dot{\Gamma}+\mathscr{f}_{c}\ddot{\Gamma}\right)\mathscr{U}-2\left(\dot{\mathscr{f}}_{c}\Gamma+\mathscr{f}_{c}\dot{\Gamma}\right)\dot{\mathscr{U}}-\mathscr{f}_{c}\Gamma\ddot{\mathscr{U}}.\end{split}$
(39)

All the elements that appear in Eq. (39) are computed with the following
formulation. Equation (40) computes the second time derivative of $\hat{z}$.

$\displaystyle\ddot{\hat{z}}$
$\displaystyle=\ddot{\mathcal{C}}\mathscr{f}_{\mathcal{D}}+2\dot{\mathcal{C}}\dot{\mathscr{f}}_{\mathcal{D}}+\mathcal{C}\ddot{\mathscr{f}}_{\mathcal{D}},$
$\displaystyle\ddot{\mathcal{C}}$
$\displaystyle=2\left(\ddot{Y}q_{2}+2\dot{Y}\dot{q}_{2}+Y\dot{q}_{2}-\ddot{X}q_{1}-2\dot{X}\dot{q}_{1}-X\ddot{q}_{1}\right),$
(40) $\displaystyle\ddot{\mathcal{D}}$
$\displaystyle=\ddot{r}q_{s}+2\dot{r}\dot{q}_{s}+r\ddot{q}_{s}\rightarrow\ddot{\mathscr{f}}_{\mathcal{D}},$
$\displaystyle\ddot{q}_{s}$
$\displaystyle=2\left(\dot{q}_{1}^{2}+q_{1}\ddot{q}_{1}+\dot{q}_{2}^{2}+q_{2}\ddot{q}_{2}\right).$

The second order time derivatives of $X$ and $Y$ are computed in Equation
(41).

$\displaystyle\dot{X}$
$\displaystyle=\dot{r}\cos{L}-\frac{h}{r}\sin{L}=\dot{r}\mathscr{c}_{L}-\mathscr{w}\mathscr{s}_{L},$
$\displaystyle\ddot{X}$
$\displaystyle=\ddot{r}\mathscr{c}_{L}+\dot{r}\dot{\mathscr{c}}_{L}-\dot{\mathscr{w}}\mathscr{s}_{L}-\mathscr{w}\dot{\mathscr{s}}_{L},$
(41) $\displaystyle\dot{Y}$
$\displaystyle=\dot{r}\sin{L}+\frac{h}{r}\cos{L}=\dot{r}\mathscr{s}_{L}+\mathscr{w}\mathscr{c}_{L},$
$\displaystyle\ddot{Y}$
$\displaystyle=\ddot{r}\mathscr{s}_{L}+\dot{r}\dot{\mathscr{s}}_{L}+\dot{\mathscr{w}}\mathscr{c}_{L}+\mathscr{w}\dot{\mathscr{c}}_{L},$
$\displaystyle\mathscr{w}$ $\displaystyle=\frac{h}{r}=h\mathscr{f}_{r},$
$\displaystyle\dot{\mathscr{w}}$
$\displaystyle=\dot{h}\mathscr{f}_{r}+h\dot{\mathscr{f}}_{r},$
$\displaystyle\dot{r}$
$\displaystyle=\frac{\mu}{c}\left(p_{2}\sin{L}-p_{1}\cos{L}\right)=\mu\mathscr{f}_{c}r_{pl},$
$\displaystyle\ddot{r}$
$\displaystyle=\mu\left(\dot{\mathscr{f}}_{c}r_{pl}+\mathscr{f}_{c}\dot{r}_{pl}\right),$
$\displaystyle r_{pl}$
$\displaystyle=p_{2}\mathscr{s}_{L}-p_{1}\mathscr{c}_{L},$
$\displaystyle\dot{r}_{pl}$
$\displaystyle=\dot{p}_{2}\mathscr{s}_{L}+p_{2}\dot{\mathscr{s}}_{L}-\dot{p}_{1}\mathscr{c}_{L}-p_{1}\dot{\mathscr{c}}_{L}.$

Equation (42) computes the second time derivative of $\mathscr{U}$ and
$U_{z}$.

$\displaystyle\ddot{\mathscr{U}}$
$\displaystyle=-A\left(\ddot{U}_{z}\mathscr{f}_{r^{3}}+2\dot{U}_{z}\dot{\mathscr{f}}_{r^{3}}+U_{z}\ddot{\mathscr{f}}_{r^{3}}\right),$
$\displaystyle\ddot{U}_{z}$
$\displaystyle=-6\left(\dot{\hat{z}}^{2}+\hat{z}\ddot{\hat{z}}\right),$
$\displaystyle\frac{d^{2}(r^{3})}{dt^{2}}=\ddot{r}_{c}=3\left(2r\dot{r}^{2}+r^{2}\ddot{r}\right).$
(42)

Equation (43) includes the second order derivatives of $\Gamma$ and $\beta$.

$\displaystyle\ddot{\Gamma}$
$\displaystyle=\ddot{\mathscr{f}}_{\alpha}+\ddot{\alpha}\left(1-\frac{r}{a}\right)-2\dot{\alpha}\frac{\dot{r}}{a}-\alpha\frac{\ddot{r}}{a},$
$\displaystyle\ddot{\alpha}$ $\displaystyle=\ddot{\mathscr{f}}_{1+\beta},$
(43) $\displaystyle\ddot{\beta}$
$\displaystyle=-\frac{1}{2}\left(\ddot{p}_{s}\mathscr{f}_{\beta}+\dot{p}_{s}\dot{\mathscr{f}}_{\beta}\right)\rightarrow\ddot{\mathscr{f}}_{1+\beta},\ddot{\mathscr{f}}_{\beta},$
$\displaystyle\ddot{p}_{s}$
$\displaystyle=2\left(\dot{p}_{1}^{2}+p_{1}\ddot{p}_{1}+\dot{p}_{2}^{2}+p_{2}\ddot{p}_{2}\right).$

Equation (44) computes the second time derivative of $I$.

$\displaystyle\ddot{I}$
$\displaystyle=3A\left[\left(\ddot{\hat{z}}\delta+2\dot{\hat{z}}\dot{\delta}+\hat{z}\ddot{\delta}\right)\mathscr{f}_{\mathscr{h}}+2\left(\dot{\hat{z}}\delta+\hat{z}\dot{\delta}\right)\dot{\mathscr{f}}_{\mathscr{h}}+\hat{z}\delta\ddot{\mathscr{f}}_{\mathscr{h}}\right],$
(44) $\displaystyle\ddot{\mathscr{h}}$
$\displaystyle=\ddot{h}r_{c}+2\dot{h}\dot{r}_{c}+h\ddot{r}_{c},$
$\displaystyle\ddot{h}$
$\displaystyle=\left(\dot{\mathscr{g}}_{c}-2\dot{\mathscr{g}}_{r}\mathscr{U}-4\mathscr{g}_{r}\dot{\mathscr{U}}-r^{2}\ddot{\mathscr{U}}\right)\mathscr{f}_{h}+\left(\mathscr{g}_{c}-2\mathscr{g}_{r}\mathscr{U}-r^{2}\dot{\mathscr{U}}\right)\dot{\mathscr{f}}_{h},$
$\displaystyle\ddot{c}$
$\displaystyle=\left(\frac{\mu^{2}}{\nu}\right)^{1/3}\ddot{\beta},$
$\displaystyle\ddot{\delta}$ $\displaystyle=-\ddot{q}_{s}.$

Equation (45) computes the second time derivative of $\mathscr{d}$.

$\displaystyle\ddot{\mathscr{d}}$
$\displaystyle=\left(\ddot{h}-\ddot{c}\right)\mathscr{f}_{r^{2}}+2\left(\dot{h}-\dot{c}\right)\dot{\mathscr{f}}_{r^{2}}+\left(h-c\right)\ddot{\mathscr{f}}_{r^{2}},$
(45) $\displaystyle\frac{d^{2}(r^{2})}{dt^{2}}$
$\displaystyle=\ddot{r}_{s}=2\dot{\mathscr{g}}_{r}.$

The second order time derivatives of $w_{h}$, $\xi_{1}$, $\xi_{2}$,
$\mathscr{s}_{L}$ and $\mathscr{c}_{L}$ are in Equation (46).

$\displaystyle\ddot{w}_{h}$
$\displaystyle=\ddot{I}\hat{z}+2\dot{I}\dot{\hat{z}}+I\ddot{\hat{z}},$ (46)
$\displaystyle\ddot{\xi}_{1}$
$\displaystyle=\frac{\ddot{X}}{a}+2\ddot{p}_{2},$
$\displaystyle\ddot{\xi}_{2}$
$\displaystyle=\frac{\ddot{Y}}{a}+2\ddot{p}_{1},$
$\displaystyle\ddot{\mathscr{s}}_{L}$
$\displaystyle=\ddot{Y}\mathscr{f}_{r}+2\dot{Y}\dot{\mathscr{f}}_{r}+Y\ddot{\mathscr{f}}_{r},$
$\displaystyle\ddot{\mathscr{c}}_{L}$
$\displaystyle=\ddot{X}\mathscr{f}_{r}+2\dot{X}\dot{\mathscr{f}}_{r}+X\ddot{\mathscr{f}}_{r}.$

####  A.4.3 Fourth order derivatives

The third order time derivatives of the GEqOE can be computed with Equation
(47).

$\displaystyle\overset{\scriptscriptstyle(4)}{\nu}=$ $\displaystyle 0,$ (47)
$\displaystyle\overset{\scriptscriptstyle(4)}{q}_{1}=$
$\displaystyle-\dddot{I}\mathscr{s}_{L}-3\ddot{I}\dot{\mathscr{s}}_{L}-3\dot{I}\ddot{\mathscr{s}}_{L}-I\dddot{\mathscr{s}}_{L},$
$\displaystyle\overset{\scriptscriptstyle(4)}{q}_{2}=$
$\displaystyle-\dddot{I}\mathscr{c}_{L}-3\ddot{I}\dot{\mathscr{c}}_{L}-3\dot{I}\ddot{\mathscr{c}}_{L}-I\dddot{\mathscr{c}}_{L},$
$\displaystyle\overset{\scriptscriptstyle(4)}{p}_{1}=$
$\displaystyle\dddot{p}_{2}\left(\mathscr{d}-w_{h}\right)+3\ddot{p}_{2}\left(\dot{\mathscr{d}}-\dot{w}_{h}\right)+3\dot{p}_{2}\left(\ddot{\mathscr{d}}-\ddot{w}_{h}\right)+p_{2}\left(\dddot{\mathscr{d}}-\dddot{w}_{h}\right)$
$\displaystyle-\left(\dddot{\mathscr{f}}_{c}\xi_{1}+3\ddot{\mathscr{f}}_{c}\dot{\xi}_{1}+3\dot{\mathscr{f}}_{c}\ddot{\xi}_{1}+\mathscr{f}_{c}\dddot{\xi}_{1}\right)\mathscr{U}-3\left(\ddot{\mathscr{f}}_{c}\xi_{1}+2\dot{\mathscr{f}}_{c}\dot{\xi}_{1}+\mathscr{f}_{c}\ddot{\xi}_{1}\right)\dot{\mathscr{U}}-3\left(\dot{\mathscr{f}}_{c}\xi_{1}+\mathscr{f}_{c}\dot{\xi}_{1}\right)\ddot{\mathscr{U}}-\mathscr{f}_{c}\xi_{1}\dddot{\mathscr{U}},$
$\displaystyle\overset{\scriptscriptstyle(4)}{p}_{2}=$
$\displaystyle\dddot{p}_{1}\left(w_{h}-\mathscr{d}\right)+3\ddot{p}_{1}\left(\dot{w}_{h}-\dot{\mathscr{d}}\right)+3\dot{p}_{1}\left(\ddot{w}_{h}-\ddot{\mathscr{d}}\right)+p_{1}\left(\dddot{w}_{h}-\dddot{\mathscr{d}}\right)$
$\displaystyle+\left(\dddot{\mathscr{f}}_{c}\xi_{2}+3\ddot{\mathscr{f}}_{c}\dot{\xi}_{2}+3\dot{\mathscr{f}}_{c}\ddot{\xi}_{2}+\mathscr{f}_{c}\dddot{\xi}_{2}\right)\mathscr{U}+3\left(\ddot{\mathscr{f}}_{c}\xi_{2}+2\dot{\mathscr{f}}_{c}\dot{\xi}_{2}+\mathscr{f}_{c}\ddot{\xi}_{2}\right)\dot{\mathscr{U}}+3\left(\dot{\mathscr{f}}_{c}\xi_{2}+\mathscr{f}_{c}\dot{\xi}_{2}\right)\ddot{\mathscr{U}}+\mathscr{f}_{c}\xi_{2}\dddot{\mathscr{U}},$
$\displaystyle\overset{\scriptscriptstyle(4)}{\mathcal{L}}=$
$\displaystyle\dddot{\mathscr{d}}-\dddot{w}_{h}-\left(\dddot{\mathscr{f}}_{c}\Gamma+3\ddot{\mathscr{f}}_{c}\dot{\Gamma}+3\dot{\mathscr{f}}_{c}\ddot{\Gamma}+\mathscr{f}_{c}\dddot{\Gamma}\right)\mathscr{U}-3\left(\ddot{\mathscr{f}}_{c}\Gamma+2\dot{\mathscr{f}}_{c}\dot{\Gamma}+\mathscr{f}_{c}\ddot{\Gamma}\right)\dot{\mathscr{U}}-3\left(\dot{\mathscr{f}}_{c}\Gamma+\mathscr{f}_{c}\dot{\Gamma}\right)\ddot{\mathscr{U}}-\mathscr{f}_{c}\Gamma\dddot{\mathscr{U}}.$

All the elements that appear in Eq. (47) are computed with the following
formulation. Equation (48) computes the second time derivative of $\hat{z}$.

$\displaystyle\dddot{\hat{z}}$
$\displaystyle=\dddot{\mathcal{C}}\mathscr{f}_{\mathcal{D}}+3\ddot{\mathcal{C}}\dot{\mathscr{f}}_{\mathcal{D}}+3\dot{\mathcal{C}}\ddot{\mathscr{f}}_{\mathcal{D}}+\mathcal{C}\dddot{\mathscr{f}}_{\mathcal{D}},$
$\displaystyle\dddot{\mathcal{C}}$
$\displaystyle=2\left(\dddot{Y}q_{2}+3\ddot{Y}\dot{q}_{2}+3\dot{Y}\ddot{q}_{2}+Y\dddot{q}_{2}-\dddot{X}q_{1}-3\ddot{X}\dot{q}_{1}-3\dot{X}\ddot{q}_{1}-X\dddot{q}_{1}\right),$
(48) $\displaystyle\dddot{\mathcal{D}}$
$\displaystyle=\dddot{r}q_{s}+3\ddot{r}\dot{q}_{s}+3\dot{r}\ddot{q}_{s}+r\dddot{q}_{s},$
$\displaystyle\dddot{q}_{s}$
$\displaystyle=2\left(3\dot{q}_{1}\ddot{q}_{1}+q_{1}\dddot{q}_{1}+3\dot{q}_{2}\ddot{q}_{2}+q_{2}\dddot{q}_{2}\right).$

The third order time derivatives of $X$ and $Y$ are computed in Equation (49).

$\displaystyle\dddot{X}$
$\displaystyle=\dddot{r}\mathscr{c}_{L}+2\ddot{r}\dot{\mathscr{c}}_{L}+\dot{r}\ddot{\mathscr{c}}_{L}-\ddot{\mathscr{w}}\mathscr{s}_{L}-2\dot{\mathscr{w}}\dot{\mathscr{s}}_{L}-\mathscr{w}\ddot{\mathscr{s}}_{L},$
(49) $\displaystyle\dddot{Y}$
$\displaystyle=\dddot{r}\mathscr{s}_{L}+2\ddot{r}\dot{\mathscr{s}}_{L}+\dot{r}\ddot{\mathscr{s}}_{L}+\ddot{\mathscr{w}}\mathscr{c}_{L}+2\dot{\mathscr{w}}\dot{\mathscr{c}}_{L}+\mathscr{w}\ddot{\mathscr{c}}_{L},$
$\displaystyle\ddot{\mathscr{w}}$
$\displaystyle=\ddot{h}\mathscr{f}_{r}+2\dot{h}\dot{\mathscr{f}}_{r}+h\ddot{\mathscr{f}}_{r},$
$\displaystyle\dddot{r}$
$\displaystyle=\mu\left(\ddot{\mathscr{f}}_{c}r_{pl}+2\dot{\mathscr{f}}_{c}\dot{r}_{pl}+\mathscr{f}_{c}\ddot{r}_{pl}\right),$
$\displaystyle\ddot{r}_{pl}$
$\displaystyle=\ddot{p}_{2}\mathscr{s}_{L}+2\dot{p}_{2}\dot{\mathscr{s}}_{L}+p_{2}\ddot{\mathscr{s}}_{L}-\ddot{p}_{1}\mathscr{c}_{L}-2\dot{p}_{1}\dot{\mathscr{c}}_{L}-p_{1}\ddot{\mathscr{c}}_{L}.$

Equation (50) computes the third order time derivative of $\mathscr{U}$ and
$U_{z}$.

$\displaystyle\dddot{\mathscr{U}}$
$\displaystyle=-A\left(\dddot{U}_{z}\mathscr{f}_{r^{3}}+3\ddot{U}_{z}\dot{\mathscr{f}}_{r^{3}}+3\dot{U}_{z}\ddot{\mathscr{f}}_{r^{3}}+U_{z}\dddot{\mathscr{f}}_{r^{3}}\right),$
(50) $\displaystyle\dddot{U}_{z}$
$\displaystyle=-6\left(3\dot{\hat{z}}\ddot{\hat{z}}+\hat{z}\dddot{\hat{z}}\right),$
$\displaystyle\frac{d^{3}(r^{3})}{dt^{3}}$
$\displaystyle=\dddot{r}_{c}=3\left(2\dot{r}^{3}+6r\dot{r}\ddot{r}+r^{2}\dddot{r}\right).$

Equation (51) includes the third order derivatives of $\Gamma$ and $\beta$.

$\displaystyle\dddot{\Gamma}$
$\displaystyle=\dddot{\mathscr{f}}_{\alpha}+\ddot{\alpha}\left(1-\frac{r}{a}\right)-3\ddot{\alpha}\frac{\dot{r}}{a}-3\dot{\alpha}\frac{\ddot{r}}{a}-\alpha\frac{\dddot{r}}{a},$
$\displaystyle\dddot{\alpha}$ $\displaystyle=\dddot{\mathscr{f}}_{1+\beta},$
(51) $\displaystyle\dddot{\beta}$
$\displaystyle=-\frac{1}{2}\left(\dddot{p}_{s}\mathscr{f}_{\beta}+2\ddot{p}_{s}\dot{\mathscr{f}}_{\beta}+\dot{p}_{s}\ddot{\mathscr{f}}_{\beta}\right),\dddot{\mathscr{f}}_{\beta},$
$\displaystyle\dddot{p}_{s}$
$\displaystyle=2\left(3\dot{p}_{1}\ddot{p}_{1}+p_{1}\dddot{p}_{1}+3\dot{p}_{2}\ddot{p}_{2}+p_{2}\dddot{p}_{2}\right).$

Equation (52) computes the third order time derivative of $I$.

$\displaystyle\dddot{I}$
$\displaystyle=3A\left[\left(\dddot{\hat{z}}\delta+3\ddot{\hat{z}}\dot{\delta}+3\dot{\hat{z}}\ddot{\delta}+\hat{z}\dddot{\delta}\right)\mathscr{f}_{\mathscr{h}}+3\left(\ddot{\hat{z}}\delta+2\dot{\hat{z}}\dot{\delta}+\hat{z}\ddot{\delta}\right)\dot{\mathscr{f}}_{\mathscr{h}}+3\left(\dot{\hat{z}}\delta+\hat{z}\dot{\delta}\right)\ddot{\mathscr{f}}_{\mathscr{h}}+\hat{z}\delta\dddot{\mathscr{f}}_{\mathscr{h}}\right],$
(52) $\displaystyle\dddot{\mathscr{h}}$
$\displaystyle=\dddot{h}r_{c}+3\ddot{h}\dot{r}_{c}+3\dot{h}\ddot{r}_{c}+h\dddot{r}_{c},$
$\displaystyle\dddot{h}$
$\displaystyle=\left(\ddot{\mathscr{g}}_{c}-2\ddot{\mathscr{g}}_{r}\mathscr{U}-6\dot{\mathscr{g}}_{r}\dot{\mathscr{U}}-6\mathscr{g}_{r}\ddot{\mathscr{U}}-r^{2}\dddot{\mathscr{U}}\right)\mathscr{f}_{h}+2\left(\dot{\mathscr{g}}_{c}-2\dot{\mathscr{g}}_{r}\mathscr{U}-4\mathscr{g}_{r}\dot{\mathscr{U}}-r^{2}\ddot{\mathscr{U}}\right)\dot{\mathscr{f}}_{h}+\left(\mathscr{g}_{c}-2\mathscr{g}_{r}\mathscr{U}-r^{2}\dot{\mathscr{U}}\right)\ddot{\mathscr{f}}_{h},$
$\displaystyle\dddot{c}$
$\displaystyle=\left(\frac{\mu^{2}}{\nu}\right)^{1/3}\dddot{\beta},$
$\displaystyle\dddot{\delta}$ $\displaystyle=-\dddot{q}_{s}.$

Equation (53) computes the third time derivative of $\mathscr{d}$.

$\displaystyle\dddot{\mathscr{d}}$
$\displaystyle=\left(\dddot{h}-\dddot{c}\right)\mathscr{f}_{r^{2}}+3\left(\ddot{h}-\ddot{c}\right)\dot{\mathscr{f}}_{r^{2}}+3\left(\dot{h}-\dot{c}\right)\ddot{\mathscr{f}}_{r^{2}}+\left(h-c\right)\dddot{\mathscr{f}}_{r^{2}},$
$\displaystyle\frac{d^{3}(r^{2})}{dt^{3}}=\dddot{r}_{s}=2\ddot{\mathscr{g}}_{r}.$
(53)

The third order time derivatives of $w_{h}$, $\xi_{1}$, $\xi_{2}$,
$\mathscr{s}_{L}$ and $\mathscr{c}_{L}$ are in Equation (54).

$\displaystyle\dddot{w}_{h}$
$\displaystyle=\dddot{I}\hat{z}+3\ddot{I}\dot{\hat{z}}+3\dot{I}\ddot{\hat{z}}+I\dddot{\hat{z}},$
(54) $\displaystyle\dddot{\xi}_{1}$
$\displaystyle=\frac{\dddot{X}}{a}+2\dddot{p}_{2},$
$\displaystyle\dddot{\xi}_{2}$
$\displaystyle=\frac{\dddot{Y}}{a}+2\dddot{p}_{1},$
$\displaystyle\dddot{\mathscr{s}}_{L}$
$\displaystyle=\dddot{Y}\mathscr{f}_{r}+3\ddot{Y}\dot{\mathscr{f}}_{r}+3\dot{Y}\ddot{\mathscr{f}}_{r}+Y\dddot{\mathscr{f}}_{r},$
$\displaystyle\dddot{\mathscr{c}}_{L}$
$\displaystyle=\dddot{X}\mathscr{f}_{r}+3\ddot{X}\dot{\mathscr{f}}_{r}+3\dot{X}\ddot{\mathscr{f}}_{r}+X\dddot{\mathscr{f}}_{r}.$

###  A.5 Derivatives of the $\mathscr{f}_{s}$ and $\mathscr{g}_{s}$ time
derivatives with respect to the element $x$

When $\overset{\scriptscriptstyle(k)}{\mathscr{f}}_{s}$ is derived with
respect to $x$, then all the derivatives up to
$\overset{\scriptscriptstyle(k)}{s}$ with respect to $x$ (with standard
notation $\overset{\scriptscriptstyle(k)}{s}_{x}$) are needed (or up to
$\overset{\scriptscriptstyle(k+1)}{s}$ in the case of
$\overset{\scriptscriptstyle(k)}{\mathscr{g}}_{s}$). The notation for the
cross derivative used here meets
$\partial\overset{\scriptscriptstyle(k)}{\mathscr{f}_{s}}/\partial
x=\overset{\scriptscriptstyle(k),x}{\mathscr{f}_{s}}$, chose on purpose for
the compact expressions developed in this appendix section.

Thus, the function that computes
$\overset{\scriptscriptstyle{(k),x}}{\mathscr{f}}_{s}$ needs not only the
vector $\left[s,\dot{s},\cdots,\overset{\scriptscriptstyle(k)}{s}\right]$, but
also
$\left[s_{x},\dot{s}_{x},\cdots,\overset{\scriptscriptstyle(k)}{s}_{x}\right]$
as input. Here are the expressions of the function that implements these
derivatives:

$\displaystyle\frac{\partial\mathscr{f}_{s}}{\partial
x}=\overset{\scriptscriptstyle{x}}{\mathscr{f}}_{s}=$
$\displaystyle-\frac{s_{x}}{s^{2}},$ (55)
$\displaystyle\frac{\partial\dot{\mathscr{f}}_{s}}{\partial
x}=\overset{\scriptscriptstyle{(1),x}}{\mathscr{f}}_{s}=$
$\displaystyle-\frac{\dot{s}_{x}}{s^{2}}+2\frac{\dot{s}s_{x}}{s^{3}},$
$\displaystyle\frac{\partial\ddot{\mathscr{f}}_{s}}{\partial
x}=\overset{\scriptscriptstyle{(2),x}}{\mathscr{f}}_{s}=$
$\displaystyle-\frac{\ddot{s}_{x}}{s^{2}}+2\frac{\ddot{s}s_{x}+2\dot{s}\dot{s}_{x}}{s^{3}}-6\frac{\dot{s}^{2}s_{x}}{s^{4}},$
$\displaystyle\frac{\partial\dddot{\mathscr{f}}_{s}}{\partial
x}=\overset{\scriptscriptstyle{(3),x}}{\mathscr{f}}_{s}=$
$\displaystyle-\frac{\dddot{s}_{x}}{s^{2}}+2\frac{\dddot{s}s_{x}+3\left(\ddot{s}_{x}\dot{s}+\ddot{s}\dot{s}_{x}\right)}{s^{3}}-18\frac{\ddot{s}\dot{s}s_{x}+\dot{s}^{2}\dot{s}_{x}}{s^{4}}+24\frac{\dot{s}^{3}s_{x}}{s^{5}},$
$\displaystyle\frac{\partial\overset{\scriptscriptstyle{(4)}}{\mathscr{f}}_{s}}{\partial
x}=\overset{\scriptscriptstyle{(4),x}}{\mathscr{f}}_{s}=$
$\displaystyle-\frac{\overset{\scriptscriptstyle(4)}{s}_{x}}{s^{2}}+2\frac{\overset{\scriptscriptstyle(4)}{s}s_{x}+6\ddot{s}\ddot{s}_{x}+4\left(\dddot{s}_{x}\dot{s}+\overset{\scriptscriptstyle(3)}{s}\dot{s}_{x}\right)}{s^{3}}-6\frac{4\overset{\scriptscriptstyle(3)}{s}\dot{s}s_{x}+6\left(2\dot{s}\dot{s}_{x}\ddot{s}+\dot{s}^{2}\ddot{s}_{x}\right)+3\ddot{s}^{2}s_{x}}{s^{4}}$
$\displaystyle+12\frac{12\dot{s}^{2}\ddot{s}s_{x}-8\dot{s}^{3}\dot{s}_{x}}{s^{5}}+120\frac{\dot{s}^{4}s_{x}}{s^{6}}.$

And similarly for $\overset{\scriptscriptstyle{(k),x}}{\mathscr{g}}_{s}$ the
derivatives are:

$\displaystyle\frac{\partial\mathscr{g}_{s}}{\partial
x}=\overset{\scriptscriptstyle{x}}{\mathscr{g}}_{s}=$ $\displaystyle
s_{x}\dot{s}+s\dot{s}_{x},$ (56)
$\displaystyle\frac{\partial\dot{\mathscr{g}}_{s}}{\partial
x}=\overset{\scriptscriptstyle{(1),x}}{\mathscr{g}}_{s}=$ $\displaystyle
2\dot{s}\dot{s}_{x}+s_{x}\ddot{s}+s\ddot{s}_{x},$
$\displaystyle\frac{\partial\ddot{\mathscr{g}}_{s}}{\partial
x}=\overset{\scriptscriptstyle{(2),x}}{\mathscr{g}}_{s}=$ $\displaystyle
3\dot{s}_{x}\ddot{s}+3\dot{s}\ddot{s}_{x}+s_{x}\dddot{s}+s\dddot{s}_{x},$
$\displaystyle\frac{\partial\dddot{\mathscr{g}}_{s}}{\partial
x}=\overset{\scriptscriptstyle{(3),x}}{\mathscr{g}}_{s}=$ $\displaystyle
6\ddot{s}\ddot{s}_{x}+4\dot{s}_{x}\overset{\scriptscriptstyle(3)}{s}+4\dot{s}\overset{\scriptscriptstyle(3)}{s}_{x}+s_{x}\overset{\scriptscriptstyle(4)}{s}+s\overset{\scriptscriptstyle(4)}{s}_{x},$
$\displaystyle\frac{\partial\overset{\scriptscriptstyle{(4)}}{\mathscr{g}}_{s}}{\partial
x}=\overset{\scriptscriptstyle{(4),x}}{\mathscr{g}}_{s}=$ $\displaystyle
10\ddot{s}_{x}\overset{\scriptscriptstyle(3)}{s}+10\ddot{s}\overset{\scriptscriptstyle(3)}{s}_{x}+5\dot{s}_{x}\overset{\scriptscriptstyle(4)}{s}+5\dot{s}\overset{\scriptscriptstyle(4)}{s}_{x}+s_{x}\overset{\scriptscriptstyle(5)}{s}+s\overset{\scriptscriptstyle(5)}{s}_{x}.$

###  A.6 Generalized equinoctial orbital elements auxiliary functions

Algorithms 1 and 2 are basic conversion functions needed to work with the
generalized equinoctial orbital elements. The analytical jacobians of these
conversions are given explicitly in Baù et al. [2021].

Algorithm 1
$\>\bm{\chi}=\left[\nu\>p_{1}\>p_{2}\>q_{1}\>q_{2}\>\mathcal{L}\right]^{\intercal}=\text{RV2GEQOE}_{fun}(\bm{r}\rfloor_{ECI},\bm{v}\rfloor_{ECI},\mu,J_{2},R_{\oplus})$

1:$r=\lVert\bm{r}\lVert$, $\mathscr{v}=\lVert\bm{v}\lVert$

2:$\bm{h}=\bm{r}\times\bm{v},\>h=\lVert\bm{h}\lVert$,
$\dot{r}=\frac{\bm{r}\cdot\bm{v}}{r}$

3:$\hat{z}=\bm{r}(3)/r$

4:$A=\frac{\mu\text{J}_{2}R_{\oplus}^{2}}{2}$,
$\mathscr{U}=-\frac{A}{r^{3}}(1-3\hat{z}^{2})$,
$\mathscr{U}_{\text{eff}}=\frac{h^{2}}{2r^{2}}+\mathscr{U}$

5:$\bm{e}_{r}=\bm{r}/r$, $\bm{e}_{h}=\bm{h}/h$,
$\bm{e}_{f}=\bm{e}_{h}\times\bm{e}_{r}$

6:$\bm{e}_{x}=\left[1,0,0\right]^{\intercal}$,
$\bm{e}_{y}=\left[0,1,0\right]^{\intercal}$,
$\bm{e}_{z}=\left[0,0,1\right]^{\intercal}$

7:Total energy:
$\varepsilon=\varepsilon_{K}+\mathscr{U}=\frac{\mathscr{v}^{2}}{2}-\frac{\mu}{r}+\mathscr{U}$

8:Generalized mean motion: $\nu=\frac{1}{\mu}\left(-2\varepsilon\right)^{3/2}$

9:$q_{1}=\frac{e_{h}\dot{e}_{x}}{1+e_{h}\dot{e}_{z}}$,
$q_{2}=\frac{-e_{h}\dot{e}_{x}}{1+e_{h}\dot{e}_{z}}$

10:$\bm{e}_{X}=\frac{1}{1+q_{1}^{2}+q_{2}^{2}}\left[1-q_{1}^{2}+q_{2}^{2},2q_{1}q_{2},-2q_{1}\right]^{\intercal}$

11:$\bm{e}_{Y}=\frac{1}{1+q_{1}^{2}+q_{2}^{2}}\left[2q_{1}q_{2},1+q_{1}^{2}-q_{2}^{2},2q_{2}\right]^{\intercal}$

12:Generalized angular momentum: $c=\sqrt{2r^{2}\mathscr{U}_{\text{eff}}}$

13:Generalized velocity vector:
$\bm{\upsilon}=\dot{r}\bm{e}_{r}+\frac{c}{r}\bm{e}_{f}$

14:Generalized eccentricity vector:
$\bm{g}=\frac{1}{\mu}\bm{\upsilon}\left(\bm{r}\times\bm{\upsilon}\right)-\bm{e}_{r}$

15:$p_{1}=\bm{g}\cdot\bm{e}_{Y}$, $p_{2}=\bm{g}\cdot\bm{e}_{X}$

16:$X=\bm{r}\cdot\bm{e}_{X}$, $Y=\bm{r}\cdot\bm{e}_{Y}$

17:$\beta=\sqrt{1-p1^{2}-p_{2}^{2}}$, $\alpha=\frac{1}{1+\beta}$

18:Generalized semi-major axis:
$a=-\frac{\mu}{2\varepsilon}=\left(\frac{\mu}{\nu^{2}}\right)^{1/3}$

19:$\cos{\mathcal{K}}=p_{2}+\frac{1}{a\beta}\left[\left(1-\alpha
p_{2}^{2}\right)X-\alpha p_{1}p_{2}Y\right]$

20:$\sin{\mathcal{K}}=p_{1}+\frac{1}{a\beta}\left[\left(1-\alpha
p_{1}^{2}\right)Y-\alpha p_{1}p_{2}X\right]$

21:$\mathcal{L}=\text{atan2}(\sin{\mathcal{K}},\cos{\mathcal{K}})+\frac{1}{a\beta}\left(Xp_{1}-Yp_{2}\right)$

Algorithm 2
$\>\bm{X}\rfloor_{ECI}=\left[\bm{r}^{\intercal},\bm{v}^{\intercal}\right]^{\intercal}=\text{GEQOE2RV}_{fun}(\nu,p_{1},p_{2},q_{1},q_{2},\mathcal{L},\mu,\text{J}_{2},R_{\oplus})$

1:Numerically solve Kepler's equation:
$\mathcal{K}\leftarrow\mathcal{L}=\mathcal{K}+p_{1}\cos{\mathcal{K}}-p_{2}\sin{\mathcal{K}}$

2:$\bm{e}_{X}=\frac{1}{1+q_{1}^{2}+q_{2}^{2}}\left[1-q_{1}^{2}+q_{2}^{2},2q_{1}q_{2},-2q_{1}\right]^{\intercal}$

3:$\bm{e}_{Y}=\frac{1}{1+q_{1}^{2}+q_{2}^{2}}\left[2q_{1}q_{2},1+q_{1}^{2}-q_{2}^{2},2q_{2}\right]^{\intercal}$

4:$\beta=\sqrt{1-p1^{2}-p_{2}^{2}}$, $\alpha=\frac{1}{1+\beta}$

5:$a=\left(\frac{\mu}{\nu^{2}}\right)^{1/3}$

6:$X=a\left[\alpha p_{1}p_{2}\sin{\mathcal{K}}+(1-\alpha
p_{1}^{2})\cos{\mathcal{K}}-p_{2}\right]$

7:$Y=a\left[\alpha p_{1}p_{2}\cos{\mathcal{K}}+(1-\alpha
p_{2}^{2})\sin{\mathcal{K}}-p_{1}\right]$

8:$\bm{r}=X\bm{e}_{X}+Y\bm{e}_{Y}$

9:$r=a\left(1-p_{1}\sin{\mathcal{K}}-p_{2}\cos{\mathcal{K}}\right)$

10:$\hat{z}=\bm{r}(3)/r$

11:$A=\frac{\mu\text{J}_{2}R_{\oplus}^{2}}{2}$,
$\mathscr{U}=-\frac{A}{r^{3}}(1-3\hat{z}^{2})$

12:$\dot{r}=\frac{\sqrt{\mu
a}}{r}\left(p_{2}\sin{\mathcal{K}}-p_{1}\cos{\mathcal{K}}\right)$

13:$\cos{L}=\frac{X}{r}$, $\sin{L}=\frac{Y}{r}$

14:$c=\left(\frac{\mu^{2}}{\nu}\right)^{1/3}\sqrt{1-q_{1}^{2}-q_{2}^{2}}$,
$h=\sqrt{c^{2}-2r^{2}\mathscr{U}}$

15:$\dot{X}=\dot{r}\cos{L}-\frac{h}{r}\sin{L}$

16:$\dot{Y}=\dot{r}\sin{L}+\frac{h}{r}\cos{L}$

17:$\bm{v}=\dot{X}\bm{e}_{X}+\dot{Y}\bm{e}_{Y}$