text:
\section{Introduction}
The Metal Oxide Semiconductor Field Effect Transistor (MOSFET) is the backbone of the modern integrated circuits. In the case  the active area is made of traditional materials like, for example, silicon or gallium arsenide, a lot of analysis and simulations have been performed in order to optimize the design. Lately a great attention has been devoted to graphene  on account of its peculiar features, and in particular, from the point of view of nano-electronics, for the high
electrical conductivity. It is highly tempting to try to replace the traditional semiconductors with graphene in the active area of electron devices like the MOSFETs \cite{Sch}   even if many aspect  about the actual performance  in real applications
remain unclear.

Scaling
theory predicts that a FET with a thin barrier and a thin gate-controlled
region will be robust
against short-channel effects down to very short gate lengths. The possibility of having
channels that are just one atomic layer thick is perhaps the most
attractive feature of graphene for its use in transistors.
Main drawbacks of a large-area single monolayer graphene are the  zero gap and, for graphene on substrate, the degradation of the mobility. Therefore accurate simulation are warranted for the set up of a viable
 graphene field effect transistor.

The standard mathematical model is  given by the drift-diffusion-Poisson system.
Usually the GFETs are investigated by adopting reduced one dimensional models of the Poisson equation with some averaging procedure \cite{Jimenez,Upadhyay,Meric}. Here a full two-dimensional simulation is presented.

A crucial point is the determination of the mobilities entering the drift-diffusion equations. A rather popular model is that proposed in \cite{Dor}. Here a different approach is adopted. Thanks to the discontinuous Galerkin (DG) scheme developed in \cite{RoMajCo,CoMajRo,ART:CoMajNaRo_AAPP,ART:MajNaRo_CICP}, we have performed, for graphene on a substrate, an extensive numerical simulation based on the semiclassical Boltzmann equations, including electron-phonon and electron-impuritiy scatterings along with the scattering with remote phonons of the substrate. Both  intra and inter-band scatterings  have been taken into account. For other simulation approaches the interested reader is referred to \cite{LMS,MuWa}. Quantum effects have been also introduced in \cite{Bar,Morandi11,LuRo1,MaRo2,MaRo_entropy}. An alternative approach could be resorting to hydrodynamical models
(see \cite{CaRo,LuRo2}).

From the numerical solutions of the semiclassical Boltzmann equation a model for the mobility functions has been deduced, similarly to  what already done in \cite{MajMaRo} and in \cite{ART:NaRo_CAIM} in the case of suspended monolayer graphene.

In the present paper we propose a slight different geometry of GFET which leads to a clear and sizable off region, avoiding the problem of triggering the minority charge carriers and producing characteristic curves which seem very promising for the use of large-area graphene  in electron devices.

The plan of the paper is as follows. In Sec. 2 the structure of the proposed device is presented along with the  mathematical model. In Sec. 3 the mobilities are sketched  and in Sec. 4  the numerical simulations of the proposed GFET are shown.

\section{Device structure and mathematical model}
We propose the device with the geometry  depicted in Fig.~\ref{FIG:GMOSFET}. The active zone is made of a single layer of graphene which is between two strips of insulator, both of  them being  SiO$_2$. The source and drain contacts are directly attached to the graphene. The two gate contacts (up and down) are attached to the oxide.  The choice of the type of oxide is not crucial at this stage. We have considered SiO$_2$ as an example. Instead it is crucial to put the source and drain contacts along all the lateral edges for a better control of the electrostatic potential as will be clear from the simulations. In the direction orthogonal to the section the device is considered as infinitely long.

We solve a 2D Poisson equation for the electrostatic potential, assuming that the charge is concentrated on the volume occupied by the atoms composing the graphene layer. Accordingly, the surface charge density of graphene is supposed to be spanned on its thickness which experimental measurements refer between 0.4 nm and 1.7 nm \cite{Shearer}. In order to simulate the current flowing in the channel we adopt the 1D bipolar drift-diffusion model, coupled to the Poisson equation for the electrostatic potential in the whole section. A special attention is required by the initial carrier density profiles that have to be determined compatibly with the electric potential, leading to a nonlinear Poisson equation as  discussed in the following.
\begin{figure}[!t]
\centerline{\includegraphics[width=0.8\columnwidth]{Fig/DG_GFET.pdf}}
\caption{Schematic representation of the GFET investigated in the paper.}\label{FIG:GMOSFET}
\end{figure}
The bipolar drift-diffusion model is adopted in $[x_1,x_4]\times\lbrace y_{gr}\rbrace$ and it reads
\begin{eqnarray*}
\frac{\partial n}{\partial t} - \frac{\partial}{\partial x}\left( \mu_n U_T \frac{\partial n}{\partial x}- n\mu_n\frac{\partial \phi}{\partial x} \right)=0,\\
\frac{\partial p}{\partial t} + \frac{\partial}{\partial x}\left( -\mu_p U_T\frac{\partial p}{\partial x}- p\mu_p\frac{\partial \phi}{\partial x} \right)=0,
\end{eqnarray*}
where $n(t,x)$, $p(t,x)$ are the electron and hole densities in  graphene respectively, at time $t$ and position $x$, $U_T=k_B T/e$ is the thermal voltage, being $e$ the positive elementary charge, $k_B$ is the Boltzmann constant, $T$ is the lattice temperature (kept constant). The functions $\mu_n$ and $\mu_p$ are the mobilities for electrons and holes respectively and $\phi(x,y)$ is the electric potential, here evaluated on $y=y_{gr}$, $y_{gr}$ being  the average $y$-coordinate of the graphene sheet (see Fig. \ref{FIG:GMOSFET}). The generation and recombination terms are set equal to zero \cite{ART:NaRo_CAIM}. Indeed, this relation is strictly valid at steady state but here will be assumed during the transient as well. We expect that the stationary solutions is not affected by such an approximation.  A typical behavior of the total generation and recombination term versus time, obtained with the DG method \cite{RoMajCo,CoMajRo,ART:CoMajNaRo_AAPP,ART:MajNaRo_CICP}, can be found in \cite{NaRo_DGFET}.

The system is solved in the interval $[x_1,x_4]$ augmented with Dirichlet boundary conditions for electron and hole densities as
will be explained below.

 The electric potential solves the 2D Poisson equation
\begin{equation}\label{EQ:Poiss}
\nabla\cdot(\epsilon\nabla\phi)=h(x,y),
\end{equation}
where
$h(x,y) = e(n(x)-p(x))/t_{gr}$ if $\, (x, y) \in [x_1, x_4] \times [y_2, y_3] $ and
$h(x,y) = 0$ otherwise,
and $\epsilon$ is given by
$\epsilon_{gr} $  if $y \in [y_2, y_3]$ and
$\epsilon_{ox}$ otherwise.
Here $\epsilon_{gr} = 3.3 \, \epsilon_0$ and $\epsilon_{ox} 3.6 \, \epsilon_0$ are the dielectric constants of the graphene and oxide (SiO$_2$) respectively,  $\epsilon_0$ being the dielectric constant in the vacuum;
 $t_{gr}$ is the width of the graphene layer
which is assumed to be 1 nm (values between 0.4 and 1.7 nm are reported in \cite{Shearer}).  The charge in the graphene layer is considered as distributed in the volume enclosed by the parallelepiped of base the area of the graphene  and height $t_{gr}$. Recall that $n$ and $p$ are areal densities. Dirichlet conditions are imposed on the  gate contacts and homogeneous Neumann conditions on the external oxide edges.
A major issue is to model the source and drain regions where metal and graphene touch. We assume that source and drain are thermal bath reservoirs of  charges which obey a Fermi-Dirac distribution.  The injection of charges is determined by a work function $W_F$.  Indeed it depends on the specific material the contacts are made of.
We set $W_F = 0.25 V$ which is appropriate for Cu within the experimentally reported range of 0.20 eV  \cite{FraVe} and 0.30 eV  \cite{WaNi}.

As summary the following boundary conditions for the electric potential  are imposed
\begin{eqnarray*}
\begin{array}{ll}
\phi = V_S = W_F & \quad \mbox{at} \quad  x= x_1, \, y \in [y_1, y_4]\\
\phi = V_D = W_F  + V_b & \quad \mbox{at} \quad  x= x_4, \, y \in [y_1, y_4]\\
\phi = V_{Gu} = W_F + V_{G} &\quad \mbox{at} \quad y=y_4, \, x \in [x_2, x_3]\\
\phi = V_{Gd} = W_F + V_{G} &\quad \mbox{at} \quad y=y_1, \, x \in [x_2, x_3]\\
\nabla_{\nu}\phi = 0 & \quad \mbox{at the remaining part of the boundary}.
\end{array}
\end{eqnarray*}
Here $V_b$ is the bias voltage, $V_{S}$ is the source potential, $V_{D}$ is the source potential, $V_{Gu}$ is the upper gate potential, $V_{Gd}$ is the down gate potential. $V_{G}$ is the gate bias potential which is considered equal at both the gate contacts.  $\nabla_{\nu}$ denotes the normal derivative.

 In graphene a sort of doping is induced by the electrostatic potential  \cite{Landauer} which leads to a shift of the Fermi energy. Assuming thermal equilibrium,  the initial carrier densities $n_0(x)$ and $ p_0(x)$ of electrons and holes respectively are related to the electric potential by
\begin{equation*}\label{EQ:pot_dens_rel}
\begin{alignedat}{3}
& n_0(x) && = \frac{2}{(2\pi)^2}\int \! f_{FD}(\mathbf{k}; e\phi(x,y_{gr})) \, d\mathbf{k}, && \quad x \in [x_1, x_4], \\
& p_0(x) && = \frac{2}{(2\pi)^2}\int \! f_{FD}(\mathbf{k}; -e\phi(x,y_{gr})) \, d\mathbf{k}, && \quad x \in [x_1, x_4],
\end{alignedat}
\end{equation*}
being $f_{FD}$ the Fermi-Dirac distribution
\begin{equation*}
f_{FD}(\mathbf{k};\varepsilon_F) = \left[1+\exp\left( \frac{\varepsilon(\mathbf{k})-\varepsilon_F}{k_B T} \right) \right]^{-1},
\end{equation*}
where $\varepsilon_F$ is the Fermi level (in pristine graphene $\varepsilon_F = 0$), $\varepsilon(\mathbf{k}) = \hbar v_F\vert \mathbf{k}\vert $ is the graphene dispersion relation (strictly valid around the Dirac points),
which is the same for electrons and holes (see \cite{BOOK:Jacoboni,Kittel,CaNe}), $\hbar$ is the reduced Planck constant and $v_F$ is the Fermi velocity. The crystal momentum of electrons and holes is assumed to vary over  $\mathbb{R}^2$.
The boundary conditions at the contacts are given by
\begin{eqnarray}
n(t, x_i) = n_0 (x_i), \quad p(t, x_i) = p_0 (x_i), \quad i = 1, 4. \label{Poisson_BC}
\end{eqnarray}

Altogether, in order to get the initial density profile as function of the electric potential, we must solve the following nonlinear Poisson equation
\begin{equation}
\nabla\cdot(\epsilon\nabla\phi)=g(\phi), \label{IniPoiss1}
\end{equation}
augmented with the boundary conditions (\ref{Poisson_BC}),
where
\begin{equation}
g(\phi(x,y)) =\left\lbrace
\begin{aligned}
&e \left[ n_0(\phi(x,y))-p_0(\phi(x,y))\right]/t_{gr} \qquad \mbox{if} \quad  (x, y) \in  [x_1, x_4]\times [y_2, y_3]\\
&0 \qquad\qquad\qquad\qquad\qquad\qquad\qquad\quad \mbox{otherwise}.
\end{aligned}
\right.
\end{equation}

\section{Mobility model}
\label{section:mobility}
We adopt the mobilities deduced in \cite{NaRo_DGFET} from a direct numerical simulation of the transport equations by using a DG scheme (the interested reader is referred to \cite{CoMajRo,ART:MajNaRo_CICP} for the details). The  behaviour is the same for holes on account of the symmetry between the hole and electron distributions.

The low field mobility has the expression
\begin{eqnarray*}
\mu_0(n) &=& \tilde{\mu}_1 - \tilde{\mu}_0 \frac{\exp\left( -\frac{(\log(n/n_{ref})-m)^2}{2 \sigma^2} \right)}{\sqrt{2\pi} \sigma n/n_{ref}}
 \left( a\left(\frac{n}{n_{ref}}\right)^2+b\frac{n}{n_{ref}}+c\right), \label{lowmu}
\end{eqnarray*}
where the fitting parameters have been estimated by the least squares method as follows
 $\tilde{\mu}_0$ = 0.2978 $\mu$m$^2$/V ps , $\tilde{\mu}_1 = $ 4.223  $\mu$m$^2$/V ps, $n_{ref}$ =  376.9 $\mu$m$^{-2}$, $m = - $0.2838, $\sigma$ = 2.216, $a$ =  4.820, $b$ = 68.34 and $c$ = 2.372.

The complete mobility model is given by (see also \cite{ART:NaRo_CAIM},\cite{MajMaRo})
\begin{equation}
\mu(E,n) = \frac{\mu_0(n)+\tilde{\mu}\left( \frac{E}{E_{ref}} \right)^{\beta_1}}{1+\left( \frac{E}{E_{ref}} \right)^{\beta_2}+\gamma\left( \frac{E}{E_{ref}} \right)^{\beta_3}}, \label{highmu}
\end{equation}
where $E_{ref}$, $\tilde{\mu}$, $\beta_1$, $\beta_2$, $\beta_3$ and $\gamma$ are fitting parameters.

We have calculated the coefficients $E_{ref}$, $\beta_1$, $\beta_2$, $\beta_3$, $\gamma$ and $\tilde{\mu}$ by means of least square method for several value of the electron density $n_i$, obtaining the data reported in Tab.~\ref{TAB:High_field_mob_param}. The values of the density correspond to the Fermi energies 0.1, 0.2, 0.3, 0.4, 0.5 eV.
\begin{table}[htp]
\centering
\begin{tabular}{c|c|c|c|c|c|c}
$n_i$ & $E_{ref}$ & $\beta_1$ & $\beta_2$ & $\beta_3$ & $\gamma$ & $\tilde{\mu}$\\
\hline
4471.0 & 0.05265 & 1.034 & 2.135 & 1.059 & 14.53 & 12.78\\
15500.0 & 0.02126 & 0.4615 & 1.584 & 0.4276 & 21.77 & 33.3\\
33877.0 & 0.1096 & 1.344 & 2.52 & 1.457 & 4.595 & 6.579\\
59588.0 & 0.1776 & 2.304 & 3.099 & 1.335 & 1.395 & 1.251\\
92644.0  & 0.05047 & 1.988 & 2.661 & 1.109 & 0.8816 & 1.915
\end{tabular}
\caption{High field mobility parameters.}
\label{TAB:High_field_mob_param}
\end{table}
In each interval $[n_i, n_{i+1}]$  a third degree polynomial interpolation has been adopted for the parameters $E_{ref}$, $\beta_1$, $\beta_2$, $\beta_3$, $\gamma$ and $\tilde{\mu}$ \cite{NaRo_DGFET}.

\section{Numerical simulations}
The GFET of Fig.~\ref{FIG:GMOSFET} has been intensively simulated. We have set
the length 100 nm, the  width of the lower and upper oxide (SiO$_2$) 10 nm. The source and drain contacts are  long 50 nm. The lateral contacts are long 21 nm.
The two gate potentials are set equal in the simulated cases. We considered a mesh of 40  grid points along the $x$-direction and 23 grid points along the y-direction. In the graphene layer a single row of 40 nodes has been employed.
In order to get the solution the following strategy has been adopted:
first the Poisson equation is solved by keeping the charge in the graphene layer equal to $n_0(x)$ and $p_0(x)$ for electron and holes, respectively; then the nonlinear Poisson problem (\ref{IniPoiss1}) is solved with an iterative scheme, by taking as initial guess the solution of  the previous step;
once the initial data for the electron and hole density have been determined, the full transient drift-diffusion-Poisson system is solved until the steady state is reached. One gets the stationary regime in about three picoseconds.
\begin{figure}[ht]
\centerline{\includegraphics[width=0.49\textwidth]{Fig/tot_curr_vs_VG.pdf}\includegraphics[width=0.49\textwidth]{Fig/tot_curr_vs_VG_log.pdf}}
\caption{Total current versus gate voltage at fixed source-drain bias in a linear scale (left) and a semi-logaritmic scale (right).}\label{FIG:tot_curr_vs_VG}
\end{figure}
\begin{figure}[ht]
\centerline{\includegraphics[width=0.49\textwidth]{Fig/den_1VG_04Vb.pdf}\includegraphics[width=0.49\textwidth]{Fig/den_m1VG_04Vb.pdf}}
\caption{Carrier densities at $V_b$ = 0.4 V. Left: $V_G$ = 1 V (case on). Right:  $V_G =  -$1 V (case off). }\label{FIG:densities}
\end{figure}
\begin{figure}[ht!]
\centerline{\includegraphics[width=0.49\textwidth]{Fig/pot_1VG_04Vb.pdf}\includegraphics[width=0.49\textwidth]{Fig/pot_m1VG_04Vb.pdf}}
\caption{Electrostatic potential with $V_{b}$ = 0.4 V. On the left the case of $V_G = 1$ V. On the right  the case of $V_G = - 1$ V.}\label{FIG:potential}
\end{figure}
The characteristic curves are shown in  Fig.s~\ref{FIG:tot_curr_vs_VG}.
It is evident that an acceptable field effect transistor is obtained. By decreasing the gate voltage the current settles at a very low value. The plot in logarithmic scale shows a current-on/current-off ratio of five orders of magnitudo at least.  The minority charges are only slightly activated and the inversion gate voltage observed with other geometries \cite{Sch} is practically not present. From  Fig.s~\ref{FIG:densities} it is evident that the hole density remains negligible with respects to the majority one. The main reason for such an effect is the peculiar position and breadth of the source and drain contacts which create an electrostatic potential such that the electron total energy is always of constant sign.
Devising a working structure has required a full 2D numerical solution of the Poisson equation. The common use of lumped 1D models for the electrostatic potential  has the drawback of hiding the effects related to the source and drain contact position and breadth. The simulations have been also performed adopting the  mobility model in \cite{Dor}.

The results are qualitatively  in good agreement with those obtained by using the  mobility model of Sec. \ref{section:mobility}. Therefore we are rather confident that the efficient FET performance of the proposed device is  not an artifact related to the
adopted specific mobility model but a general feature stemming from the chosen geometry. The proposed device overcomes the drawback related to the zero gap in pristine large area graphene and  represent a viable prototype for for the design of GFETs.
\section{Acknowledgements}
The authors  acknowledge the support from INdAM (GNFM). This work has been supported by the Universit\`{a} degli Studi di Catania, {\em Piano della Ricerca 2016/2018 Linea di intervento 2}.


text:
\section{Introduction}
\begin{wrapfigure}{r}{0.60\textwidth}
  \begin{center}
    \includegraphics[width=0.60\textwidth]{Figures/schematics.pdf}
  \end{center}
  \caption{A set of conferrable adversarial examples used as a fingerprint to identify surrogate models. \label{fig:schematics}}
\end{wrapfigure}

Deep neural network (DNN) classifiers have become indispensable tools for addressing practically relevant problems, such as autonomous driving~\citep{tian2018deeptest}, natural language processing \citep{young2018recent} and health care predictions~\citep{esteva2019guide}.
While a DNN provides substantial utility, training a DNN is costly because of data preparation (collection, organization, and cleaning) and computational resources required for validation of a model~\citep{press2016}.
For this reason, DNNs are often provided by a single entity and consumed by many, such as in the context of Machine Learning as a Service (MLaaS).
A threat to the provider is \emph{model stealing}, in which an adversary derives a \emph{surrogate model} from only API access to a \emph{source model}.
We refer to an independently trained model for the same task as a \emph{reference model}.

Consider a MLaaS provider that wants to protect their service and hence restrict its redistribution, e.g., through a contractual usage agreement because trained models constitute their intellectual property.
A threat to the model provider is an attacker who derives surrogate models and publicly deploys them.
Since access to the source model has to be provided, users cannot be prevented from deriving surrogate models.
\citet{krishna2019thieves} have shown that model stealing is (i) effective, because even high-fidelity surrogates of large models like BERT can be stolen, and (ii) efficient, because surrogate models can be derived for a fraction of the costs with limited access to domain data.

This paper proposes a DNN fingerprinting method to predict whether a model is a (stolen) surrogate or a (benign) reference model relative to a source model.
DNN fingerprinting is a new area of research that extracts a persistent, identifying code (fingerprint) from an already trained model.
Model stealing can be categorized into \emph{model modification}, such as weight pruning~\citep{zhu2017prune}, or \emph{model extraction} that uses some form of knowledge distillation~\citep{hinton2015distilling} to derive a surrogate from scratch.
Claimed security properties of existing defenses~(\citep{adi2018turning, zhang2018protecting}), have been broken by model extraction attacks~\citep{shafieinejad2019robustness}.
Our fingerprinting method is the first passive defense that is specifically designed towards withstanding model extraction attacks, which extends to robustness against model modification attacks.

Our research provides new insight into the transferability of adversarial examples.
In this paper, we hypothesize that there exists a subclass of targeted, transferable, adversarial examples that transfer exclusively to surrogate models, but not to reference models.
We call this subclass {\em conferrable}.
Any conferrable example found in the source model should have the same misclassification in a surrogate model, but a different one in reference models.
We propose a metric to measure conferrability and an ensemble adversarial attack that optimizes this new metric.
We generate conferrable examples as the source model's fingerprint.

Retrained CIFAR-10 surrogate models can be verified with a perfect ROC AUC of $1.0$ using our fingerprint, compared to an ROC AUC of $0.63$ for related work~\citep{cao2019ipguard}.
While our fingerprinting scheme is robust to almost all derivation and extraction attacks, we show that some adapted attacks may remove our fingerprint.
Specifically, our fingerprint is not robust to transfer learning when the attacker has access to a model pre-trained on ImageNet32 and access to CIFAR-10 domain data.
Our fingerprint is also not robust against adversarial training~\citep{madry2017towards} from scratch.
Adversarial training is an adapted model extraction attack specifically designed to limit the transferability of adversarial examples.
We hypothesize that incorporating adversarial training into the generation process of conferrable adversarial examples may lead to higher robustness against this attack.

\section{Related Work}

In black-box adversarial attacks~\citep{papernot2017practical, tramer2016stealing, madry2017towards}, access to the target model is limited, meaning that the target architecture is unknown and computing gradients directly is not possible.
Transfer-based adversarial attacks~\citep{papernot2016transferability, papernot2017practical} exploit the ability of an adversarial example to \emph{transfer} across models with similar decision boundaries.
\emph{Targeted} transferability additionally specifies the target class of the adversarial example.
Our proposed adversarial attack is a targeted, transfer-based attack with white-box access to a source model (that should be defended), but black-box access to the stolen model derived by the attacker.

\citet{liu2016delving} and \citet{tramer2017ensemble} show that (targeted) transferability can be boosted by optimizing over an ensemble of models.
Our attack also optimizes over an ensemble of models to maximize transferability to stolen surrogate models, while minimizing transferability to independently trained models, called reference models.
We refer to this special subclass of targeted transferability as \emph{conferrable}.
\citet{tramer2017ensemble} empirically study transferability and find that transferable adversarial examples are located in the intersection of high-dimensional "adversarial subspaces" across models.
We further their studies and show that (i) stolen models apprehend adversarial vulnerabilities from the source model and (ii) parts of these subspaces, in which conferrable examples are located, can be used in practice to predict whether a model has been stolen.

Watermarking of DNNs is a related method to DNN fingerprinting where an identifying code is embedded into a DNN, thereby potentially impacting the model's utility.
\citet{uchida2017embedding} embed a secret message into the source model's weight parameters, but require white-box access to the model's parameters for the watermark verification.
\citet{adi2018turning} and \citet{zhang2018protecting} propose backdooring the source model on a set of unrelated or slightly modified images.
Their approaches allow black-box verification that only requires API access to the watermarked model.
Frontier-Stitching~\citep{merrer2017adversarial} and BlackMarks~\citep{dong2018boosting} use (targeted) adversarial examples as watermarks.
These watermarks have been evaluated only against model modification attacks, but not against model extraction attacks that train a surrogate model from scratch.
At least two of these watermarking schemes~\citep{adi2018turning, zhang2018protecting} are not robust to model extraction attacks~\citep{shafieinejad2019robustness}.
\citet{cao2019ipguard} recently proposed a fingerprinting method with adversarial examples close to the source model's decision boundary.
We show that their fingerprint does not withstand retraining as a model extraction attack and propose a fingerprint with improved robustness to model extraction attacks.

\section{DNN Fingerprinting}
\label{sec:requirements}

\textbf{Threat Model. }
The attacker's goal is to derive a surrogate model from access to the defender's source model that (i) has a similar performance (measured by test accuracy) and (ii) is not verified as a surrogate of the source model by the defender.
We protect against an informed attacker that can have (i) white-box access to the source model, (ii) unbounded computational capacity and (iii) access to domain data from the same distribution.
A more informed attacker can drop information and invoke all attacks of a less informed attacker.
Robustness against a more informed attacker implies robustness against a less informed attacker.
Our attacker is limited in their access to ground-truth labeled data, otherwise they could train their own model and there would be no need to steal a model.
In our evaluation, we experiment with attackers that have up to 80\% of ground-truth labels for CIFAR-10~\citep{cifar10}.
A practical explanation for the limited accessibility of ground truth labels could be that attackers do not have access to a reliable oracle (e.g., for medical applications), or acquiring labels may be associated with high expenses (e.g., Amazon Mechanical Turk\footnote{\url{https://www.mturk.com/}}).

The defender's goal is to identify stolen surrogate models that are remotely deployed by the attacker.
Their capabilities are (i) white-box access to the source model, (ii) black-box access to the target model deployed by the attacker and (iii) a limited set of $n$ queries to verify a remotely deployed surrogate model.
Black-box access to the suspect model excludes knowledge of the suspect model's architecture or attack used to distill the model.
The defender also does not have knowledge of, nor control over the attacker's dataset used to steal the model.

\textbf{Fingerprinting Definitions. }
A fingerprinting method for DNNs consists of two algorithms: (i) A fingerprint \emph{generation} algorithm that generates a secret fingerprint $\mathcal{F} \subseteq \mathcal{X}^n$ of size $n$, and a fingerprint verification key $\mathcal{F}_y \subseteq \mathcal{Y}^n$; (ii) A fingerprint \emph{verification} algorithm that verifies surrogates of the source model.
These algorithms can be summarized as follows.
\begin{itemize}
	\item \textbf{Generate}$(M, D)$: Given white-box access to a source model $M$ and training data $D\in \mathcal{D}$. Outputs a fingerprint $\mathcal{F}$ and the verification keys $\mathcal{F}_y=\{M(x) | x \in \mathcal{F}\} $.
	\item $\textbf{Verify}(\hat{M}(\mathcal{F}), \mathcal{F}_y)$: Given black-box access to a suspect model $\hat{M}$, a fingerprint $\mathcal{F}$ and a verification key $\mathcal{F}_y$. Outputs $1$ if $\hat{M}$ is verified by the fingerprint and $0$ otherwise.
\end{itemize}
The verification algorithm computes an error rate between the outputs of the source and target model on the fingerprint.
We empirically measure the error rate separately for surrogate and reference models, which allows to choose a decision threshold $\rho\in[0,1]$.
If the error-rate of a target model exceeds $1-\rho$, the verification predicts the target model to be a reference model, otherwise the prediction is a surrogate model.
We define that a fingerprint must be \emph{irremovable}, i.e., surrogate and reference models are correctly verified, despite an attacker's removal attempt.
It must also be \emph{non-evasive}, meaning that an attacker cannot evade black-box verification by detecting fingerprint queries.
We refer to Appendix~\ref{sec:appendix_securitygames} for security games of DNN fingerprinting.
\begin{figure*}
	\centering
	\subfloat[]{\raisebox{-0.5\height}{\includegraphics[width=.33\linewidth]{Figures/venn_svg.pdf}}\label{fig:intuition1}}\hspace{5ex}
	\subfloat[]{\raisebox{-0.5\height}{\includegraphics[width=.55\linewidth]{Figures/confer_graphics_updated.pdf}}\label{fig:intuition2}}
	\caption{(a) A summary of the relationship between transferability and conferrability as intersections between the set of all adversarial examples per model type. (b) A representation of transferable and conferrable examples in the decision space, relative to the ground-truth provided by an oracle.\label{fig:intuition}}
\end{figure*}
\section{Conferrable Adversarial Examples}

\textbf{Motivating Conferrability}
Conferrability is a new property for adversarial examples, in which  targeted transferability occurs only from a source model to its surrogates, but not to independently trained reference models.
Intuitively, surrogate models are expected to be more similar to the source model than any reference model, but quantifying this similarity is non-trivial.
Conferrable examples are an attempt to quantify this similarity by shared adversarial vulnerabilities of the source model and its surrogates.
Surrogate models differ from reference models in the objective function which they optimize.
Reference models maximize fidelity to the ground truth labels, whereas surrogate models maximize fidelity to the source model's labels, which do not always coincide.

Fig.~\ref{fig:intuition1} shows the relation between targeted adversarial, transferable, and conferrable examples for the source model, its surrogates, and all reference models.
An example is transferable if it is adversarial to any model.
It is conferrable when it is adversarial only to surrogate and source models.
In that sense, conferrable adversarial examples are a subclass of targeted transferable adversarial examples.
Fig.~\ref{fig:intuition2} shows the relation between transferable and conferrable adversarial examples in the decision space, simplified for illustrative purposes to binary classification.
Transferable adversarial examples occur in those adversarial subspaces where the decision boundary of surrogate and reference models coincide~\citep{tramer2017space}.
Conferrable adversarial examples occur in those adversarial subspaces where the decision boundary of surrogate and reference models differ.

Targeted transferability for a class $t$ and a set of models $\mathcal{M}$ can be computed as follows.
\begin{align}
     \text{Transfer}(\mathcal{M}, x;t) = \underset{M\in \mathcal{M}}{Pr}[M(x)=t]
\end{align}

\textbf{Objective} We want to find adversarial examples that maximize the output activation difference between surrogate models $\mathcal{S}$ and reference models $\mathcal{R}$.
This difference is quantified by our \emph{conferrability} score, which measures the example's transferability to surrogate and reference models.
\begin{align}
\label{eq:confer}
\text{Confer}(\mathcal{S}, \mathcal{R}, x;t) &= \text{Transfer}(\mathcal{S}, x;t) (1-\text{Transfer}(\mathcal{R}, x;t))
\end{align}
The central challenge for optimizing conferrability directly is that we require access to a differentiable function \textit{Transfer} that estimates an example's transferability score.
To the best of our knowledge, the only known method to evaluate transferability is to train a representative set of DNNs and estimating the example's transferability itself.
We use this method to evaluate transferability.
In the case of conferrable examples, we evaluate transferability on a set of surrogate and reference models trained locally by the defender and use Equation~\ref{eq:confer} to obtain the conferrability score.
Our hypothesis is that conferrability generalizes, i.e., examples that are conferrable to a representative set of surrogate and reference DNNs are also conferrable to other, unseen DNNs.

\textbf{Conferrable Ensemble Method. }
In this section, we describe our adversarial attack that generates conferrable adversarial examples, called the \emph{Conferrable Ensemble Method} (CEM).
CEM operates on the source model $M$, a set of its surrogates $\mathcal{S}_M$, and a set of reference models $\mathcal{R}$.
The attack constructs an ensemble model $M_E$ with a single shared input that outputs conferrability scores for all output classes $y\in \mathcal{Y}$.
CEM generates highly conferrable adversarial examples by finding a perturbation $\delta$ so that $x'=x_0+\delta$ maximizes the output of the ensemble model for a target class $t$.

We now present the construction of the ensemble model $M_E$.
The ensemble model produces two intermediate outputs for an input $x\in \mathcal{X}$, representing the average predictions of all surrogate and all reference models on $x$.
We use Dropout~\citep{srivastava2014dropout} with a drop ratio of $d=0.3$.
\begin{align}
\text{Surr}(\mathcal{S}_M, x) &= \frac{1}{|\mathcal{S}_M|} \sum_{S\in \mathcal{S}_M} \text{Dropout}(S(x); d) \\
\text{Ref}(\mathcal{R}, x) &= \frac{1}{|\mathcal{R}|} \sum_{R \in \mathcal{R}} \text{Dropout}(R(x); d)
\end{align}

The output of the ensemble model is the conferrability score of the input for each class.
The Softmax activation function is denoted by $\sigma$.
We refer to Appendix~\ref{sec:appendix_conferrabilityscore} for details about the optimization.
\begin{align}
M_E(x; \mathcal{S}_M, \mathcal{R}) = \sigma(\text{Surr}(\mathcal{S}_M, x)(\mathbf{1}-\text{Ref}(\mathcal{R}, x)))
\end{align}
The loss function consists of three summands and is computed over the benign, initial example $x_0$, and the example at an intermediate iteration step $x'=x_0+\delta$.
We denote the cross-entropy loss by $H(\cdot, \cdot)$.
The first summand function maximizes the output of the ensemble model for some target class $t$.
The second summand maximizes the categorical cross-entropy between the current and initial prediction for the source model.
The third summand minimizes the categorical cross-entropy between the source model's prediction and the predictions of its surrogates.
The total loss $\mathcal{L}$ is the weighted sum over all individual losses.
\begin{align}
\label{eq:optimizationtarget}
\mathcal{L}(x_0, x') &= \alpha H(1, \max_t [M_E(x')_t])-\beta H(M(x_0), M(x')) + \gamma H(M(x'), \text{Surr}(\mathcal{S}_M, x'))
\end{align}
In all our experiments, we use weights $\alpha=\beta=\gamma=1$ and refer to Appendix~\ref{sec:appendix_sensitivty} for an empirical sensitivity analysis of the hyperparameters.
We address the box constraint, i.e., $||\delta|| \leq \epsilon$, similarly to PGD~\citep{madry2017towards} by clipping intermediate outputs around the $\epsilon$ ball of the original input $x_0$.
For the optimization of an input with respect to the loss we use the Adam optimizer~\citep{kingma2014adam}.
Note that we use an untargeted attack to generate targeted adversarial examples.
Targeted attacks are harder to optimize than their untargeted counterparts~\citep{wu2019universal}.
We assign the source model's predicted label for the generated adversarial example as the target label and ensure that target classes are balanced in the fingerprint verification key.

\textbf{Fingerprinting Algorithms. }
We now describe our fingerprinting generation and verification algorithms.
For the generation, the defender locally trains a set of $c_1$ surrogate and $c_2$ reference models ($c_1=c_2=18$) on their training data prior to executing CEM.
Surrogate models are trained on data labeled by the source model, whereas reference models are trained on ground-truth labels.
The defender composes the ensemble model $M_E$ as described in the previous paragraph and optimizes for a perturbation $\delta$ given an input $x_0$, so that $\mathcal{L}(x_0, x_0+\delta)$ is minimized.
The optimization returns a set of adversarial examples that are filtered by their conferrability scores on the locally trained models.
If their conferrability score (see Equation~\ref{eq:confer}) exceeds a minimum threshold ($\tau \geq 0.95$), they are added to the fingerprint.
For the fingerprint verification, we compute the error rate between the source model's prediction on the fingerprint and a target model's predictions.
If the error rate is greater than $1-\rho$, which we refer to as the \emph{decision threshold}, the target model is predicted to be a reference model and a surrogate model otherwise.

\section{Experimental Setup}
\label{sec:experimental_setup}
\textbf{Experimental Procedure.} We evaluate the irremovability and non-evasiveness of our fingerprint.
We show that our proposed adversarial attack CEM improves upon the mean conferrability scores compared to other adversarial attacks, such as FGM~\citep{goodfellow2014explaining}, BIM~\citep{kurakin2016adversarial}, PGD~\citep{madry2017towards}, and CW-$L_\infty$~\citep{carlini2017towards}.
We demonstrate the non-evasiveness of our fingerprint by evaluating against the evasion algorithm proposed by \citet{hitaj2019evasion}.
Our study on irremovability is the most extensive study conducted on DNN fingerprints or DNN watermarks compared to related work~\citep{cao2019ipguard, uchida2017embedding, zhang2018protecting, adi2018turning, merrer2017adversarial, dong2018boosting, szyller2019dawn}.
We compare our DNN fingerprinting to another proposed DNN fingerprint called IPGuard~\citep{cao2019ipguard}.
We refer to Appendix~\ref{sec:cifar10_additional} for more details on the CIFAR-10 experiments and to Appendix~\ref{sec:imagenet32exp} for results on ImageNet32.

\textbf{Evaluation Metrics. }
We measure the fingerprint retention as the success rate of conferrable examples in a target model $\hat{M}$.
An example is successful if the target model predicts the target label given in $\mathcal{F}_y$.
We refer to this metric as the \emph{ Conferrable Adversarial Example Accuracy} (CAEAcc).
\begin{align}
    \text{CAEAcc}(\hat{M}(\mathcal{F}), \mathcal{F}_y) = \underset{(y, y^*) \in \hat{M}(\mathcal{F}), \mathcal{F}_y}{Pr}[\bm{1}_\mathrm{y=y^*}]
\end{align}

\textbf{CIFAR-10 models.}
We train popular CNN architectures on CIFAR-10 without any modification to the standard training process\footnote{\url{https://keras.io/examples/cifar10_resnet}}.
The defender's source model has a ResNet20~\citep{resnet} architecture.
All surrogate and reference models required by the defender for CEM are ResNet20 models trained on CIFAR-10 using retraining.
The attacker has access to various model architectures, such as DenseNet~\citep{densenet}, VGG-16/VGG-19~\citep{vgg} and ResNet20~\citep{resnet}, to show that conferrability is maintained across model architectures.

\textbf{Removal Attacks.}
Removal attacks are successful if (i) the surrogate model has a high test accuracy (at least $85.55\%$ for CIFAR-10, see Fig.~\ref{fig:caevstest}) and (ii) the stolen surrogate's CAEAcc is lower than the verification threshold $\rho$.
Stolen models are derived with a wide range of fingerprint removal attacks, which we categorize into (i) model modification, (ii) model extraction and (iii) adapted model extraction attacks.
Model modification attacks modify the source model directly, such as fine-tuning or parameter pruning.
Model extraction distills knowledge from a source model into a fresh surrogate model, such as retraining or model extraction attacks from related work~\citep{papernot2017practical, jagielski2019high, orekondy2019knockoff}.
A defense against model extraction attacks is to restrict the source model's output to the top-1 predicted label.
We refer to an extraction attack that retrains on the top-1 predicted label by the source model as \textit{Black-Box} attack.
We also evaluate transfer learning as a model extraction attack, where a pre-trained model from a different domain (ImageNet32) is transfer learned onto the source model's domain using the source model's labels.

Adapted model extraction attacks limit the transferability of adversarial examples in the surrogate model, such as adversarial training~\citep{madry2017towards}.
In our experiment, we use multi-step adversarial training with PGD.
We design another adapted attacks against our fingerprint called the \textit{Ground-Truth} attack, where the attacker has a fraction of $p\in[0.6, 0.7, 0.8]$ ground-truth labels.
We trained surrogate models with Differential Privacy (DP), using DP-SGD~\citep{abadi2016deep}, but even for large $\epsilon$ the surrogate models had unacceptably low CIFAR-10 test accuracy of about $76\%$.
The test accuracies of all surrogate models can be seen in Table~\ref{tab:test_acc_cifar10}.
We repeat all removal attacks three times and report the mean values and the standard deviation as error bars.
Our empirically chosen CIFAR-10 decision thresholds $\theta_\epsilon$ (see Table~\ref{tab:thethaeps}) are chosen relative to the perturbation threshold $\epsilon$ with which our adversarial attacks are instantiated.

\begin{table}
	\centering
	\begin{tabular}{|c||c|c|c|c|c|c|c|}
		\hline
		$\epsilon$ &  $0.01$ &$0.025$ & $0.05$ & $0.075$ & $0.1$ & $0.125$ & $0.15$ \\
		\hline
		$\theta_\epsilon$ & 63.00\% & 75.00\% & 84.00\% & 86.00\% & 87.00\% & 87.00\% & 87.00\% \\
		\hline
	\end{tabular}
	\caption{Experimentally derived values for the decision threshold $\theta_\epsilon$ for CIFAR-10. \label{tab:thethaeps}}
\end{table}

\textbf{Evasiveness Attack. } \citet{hitaj2019evasion} present an evasion attack against black-box watermark verification. Their evasion attack trains a binary DNN classifier to distinguish between benign and out-of-distribution images and reject queries by the defender.
Our attacker trains a binary classifier to classify benign images and adversarial
examples generated by FGM, CW-$L_\infty$, and PGD adversarial attacks.
The evasion attack is evaluated
by the area under the curve (AUC) of the receiver operating characteristic (ROC) curve for varying
perturbation thresholds.

\textbf{Attacker Datasets. } We distinguish between attackers with access to different datasets, which are CIFAR-10~\citep{cifar10}, CINIC~\citep{darlow2018cinic} and ImageNet32~\citep{chrabaszcz2017downsampled}.
CINIC is an extension of CIFAR-10 with downsampled images from ImageNet~\citep{deng2009imagenet}, for classes that are also defined in CIFAR-10.
ImageNet32 is a downsampled version of images from all classes defined in ImageNet, i.e., ImageNet32 is most dissimilar to CIFAR-10.
We observe that the similarity of the attacker's to the defender's dataset positively correlates with the surrogate model's test accuracy and thus boosts the effectiveness of the model stealing attack.

\section{Empirical Results}

\begin{figure*}
	\centering

	\subfloat[]{\includegraphics[trim=17 5 43 19, clip,width=.329\textwidth]{Figures/conf_scores.pdf} \label{fig:confscores}} \hfill
	\subfloat[]{\includegraphics[trim=18 5 43 19, clip,width=.329\textwidth]{Figures/detectability_roc.pdf} \label{fig:detectability}} \hfill

	\subfloat[]{\includegraphics[trim=10 5 43 19, clip, height=3.25cm, width=.329\textwidth]{Figures/fine_tuning.pdf}\label{fig:fine_tuning}} \hfill

	\subfloat[]{\includegraphics[trim=17 5 43 19, clip,width=.33\textwidth]{Figures/cifar_pruning.pdf}\label{fig:pruning}} \hfil
	\subfloat[]{\includegraphics[trim=17 5 43 19, clip,width=.33\textwidth]{Figures/fingerprint_retention_cifar10.pdf}\label{fig:retrain}} \hfill
	\subfloat[]{\includegraphics[trim=17 5 43 19, clip,width=.33\textwidth]{Figures/diff_architecture.pdf}\label{fig:diffarch}}\hfill
	\caption{
	Results for the (a) conferrability scores, (b) non-evasiveness and (c-d) irremovability to model modification attacks.
	Figures (e,f) show results against model extraction attacks.
	\label{fig:res_cifar10_row1}
	\label{fig:res_cifar10_row2}
	}
\end{figure*}
\textbf{Conferrability Scores.} Our empirical results show support in favor of the existence of conferrable adversarial examples.
CEM produces adversarial example with the highest conferrability scores compared to classical adversarial attacks such as FGM, PGD and CW-$L_\infty$.
We generate $n=100$ adversarial examples from the same set of CIFAR-10 test samples, filter those that are non-successful to the source model and evaluate conferrability scores as specified by Equation~\ref{eq:confer} across five unseen CIFAR-10 surrogate and reference models trained from scratch.
The results in Fig.~\ref{fig:confscores} show that CEM outputs adversarial examples with significantly higher conferrability scores than other attacks for small perturbation thresholds $\epsilon$, which worsen as $\epsilon$ increases.
For $\epsilon=0.01$ we measure a mean conferrability score of $0.49$, which amounts to a mean CAEAcc of $0.42$ in reference models and $0.85$ in surrogate models (see Fig.~\ref{fig:pruning}).
For $\epsilon=0.15$ we measure a conferrability score of only $0.24$, which translates to a mean CAEAcc of $0.98$ in surrogate models and $0.75$ in reference models.
These results show that with increasing $\epsilon$, the transferability of generated adversarial examples improves at the expense of lower conferrability.
We observed that the conferrability scores measured on surrogate and reference models used in CEM are nearly perfect.
This indicates that CEM could benefit from access to an even larger ensemble of models.

\textbf{Non-Evasiveness.} We evaluate non-evasiveness against the attack proposed by \citet{hitaj2019evasion}.
The results in Fig.~\ref{fig:detectability} show that (i) detectability increases with larger perturbation thresholds $\epsilon$ and (ii) for small perturbation thresholds $\epsilon \leq 0.025$ we measure a ROC AUC of only $0.67$.
This method of detection is too unreliable for the attacker to deploy in practice considering that only a fraction of requests would contain the fingerprint. Sufficiently high recall against our fingerprint comes at the cost of more false positives, which diminishes the utility of the attacker's model to other users.

\textbf{Model Modification Attacks.} Fig.~\ref{fig:fine_tuning} and \ref{fig:pruning} show that our fingerprint is robust against model modification attacks.
We show irremovability against four different types of fine-tuning attacks tried by \citet{uchida2017embedding}, which are implemented as follows.
\begin{enumerate}
	\item Fine-Tune Last Layer (\textbf{FTLL}): Freezes all layers except for the penultimate layer.
	\item Fine-Tune All Layers (\textbf{FTAL}): Fine-tuning of all layers.
	\item Retrain Last Layer (\textbf{RTLL}): Re-initializes the penultimate layer's weights and updates only the penultimate layer, while all other layers are frozen.
	\item Retrain All Layers (\textbf{RTAL}): Re-initializes the penultimate layer's weights, but all layers are updated during fine-tuning.
\end{enumerate}
Our results for iterative weight pruning~\citep{zhu2017prune} show that larger pruning rates $p\in[0.7,0.8,0.9]$ have greater impact on the CAEAcc, but the surrogate test accuracy also deteriorates significantly (see Table~\ref{tab:test_acc_cifar10}).
Note that we show robustness to much higher pruning rates than related work~\citep{cao2019ipguard}.
The attack for $p=0.9$ is unsuccessful because the surrogate test accuracy is lower than $85.55\%$, as described in Section~\ref{sec:experimental_setup}. For all remaining pruning configurations the fingerprint is not removable.

\textbf{Model Extraction Attacks.} Fig.~\ref{fig:res_cifar10_row2} and \ref{fig:experiments_results_2} show that our fingerprint is robust to model extraction attacks, except for transfer learning when the attacker has CIFAR-10 data.
Surprisingly, we find that surrogates extracted using CINIC have higher mean CAEAcc values than CIFAR-10 surrogates at a lower test accuracy ($-2.16\%$).
Similarly, we measure significantly lower CAEAcc values for reference models trained on CINIC compared to reference models trained on CIFAR-10.
This means that it is increasingly difficult for the attacker to remove our fingerprint the more dissimilar the attacker's and defender's datasets become.
In Fig.~\ref{fig:diffarch}, we show that our fingerprint is not removable even across different surrogate model architectures.
Fig.~\ref{fig:transferlearn} shows that transfer learning removes our fingerprint when the attacker has access to the defender's dataset, but our fingerprint is not removable when the attacker only has access to CINIC data.
We observed that the pre-trained ImageNet32 model exceeds 80\% test accuracy after only a single epoch with CIFAR-10, which may be too few updates for our fingerprint to transfer to the surrogate.
The irremovability of our fingerprint to a wide range of model extraction attacks shows that certain adversarial vulnerabilities (that enable conferrable examples) are consistently carried over from the source model to its surrogates.
Our data supports that the class of transferable examples can be broken down further into conferrable examples and that known findings for transferable examples extend to conferrable examples.

\begin{figure*}

	\subfloat[]{\includegraphics[trim=18 5 43 20, clip,width=.33\textwidth]{Figures/model_extraction.pdf}\label{fig:modelextraction}} \hfill
	\subfloat[]{\includegraphics[width=.33\textwidth, height=3.75cm]{Figures/transfer_learning.pdf}\label{fig:transferlearn}} \hfill

	\subfloat[]{\includegraphics[trim=17 5 43 20, clip,width=.33\textwidth]{Figures/adv_train.pdf}\label{fig:advtrain}} \hfill

	\subfloat[]{\includegraphics[trim=18 5 43 20, clip,width=.33\textwidth]{Figures/true_labels.pdf}\label{fig:truelabels}} \hfill
	\subfloat[]{\includegraphics[trim=17 5 38 10, clip, height=3.4cm, width=.33\textwidth]{Figures/cae_vs_test_both.pdf}\label{fig:caevstest}} \hfill
	\subfloat[]{\includegraphics[trim=17 5 43 20, clip,width=.33\textwidth]{Figures/fp_eval.pdf}\label{fig:fpeval}} \hfill

	\caption{Fig.~(a-d) show the robustness of our fingerprint against (adapted) model extraction attacks. Fig.~(e) shows the test accuracy in relation to the CAEAcc for surrogate and reference models during training. Fig.~(f) shows the ROC curve of our fingerprint compared to IPGuard~\citep{cao2019ipguard}.}
	\label{fig:experiments_results_2}
\end{figure*}

\textbf{Adapted Model Extraction Attacks.} We now show limitations of our fingerprint against defenses that specifically limit the transferability of adversarial examples.
Fig.~\ref{fig:advtrain} shows that adversarial training from scratch removes our fingerprint.
In CIFAR-10 surrogates, we measure a mean CAEAcc of only $15\%$ for $\epsilon=0.025$, which supports the claims by \citet{madry2017towards} that adversarial training increases robustness to transfer attacks.
We hypothesize that incorporating adversarial training into the generation process of conferrable
adversarial examples leads to higher robustness against this removal
attack.
The results from our Ground-Truth attack shown in Fig.~\ref{fig:truelabels} show that access to more ground-truth labels decreases the CAEAcc in surrogate models.
The results show that our fingerprint is not removable for attackers with up to $50\%$ ground-truth labels for CIFAR-10.

\textbf{Confidence Analysis.} In Fig.~\ref{fig:fpeval} we show the ROC curve (false positive vs true positive rate) of our fingerprint verification on ten unseen, retrained CIFAR-10 surrogate and reference models.
We compare our work with IPGuard~\citep{cao2019ipguard} and generate $n=100$ fingerprints in both cases.
Our ROC AUC is $1.0$, whereas IPGuard has a ROC AUC of only $0.63$.
These results show that IPGuard is not robust to retraining as a model extraction attack.
Inspecting adversarial examples generated by IPGuard further, we measure a mean adversarial success rate of $17.61\%$ for surrogate and $16.29\%$ for reference models.
For our fingerprint, Fig.~\ref{fig:caevstest} visualizes the difference in CAEAcc for well-trained surrogate and reference models at $\epsilon=0.025$.
The plot shows that CAEAcc positively correlates with the surrogate's test accuracy.
We measure a mean difference in CAEAcc of about $30\%$ between well-trained surrogate and reference models with our fingerprint.

\section{Conclusion}
We empirically show the existence of conferrable adversarial examples.
Our ensemble adversarial attack CEM outperforms existing adversarial attacks such as FGM~\citep{goodfellow2014explaining}, PGD~\citep{madry2017towards} and CW-$L_\infty$~\citep{carlini2017towards} in producing highly conferrable adversarial examples.
We formally define fingerprinting for DNN classifiers and use the generated conferrable adversarial examples as our fingerprint.
Our experiments on the robustness of our fingerprint show increased robustness to model modification attacks and model extraction extraction attacks.
Transfer learning is a successful removal attack when the attacker has access to CIFAR-10 data, but not when the attacker only has access to CINIC.
Adversarial training from scratch is the most effective removal attack and successfully removes our fingerprint.
We hypothesize that adding adversarial training into the generation process of conferrable adversarial examples may increase robustness against adversarial training.
Our experiments confirm the non-evasiveness of our fingerprint against a detection method proposed by \citet{hitaj2019evasion}.
We empirically find that our fingerprint is the first to perfectly verify retrained CIFAR-10 surrogates with a ROC AUC of 1.0.

\begin{table}[]
\begin{tabular}{|c|c|c|c|c|c|c|c|}
\hline
                   & \multicolumn{4}{c|}{\textbf{Fine-Tuning}}                                     & \multicolumn{3}{c|}{\textbf{Pruning}}    \\ \hline
Source             & FTLL               & FTAL                     & RTLL                 & RTAL   & p=0.7        & p=0.8       & p=0.9       \\ \hline
89.30              & 89.59              & 89.57                    & 87.90                & 87.25  & 88.00        & 86.97       & 83.36       \\ \hline
\multicolumn{5}{|c|}{\textbf{Retrain}}                                                             & \multicolumn{3}{c|}{\textbf{Extraction}} \\ \hline
ResNet20          & Densenet           & VGG16                    & VGG19                & CINIC  & Jagielski    & Papernot    & Knockoff    \\ \hline
89.22              & 90.88              & 90.91                    & 90.05                & 87.06  & 88.74        & 87.34       & 84.55       \\ \hline
\multicolumn{2}{|c|}{\textbf{AdvTrain}} & \multicolumn{2}{c|}{\textbf{Transfer Learning}} & \multicolumn{4}{c|}{\textbf{Ground-Truth}}        \\ \hline
CIFAR-10           & CINIC              & CIFAR-10                 & CINIC                & p=0.5  & p=0.6        & p=0.7       & p=0.8       \\ \hline
89.75              & 88.15              & 89.59                    & 88.15                & 89.54  & 89.13        & 89.43       & 89.33       \\ \hline
\end{tabular}
\caption{Mean CIFAR-10 test accuracies for surrogate models obtained by running the attacks with threefold repetition. Unless stated otherwise, all surrogates are trained on CIFAR-10 images.
\label{tab:test_acc_cifar10}}
\end{table}

redpajama_set_name: RedPajamaArXiv


text:
\section{
Introduction
				}
Elliptical galaxies have been known as X-ray sources since
several were detected in the Virgo cluster by the
$Einstein$ $Observatory$ (Forman et al. 1979).
Long and Van Speybroeck (1983) compiled many more $Einstein$
detections of ellipticals and Nulsen, Stewart and Fabian (1984)
showed that diffuse gas is their dominate source of X-ray emission.
Analysis of $Einstein$ $Observatory$ IPC spectra shows that the emission
from luminous ellipticals can be
fit by thermal plasma models with fairly low temperatures, 0.5 --
2.0 keV (Forman, Jones, \& Tucker 1985; Trinchieri, Fabbiano \& Canizares 1986).
More recent $Ginga$, $ROSAT$ and $ASCA$ observations of individual galaxies
confirm these characteristic temperatures (cf. Awaki et al.1991;
Forman et al. 1993; Fabbiano, Kim \& Trinchieri 1994; Awaki et al.1994;
Matsushita et al.1994).
The emission from lower luminosity ellipticals is apparently more complex,
with evidence for both harder and ``supersoft" emission according to
Fabbiano, Kim \& Trinchieri (1992); these authors attribute
the hard component to stellar X-ray binaries and suggest the supersoft emission
may be from M stars, ultrasoft X-ray binaries and/or relatively cool ISM.

To date, there has been no study of the X-ray properties of a $complete$
sample of normal elliptical galaxies.  The IPC detector on the $Einstein$
$Observatory$ was used to observe $\sim200$ normal galaxies,
including many ellipticals (see the review by Fabbiano 1989), but these data
allow only small ($N\approx14$) complete samples of ellipticals to be
constructed.
Furthermore, the limited energy resolution of the IPC prevented accurate
temperature determinations.
We have therefore turned to the longer exposures of $ROSAT$ pointed
observations to construct the largest possible complete sample of accurate
X-ray luminosities and temperatures for elliptical galaxies.
In this letter we
present new gas temperature and abundance determinations from
$ROSAT$ PSPC spectra of a subset of this complete sample
and explore their relation to optical and dynamical
properties of the host galaxies.  We will present detailed analyses
of their X-ray luminosities in a forthcoming paper (Davis \& White 1996).

				\section{
Sample Selection
				}
We used the $HEASARC$ $BROWSE$ software to select our optical magnitude-limited
sample from the $Third$ $Revised$ $Catalog$ $of$ $Galaxies$ (de Vaucouleurs
et al. 1991; hereafter {\sl RC3}), using its corrected $B$ magnitudes.
We limited our sample to elliptical galaxies (de Vaucouleurs T-types -4, -5
\& -6).
We then used $BROWSE$ to cross-correlate this sample with the $ROSAT$ and
$Einstein$ pointed observations.  Neglecting galaxies in the zone-of avoidance
(taken here to be $|b|<20^\circ$), we found that 38 of the 43 optically
brightest ellipticals have $ROSAT$ PSPC data, and 4 of the remaining 5 have
$Einstein$ IPC data (the remaining galaxy, NGC 7507, will soon have $ROSAT$ HRI
data available, but this has no spectral information);  since $ROSAT$ PSPC data
have better spectral resolution than $Einstein$ IPC data, they are preferable
when there is a choice.
The sample of 43 optically-selected ellipticals is complete to a limiting
(corrected) blue magnitude of $m_B^0=11.36$ and includes
four Local Group dwarfs.
The X-ray counterpart to our optical sample is thus 98\% complete, with X-ray
data for 42 of the 43 ellipticals.
If we press to fainter optical magnitudes, the marginal completeness in the
X-ray data deteriorates at the 50\% level for the next 10 or so optically
brightest galaxies.
In this letter we present a subset of 30 galaxies from this complete sample
which have sufficient spectral counts to accurately determine their X-ray
temperatures; 27 of these galaxies also have well-constrained abundances.

				\section{
X-Ray Data
				}
In order to tie the X-ray properties of ellipticals to their optical properties
in a uniform manner, we let the optical photometry dictate the spatial size of
the X-ray spectral extraction regions.
For each galaxy in our optically selected sample, we extract X-ray spectra from
within $6r_e$, where $r_e$ is the galaxy's optical effective radius
(derived from the {\sl RC3}).
Due to X-ray scattering in the $ROSAT$ mirror assembly, the minimum
suggested extraction radius for point sources is $1.5^\prime$, which is the
99.5\% encircled energy radius for a 1 keV spectrum.
Our minimum metric extraction radius must be $3r_e$ to meet this requirement
with the current sample.
The purpose of our larger metric extraction radius of $6r_e$ is to provide
more global X-ray properties, for better comparison with global
optical quantities such as magnitudes
($6r_e$ encloses 91\% of a de Vaucouleurs $r^{1/4}$ profile).

For the PSPC data, source spectra were selected from observation times when the
{\it Master Veto Rate} ($MVR$) was within recommended bounds: $20<MVR<170$.
The background for each spectrum was obtained from an annulus around
the galaxy extending from $6r_e$ to $12r_e$ from its center.
All detectable point sources were excluded from both source and background
spectra;
excluded points in the source regions tended to be in the outskirts so they
are unlikely to be X-ray binaries. Source and background spectra were also
corrected for vignetting effects and for the residual
particle background (in the manner prescribed by Snowden et al. 1994),
which is proportional to the $MVR$.
To ensure correct Poisson errors were used throughout, we modified the
appropriate data header as outlined in the $PROS$ QPSPEC documentation.
Finally, the extracted spectra were rebinned so that each channel has a
minimum of 25 counts, and systematic errors of 1\% were included.

				\section{
Analysis
				}
We used {\sl XSPEC~9.0} software (Schafer et al. 1990) to fit various models to
the extracted spectra.
An isothermal Raymond-Smith plasma model was our benchmark model,
with the temperature and abundance allowed to vary in the fits.
We also included a variable absorption component due to the column density of
Galactic hydrogen in the line-of-sight.
The redshift of the model spectrum was fixed to that corresponding to the
heliocentric velocity of the galaxy given in the {\sl RC3}.
The extracted spectra were fit between $\sim$0.2 and 2.0 keV, the exact energy
boundaries being set by the channel grouping.
Once a minimum in $\chi^2$ was found, the 90\% confidence
errors were determined for the free parameters.
Most of the galaxies have sufficient
counts so that the temperature, abundance and absorption could all be fit.
For a few galaxies, however, the number of free parameters had to be reduced to
allow reasonable errors to be determined.
For such galaxies we first fixed the abundance to 0.2
solar; when necessary, we also fixed the absorption to the Galactic value
(Stark et al. 1992).
The Raymond-Smith thermal model provided formally excellent fits to all
but two galaxies, with $\chi^2$ per degree of freedom $\nu$ of
$\chi^2_\nu\approx1$.
The two exceptions (NGC~4486 \& IC~1459) are known to have AGN,
so we fit composite models consisting of a Raymond-Smith isothermal
plus a power-law spectral component; these composite spectra achieved
$\chi^2_\nu\approx1$.
The detailed spectral fitting results are shown in Table~1, where the
indicated errors are 90\% confidence limits.

				\subsection{
$kT$--$\sigma$ Correlation
				}

For the 30 galaxies which have enough counts to allow accurate
temperature determinations, Fig.~1 shows a logarithmic plot of their
temperatures $kT$ against $\sigma$, the central
line-of-sight stellar velocity dispersions within the galaxies (from Faber
et al. 1989 or Dressler, Faber \& Burstein 1991).
Characteristics of the best fit linear regression to these data are
shown in the first line of Table~2.  We found that three galaxies (M32,
NGC 1052 and NGC 2768) have temperatures of $\sim3-4$ keV, which are
much greater than those of other galaxies with similar velocity dispersions.
These galaxies are likely to have their X-ray emission dominated by X-ray
binaries or AGN (cf. Eskridge, White \& Davis 1996;
Ho, Filipenko \& Sargent 1995), rather than by the cooler diffuse gas
expected from the accumulation of stellar ejecta.
Neglecting these three galaxies,
the best fit linear regression yields $kT \propto \sigma^{1.45\pm 0.20}$ and
the temperature dispersion is only 0.10 dex (see line 2 of Table~2);
the regression is shown as a solid line in Figure~1.
This correlation is shallower than what would obtain if the gas
and stars had the same kinetic temperature: $kT=\mu m_p\sigma^2$; this latter
relation is also shown in Fig.~1 as a dashed line.
The bulk of our temperature measurements indicate that warm diffuse gas
is being detected, even down to low luminosities and velocity dispersions
where some investigators expect such gas to be swept by winds (e.g. David,
Forman \& Jones 1991; Ciotti et al. 1991).
Fig.~1 also shows that the global gas temperatures are hotter than the kinetic
temperatures of luminous stars in the galaxies, since $kT>\mu m_p \sigma^2$
in all the galaxies.

				\subsection{
Temperature--Abundance Correlation
				}
For those 27 galaxies with enough counts for elemental abundances to be
constrained, we find that the measured abundances (driven by $\sim$1 keV Fe L
emission in $ROSAT$ PSPC spectra) are correlated with the gas temperatures:
hotter gas tends to have higher abundances (see Fig.~2).
Line 3 of Table~2 shows the regression results for all the data.
When we exclude M32 and NGC 1052 from the regression (for being hot point
sources, rather than diffuse gas), we find the abundance $A \propto T^{2.44\pm
0.69}$ (see line 4 of Table~2), with a dispersion of 0.48 dex in $A$.
This latter correlation is shown as a dotted line in Fig.~2 (where
M32 and NGC 1052 are also indicated).
We have checked our abundances against the several published $ASCA$ values
and find that in most cases they agree within the 90\% confidence
level and in all cases they agree within the 99\% confidence level.
Our PSPC temperatures are also consistent with the $ASCA$ results.
We find no correlation between our measured abundances and the stellar
abundances of the host galaxies, as inferred from $Mg_2$ indices (7S)
(see lines 5-6 of Table~2);
this is consistent with the $ASCA$ results of Loewenstein et al. (1994) and
Mushotzky et al. (1994).

				\section{
Summary and Discussion
				}

In this letter we have shown that gas temperatures are correlated with stellar
velocity dispersions, $T \propto \sigma^{1.45}$, and are everywhere hotter than
the stellar kinetic temperatures ($kT > \mu m_p \sigma^2$).
There are several possible mechanisms for the gas acquiring more specific
energy than the stars.
One possibility is that supernovae have heated the gas, but this is inconsistent
with the generally low metal abundances being found in such galaxies
(Matsushita et al. 1994; Loewenstein et al. 1994; and \S4.2 of this paper).
A more likely possibility is that the luminous parts of the galaxies are
embedded in dark matter halos which are dynamically hotter than the stars
(Fabian et al. 1986).
As stellar ejecta accumulates and flows inward in cooling flows,
movement through the gravitational potential heats the gas.
If dark matter halos are characterized by velocity dispersions that are
$\ga50\%$ greater than those of the luminous stars, then gravitational
heating can raise the gas temperatures to the observed values.

In the largest study of X-ray spectral properties of galaxies to date,
Kim, Fabbiano \& Trinchieri (1992, hereafter KFT) found that the X-ray
emission from elliptical galaxies is more complicated
than the simple thermal plasma found by Forman et al. (1985).
KFT found for galaxies with high X-ray to optical flux ratios ($L_X/L_B$)
that their spectra are indeed well described by Raymond-Smith thermal
plasma models with temperatures of $\sim1$ keV and solar abundances.
But for galaxies with intermediate values of $L_X/L_B$, KFT found spectral
hardening. At even lower values of $L_X/L_B$, the spectra require a
a very cool ``supersoft" 0.2 keV component in addition to the hard
($\sim$5 keV) component.
Although a few galaxies in our sample have unusually hard emission, we find
that the correlation of temperature with velocity dispersion persists down to
$kT\approx0.2$ keV, $L_B\approx4.7\times10^9h^{-2}$ $L_\odot$,
$\sigma\approx120$ km s$^{-1}$.
Our sample includes some of the galaxies deemed to have ``supersoft" spectra
by KFT, but we find that these sources are simply on the cool end of a
continuous distribution, with temperatures consistent with expectations
from ISM in galaxies with low velocity dispersions
(hence, shallower gravitational potentials).

While a few of the galaxies in our sample exhibit hard emission (due to X-ray
binaries, AGN, or both) we find no general hardening of the X-ray emission
as we go to lower luminosities and velocity dispersions, as might be expected
if low-luminosity ellipticals are blowing winds.
KFT find harder spectra for a subset of their sample (their ``Group 2"), which
they attribute to X-ray binaries, but our PSPC data for 4 of 8 galaxies in
this subsample (after removing NGC~1052, a Seyfert 2, which may have biased
the results of KFT) show no evidence for
having harder spectra attributable to X-ray binaries.
Thus, we detect gaseous X-ray emission without contamination by X-ray binaries
down to $L_B\approx4.7\times10^9h^{-2}$ $L_\odot$,
$L_X\approx10^{39}h^{-2}$ erg s$^{-1}$ and $\sigma\approx120$ km s$^{-1}$.
This is contrary to the predictions of some theoretical models for the
evolution of hot gas in ellipticals, in which the threshold for the onset
of supernovae-driven winds occurs at significantly higher optical luminosities
(Ciotti et al. 1991; David et al. 1991).  Such winds drive out the gas, leaving
only X-ray binaries as significant sources of X-ray emission.
M32 may be an example of such a wind-blowing galaxy, since it exhibits hard
X-ray emission and is ten times less luminous (both optically and in X-rays)
than the next less luminous galaxy in our sample.

Finally, we find that abundances are correlated with global temperatures:
$A \propto T^{2.4}$.
This may reflect an underlying trend in observed stellar abundances (more
luminous galaxies are more metal rich), but we find no detailed correlation
between our gaseous abundances and the stellar metallicity index $Mg_2$.
The $Mg_2$ index is measured within the inner few kpc of galaxies, so it may
not be representative of the galaxy out to our extraction radius of $6 r_e$.
For example, Davies, Sadler \& Peletier (1993) find that the $Mg_2$ index
drops by $\sim$25\% from the center to a radius of $\sim1 r_e$.
Analysis of $ASCA$ observations of the elliptical NGC 4636 shows there is a
gaseous abundance gradient evident from spatially resolved X-ray spectroscopy,
a gradient which persists to at least $5 r_e$ (Mushotzky et al. 1994).
Our global, emission-weighted measures of abundances obscure such structure.
The lack of correlation between the stellar and gaseous abundances
may alternatively indicate that more luminous galaxies have larger supernova
rates than less luminous galaxies, but not so large as to unbind the gas.

In a subsequent, more detailed analysis of this complete sample of ellipticals,
we will present their X-ray luminosities, their luminosity function, discuss
the Local Group dwarf galaxies, and assess additional correlations and
implications.

\acknowledgments
This research made use of the HEASARC, NED, and SkyView databases. We would
like to thank Keith Arnaud for his help with the IPC data.
This work was partially supported by the NSF and the State of Alabama through
EPSCoR grant EHR-9108761.  REW also acknowledges partial support from NASA
grants NAG 5-1718 and NAG 5-1973.


text:
\section{Introduction}

Retraction maps first appear in the literature as a topological construction in 1931~\cite{1931Borsuk}. They did not become a useful tool for designing optimization algorithms on matrix manifolds until the beginning of the XXIst century, see for instance~\cite{AbMaSeBookRetraction}. Recently, we have developed in~\cite{21MBLDMdD} a new notion of discretization map arisen from retraction maps that focuses on discretizing the configuration manifold, instead of the equations of motion. That provides a new approach to discrete mechanics that started with the foundational work~\cite{MW_Acta} based on discretizing the variational principles that define the equations of motion.

When Lagrangian or Hamiltonian systems are considered, the above-mentioned discretization map is built exploiting the properties of the tangent and cotangent bundle structures~\cite{TuHamilton,YaIs73}. In the Hamiltonian framework we obtain a systematic procedure to construct symplectic numerical methods~\cite{21MBLDMdD}, as in~\cite{hairer}.

In this paper we aim at constructing geometric integrators for optimal control problems. In 1958, Pontryagin's Maximum Principle~\cite{Pontryagin} provided necessary conditions for optimal solutions. Only a few years later, the discrete version of such conditions to obtain numerical methods for optimal control problems appeared~\cite{64Jordan}. Decades later, geometric integrators for optimal control problems have been studied and characterized~\cite{2009ChybaSymplIntOpControl,2011SinaJungeMa,2013FerKoDMdD}, even for singular optimal control problems~\cite{2009DelgadoIbort}.

Pontryagin's Maximum Principle~\cite{Pontryagin} provides the optimal control problem with Hamiltonian framework~\cite{1997Jurdjevic}.  The results from symplectic geometry~\cite{90GuiStern,LiMarle} together with the integrability algorithm~\cite{1995MMT} can be used to identify the final submanifold where the solutions live so that the discretization maps can be used. A symplectic integrator can be defined as a Lagrangian submanifold~\cite{Weinstein}. However, the geometric integrators for optimal control probelms in this paper will not preserve the symplectic 2-form, but a presymplectic 2-form~\cite{GuzMarr}.

The paper is organized as follows. After recalling Hamilton's equations in Section~\ref{Sec:HamEq}, we summarize how the discretization maps in~\cite{21MBLDMdD} are used to construct symplectic integrators in Section~\ref{Sec:DiscreteMap}. On the other hand, we introduce optimal control problems and how they can be associated with a Morse family that will be useful to run the integrability algorithm~\cite{15Morse}. As a result, using the results in the Appendix, we can construct presymplectic integrators for optimal control problems and provide an example in Section~\ref{Sec:Presympl}.

\section{Hamilton's equations}\label{Sec:HamEq}
As described in~\cite{LiMarle}, the cotangent bundle $T^*Q$ of a  differentiable manifold $Q$ is equipped with a canonical exact symplectic structure $\omega_Q=-d\theta_Q$, where $\theta_Q$ is the canonical 1-form on $T^*Q$. In canonical bundle coordinates $(q^i, p_i)$ on $T^*Q$,
$
\theta_Q= p_i\, \mathrm{d}q^i$ and
$\omega_Q= \mathrm{d}q^i\wedge \mathrm{d}p_i\; .
$
The Hamiltonian vector field associated with a Hamiltonian function $H: T^*Q\rightarrow {\mathbb R}$ must satisfy:
$
\imath_{X_H}\omega_Q=\mathrm{d}H\; ,
$
whose integral curves are solution to Hamilton's equations:
\begin{equation*}
	\frac{\mathrm{d}q^i}{\mathrm{d}t}=\frac{\partial H}{\partial p_i}\; ,\quad
	\frac{\mathrm{d}p_i}{\mathrm{d}t}=-\frac{\partial H}{\partial q^i}\; .
\end{equation*}
Some fundamental properties of the Hamiltonian dynamics are:
\begin{itemize}
	\item  Preservation of energy, that is, the Hamiltonian function is preserved:  $$0=\omega_Q(X_H, X_H)=dH(X_H)=X_H(H)\,.$$

	\item Preservation of the symplectic form, that is, the Lie derivative of the 2-form $\omega_Q$ with respect to the Hamiltonian vector field vanishes: $L_{X_H}\omega_Q=0$. Equivalently, if $\{\phi^t_{X_H}\}$ is the flow of $X_H$, then
	\[
	(\phi^t_{X_H})^*\omega_Q=\omega_Q\; .
	\]
\end{itemize}
Symplectic integrators~\cite{hairer} were designed to preserve the configuration manifold and preserve the canonical symplectic form.

\section{From retraction maps to discretization maps}  \label{Sec:DiscreteMap}
As described in~\cite{AbMaSeBookRetraction}, a retraction map on a manifold $M$ is a smooth map
$R\colon U\subseteq TM \rightarrow M$ where $U$ is an open subset containing the zero section  of the tangent bundle
such that the restriction map $R_x=R_{|T_xM}\colon T_xM\rightarrow M$ satisfies
\begin{enumerate}
	\item $R_x(0_x)=x$ for all $x\in M$,

	\item
	${\rm D}R_x(0_x)=T_{0_x}R_x={\rm Id}_{T_xM}$ where we identify $T_{0_x}T_xM\simeq T_xM$.
\end{enumerate}

\begin{example}
If  $(M, g)$ is a Riemannian manifold, then the exponential map $\hbox{exp}^g: U\subset TM\rightarrow M$ is a typical example of retraction map:
$
\hbox{exp}^g_x(v_x)=\gamma_{v_x}(1),
$
where $\gamma_{v_x}$ is the unique  Riemannian geodesic satisfying $\gamma_{v_x}(0)=x$ and $\gamma'_{v_x}(0)=v_x$~\cite{doCarmo}.
\end{example}

In \cite{21MBLDMdD}  retraction maps have been used to define discretization maps $R_d\colon U \subset TM  \rightarrow M\times M$, where  $U$ is an open neighbourhood of the zero section of $TM$,
\begin{eqnarray*}
	R_d\colon U \subset TM & \longrightarrow & M\times M\\
	v_x & \longmapsto & (R^1(v_x),R^2(v_x))\, .
\end{eqnarray*}
Discretization maps satisfy the following properties:
\begin{enumerate}
	\item $R_d(0_x)=(x,x)$, for all $x\in M$.
	\item $T_{0_x}R^2_x-T_{0_x}R^1_x={\rm Id}_{T_xM} \colon T_{0_x}T_xM\simeq T_xM \rightarrow T_xM$ is equal to the identity map on $T_xM$ for any $x$ in $M$.
\end{enumerate}
Thus, the discretization map $R_d$ is a local diffeomorphism.

\begin{example}
Examples of discretization maps on Euclidean vector spaces are:
\begin{itemize}
	\item Explicit Euler method:  $R_d(x,v)=(x,x+v).$
	\item Midpoint rule:  $R_d(x,v)=\left( x-\dfrac{v}{2}, x+\dfrac{v}{2}\right).$
		\item $\theta$-methods with $\theta\in [0,1]$:  \hspace{3mm} $R_d(x,v)=\left( x-\theta \, v, x+ (1-\theta)\, v\right).$
\end{itemize}
\end{example}

\subsection{Cotangent lift of discretization maps}

We want to define a discretization map on $T^*Q$, that is,
$R^{T^*}_d: TT^*Q \rightarrow T^*Q\times T^*Q$. The domain lives where the Hamiltonian vector field takes value. Such a map will be obtained by cotangently lifting a discretization map $R_d\colon TQ\rightarrow Q\times Q$ so that the construction $R^{T^*}_d$ will be a symplectomorphism. In order to do that, we need the following three symplectomorphisms (see \cite{21MBLDMdD} for more details):
\begin{itemize}
\item The cotangent lift of a diffeomorphism $F: M_1\rightarrow M_2$ defined by:
	\begin{equation*}
	\hat{F}: T^*M_1 \longrightarrow  T^*M_2 \mbox{ such that }
\hat{F}=(TF^{-1})^*.
	\end{equation*}
	\item The canonical symplectomorphism:
	\begin{equation*} \alpha_Q\colon T^*TQ  \longrightarrow  TT^*Q  \mbox{ such that }  \alpha_Q(q,v,p_q,p_v)= (q,p_v,v,p_q).
	\end{equation*}

	\item  The symplectomorphism between $(T^*(Q\times Q), \omega_{Q\times Q})$     and
	$(T^*Q\times T^*Q, \Omega_{12}=pr_2^*\omega_Q-pr^*_1\omega_Q)$:
		\begin{equation*}
	\Phi:T^*Q\times T^*Q \longrightarrow T^*(Q\times Q)\; , \;
	\Phi(q_0, p_0; q_1, p_1)=(q_0, q_1, -p_0, p_1).	\end{equation*}
	\end{itemize}
The following diagram summarizes the construction procress from $R_d$ to $R_d^{T^*}$:
	\begin{equation*}
\xymatrix{ {{TT^*Q }} \ar[rr]^{{{R_d^{T^*}}}}\ar[d]_{\alpha_{Q}} && {{T^*Q\times T^*Q }}  \\ T^*TQ \ar[d]^{\pi_{TQ}}\ar[rr]^{	\widehat{R_d}}&& T^*(Q\times Q)\ar[u]_{\Phi^{-1}}\ar[d]^{\pi_{Q\times Q}}\\ TQ \ar[rr]^{R_d} && Q\times Q }
\end{equation*}

\begin{proposition}\cite{21MBLDMdD}
	Let $R_d\colon TQ\rightarrow Q\times Q$ be a discretization map on $Q$. Then $${{R_d^{T^*}=\Phi^{-1}\circ \widehat{R_d}\circ \alpha_Q\colon TT^*Q\rightarrow T^*Q\times T^*Q}}$$
is a discretization map  on $T^*Q$.
\end{proposition}
\begin{corollary}\cite{21MBLDMdD} The discretization map
 ${{R_d^{T^*}}}=\Phi^{-1}\circ 	(TR_d^{-1})^* \circ \alpha_Q\colon T(T^*Q)\rightarrow T^*Q\times T^*Q$ is a symplectomorphism between $(T(T^*Q), {\rm d}_T \omega_Q)$ and $(T^*Q\times T^*Q, \Omega_{12})$.
\end{corollary}
\begin{example}\label{example3} On $Q={\mathbb R}^n$ the discretization map
	$R_d(q,v)=\left(q-\frac{1}{2}v, q+\frac{1}{2}v\right)$ is cotangently lifted to
		$$R_d^{T^*}(q,p,\dot{q},\dot{p})=\left( q-\dfrac{1}{2}\,\dot{q}, p-\dfrac{\dot{p}}{2}; \; q+\dfrac{1}{2}\, \dot{q}, p+\dfrac{\dot{p}}{2}\right)\, .$$
\end{example}

\subsection{Symplectic methods for Hamilton's equations}

For Hamilton's equation we automatically produce a symplectic integrator using
the discretization map $R_d^{T^*}\colon TT^*Q \rightarrow T^*Q\times T^*Q$ which is the cotangent lift of a discretization map on $Q$.

\begin{proposition}\cite{21MBLDMdD}
The numerical method defined by
	$$h\, X_H \left(\tau_{TQ}\left( \left(R^{TT^*Q}_d\right)^{-1}(q_k, p_k; q_{k+1}, p_{k+1})\right)\right)=\left(R^{TT^*Q}_d\right)^{-1}(q_k, p_k; q_{k+1}, p_{k+1})\,, $$
	is a symplectic integrator for the Hamiltonian system given by $H\colon T^*Q\rightarrow \mathbb{R}$.
	\end{proposition}
\section{Optimal control problems and Morse families}\label{Sec:OptMorse}
An optimal control problem (OCP) is given by a vector field depending on parameters called controls, a cost function and some end-point conditions. A solution of an OCP must be an integral curve of the vector field for specific controls, $\dot{q}=X(q,u)$, so that the functional $\int^{t_f}_{t_0} F(q(t),u(t))\, {\rm d }t,$ is minimized and the end-point conditions satisfied.

Typically, OCP are solved using Pontryagin's Maximum Principle~\cite{1997Jurdjevic,Pontryagin} that provides the problem with a Hamiltonian framework.
Let $U$ be the set of admissible controls, the associated Pontryagin's Hamiltonian function is:
\begin{equation*}
H \colon T^*Q\times U \rightarrow \mathbb{R}\, ,\quad
H(q,p,u)=\langle p, X(q,u) \rangle - F(q,u)\, ,\end{equation*}
where $\langle\cdot, \cdot \rangle$ denotes the natural pairing between $T_q^*Q$ and $T_qQ$.
A Morse family is another geometric object that can be used to define Lagrangian submanifolds. Such a notion was first  introduced by L. H\"ormander~\cite{Hormander}.  It is proved in~\cite{15Morse} that the Pontryagin's Hamiltonian function could be a Morse family over the projection ${\rm pr}_1\colon T^*Q\times U \rightarrow T^*Q$ onto the first factor if the image of the differential of $H$ and the conormal bundle
\begin{equation*}(\ker \, {\rm T}\pi)^0=\left\{\alpha\in T_\mu ^*(T^*Q\times U)\; |\; \langle \alpha, v\rangle=0, \hbox{ for all } v\in \ker T_{\mu}\pi\right\}\subset T^*(T^*Q\times U)
\end{equation*}
are transverse in $T^*(T^*Q\times U)$, that is,
\begin{equation*} {\rm T}_\alpha( {\rm d}H(T^*Q\times U))+T_\alpha (\ker \, {\rm T}\pi)^0=T_\alpha (T^*(T^*Q\times U)),
\end{equation*}
for all $\alpha\in (\ker \, {\rm T}\pi)^0 \cap  {\rm d}H(T^*Q\times U)  \subseteq T^*(T^*Q\times U)$.

\begin{proposition}\cite{15Morse} Pontryagin's Hamiltonian   $H\colon T^*Q\times U  \rightarrow \mathbb{R}$ defines a \textbf{Morse family over} the projection ${\rm pr}_1\colon T^*Q\times U \rightarrow T^*Q$ onto the first factor if and only if  the matrix
	\begin{equation*}
	{\rm D}_{(q,p,u)} \begin{pmatrix} \dfrac{\partial H}{\partial u}\end{pmatrix}=\begin{pmatrix}\dfrac{\partial^2 H}{\partial q^i \partial u^a} & \dfrac{\partial^2 H}{\partial p_i \partial u^a} & \dfrac{\partial^2 H}{\partial u^a \partial u^b}
	\end{pmatrix}_{(q,p,u)}
	\end{equation*}
	has maximum rank for all $(q,p,u)\in T^*Q\times U$ such that $\dfrac{\partial H}{\partial u}(q,p,u)=0$.
\end{proposition}

When the controls are in the interior of the set $U$, the necessary conditions of Pontryagin's Maximum Principle can be rewritten as the following Lagrangian submanifold ${\mathcal L}_H$ of $(T^*T^*Q, \omega_{T^*Q})$:
\[
{\mathcal L}_H= \left\{(q, p, P_q, P_p)\;\left| \exists \; u\in U \; \mbox{s. t. } \begin{array}{l}    P_q=\frac{\partial H}{\partial q}(q,p,u), P_p=\frac{\partial H}{\partial p}(q,p,u), \\  \frac{\partial H}{\partial u}(q,p,u)=0\,  \end{array}\right.\right\}\, .
\]

The OCP is regular if $(\pi_{T^*Q})_{|{\mathcal L}_H}: {\mathcal L}_H\rightarrow T^*Q$ is a local diffeomorphism, otherwise it is called singular.
Observe that in general ${\mathcal L}_H$ is not horizontal, that is, it is not transverse to the fibers of the canonical cotangent projection $\pi_{T^*Q}$. Consequently, it is not the image of the differential of a function on $T^*Q$ \cite{15Morse}.

Due to the symplectomorphism between $(TT^*Q, d_T\omega_Q)$ and $(T^*T^*Q, \omega_{T^*Q})$ described in~\cite{TuHamilton}, the dynamics of an optimal control problem can also be given as the following Lagrangian submanifold in $(TT^*Q, d_T\omega_Q)$:
\begin{equation*}
S_{H}=\sharp_{\omega_{T^*Q}}({\mathcal L}_H)=\{v\in TT^*Q \, | \; i_v\omega_{T^*Q}\in {\mathcal L}_H\}
\end{equation*}
Thus, a solution of the OCP  is a curve $\sigma$ in $T^*Q$ such that $\dot{\sigma}(t)$ lies in $S_{H}$.

In general, the solutions of the OCP are consistently defined in a submanifold of $T^*Q$ contained in $\tau_{T^*Q}(S_H)\subseteq T^*Q$. Thus, the integrability algorithm~\cite{1995MMT} can be used to obtain the integrable part of $S_{H}$ in $T^*Q$. First, we define
$S^0_H=S_{H}$ and $P_H^0=\tau_{T^*Q}(S_H)$.

The following steps of the algorithm are defined by
\begin{equation*}
S^k_H=T P_H^{k-1}\cap S_H^{k-1}, \quad P_H^k=\tau_{T^*Q}(S_H^k).
\end{equation*}
If the algorithm stabilizes at step $k_f$ of the constraint algorithm, there exists a final submanifold (possibly empty or singular) satisfying $S^{k_f}_H=T P^{k_f}_H \cap S^{k_f}_H$, that will be denoted by $S^f_{H}$.
On the base manifold $T^*Q$ and the tangent bundle $TT^*Q$ the algorithm generates the following two sequences of submanifolds in $T^*Q$ and $TT^*Q$, respectively:
\begin{equation*}
\begin{array}{c}
P^f_H \stackrel{i_{k_f}}{\hooklongrightarrow} P^{k_f-1}_H \stackrel{i_{k_f-1}}{\hooklongrightarrow}\cdots  P^1_H \stackrel{i_{1}}{\hooklongrightarrow} P^0_H \stackrel{i_{0}}{\hooklongrightarrow} T^*Q\, ,
\\
S^f_H \stackrel{j_{k_f}}{\hooklongrightarrow} S^{k_f-1}_H \stackrel{j_{k_f-1}}{\hooklongrightarrow}\cdots S^1_H \stackrel{j_{1}}{\hooklongrightarrow} S^0_H \stackrel{j_{0}}{\hooklongrightarrow} TT^*Q\, .
\end{array}
\end{equation*}

As a consequence, for every $\alpha$ in $P^f_H$ there exists $V$ in $T_{\alpha}S^f_H\subset TT^*Q$.
Hence, the original dynamical system has solution in the submanifold $P^f_H$.
Denote by $i_f: P^f_H\hookrightarrow T^*Q$ the canonical inclusion and by $\omega_f=i_f^*\omega_{T^*Q}$ the pullback of the canonical symplectic 2-form on $T^*Q$. Note that $\omega_f$ is now a presymplectic 2-form (see Appendix).

In conclusion, a solution to the $OCP$  is a curve $\sigma$ on $P^f_H$ such that there exist controls $u$ satisfying
\begin{equation*}\label{presymplectic-equation}
i_{\dot{\sigma}(t)}\omega_f (\sigma(t))=dH_{f}^{u}(\dot{\sigma}(t))\, ,
\end{equation*}
 where $H_f^u: S^f_H\rightarrow {\mathbb R}$ is given by $H_f^u(\dot{\sigma}(t))=H(j_f(\dot{\sigma}(t)),u(t))$. Therefore, the dynamics that we need to preserve with our numerical methods is presymplectic instead of the most classical symplectic preservation property.

 \begin{proposition}\label{proppre} The submanifold
 $S^f_H$ is a Lagrangian submanifold of the presymplectic manifold  $(TP^f_H, d_T\omega_f)$.
 \end{proposition}
 \begin{proof}
 The results follows because $\omega_f=i_f^*\omega_{T^*Q}$, ${S}_H$ is a Lagrangian submanifold of $(TT^*Q, d_T\omega_Q)$ and
 $
 d_T\omega_f=d_Ti_f^*\omega_{T^*Q}=(Ti_f)^*d_T\omega_{T^*Q}
 $.
 \end{proof}

\section{Presymplectic integrators for optimal control problems}\label{Sec:Presympl}
In this section we will use the cotangent lift of a discretization map to define a presymplectic integrator for optimal control problems, once we have run the integrability algorithm and know the final submanifold $P^f_H$ of $T^*Q$.

We restrict the cotangent lift of a discretization map $R^{T^*}_d\colon T(T^*Q)\rightarrow T^*Q\times T^*Q$ to the submanifold $TP^f_H$ and define the submanifold $P^f_{H, d}$  of $T^*Q\times T^*Q$ by
\[
P^f_{H, d}\colon =R_d^{T^*}(TP^f_H)\,.
\]
Introducing the inclusion $j_f^d:P^f_{H, d}\hookrightarrow T^*Q\times T^*Q$ the following diagram summarizes the construction process:
\begin{equation*}
\xymatrix{ {{T^*Q\times T^*Q }} && {TT^*Q }\ar[ll]_{{{R_d^{T^*}}}} \ar[rd]^{\tau_{T^*Q}}&&
	T^*T^*Q\ar[ll]_{\sharp_{\omega_Q}}\ar[ld]^{\pi_{T^*Q}}&\ar@{_(->}[l] {\mathcal  L}_H
	\\
	P^f_{H, d}\ar[u]_{j^d_f}&&TP^f_H\ar[u]_{j_f}\ar[rd]_{\tau_{P^f_H}} \ar[ll]_{{{R_d^{T^*}|_{TP^f_H}}}}&T^*Q&\\
	&&&P^f_H\ar[u]_{i_f} &
}
\end{equation*}

\begin{definition}\label{method}
We define the OCP geometric integrator as
\[
\left\{\begin{array}{l}
\frac{1}{h}\left(R_d^{T^*}\right)^{-1}(q_k, p_k; q_{k+1}, p_{k+1})\in  \left(S^f_H\right)_{\tau_{T^*Q}\left(\left(R_d^{T^*}\right)^{-1}(q_k, p_k; q_{k+1}, p_{k+1})\right)}\, ,\\
 (q_k, p_k; q_{k+1}, p_{k+1})\in P^f_{H, d}
\end{array}\right.
\]
\end{definition}
Propositions \ref{proppre}  and \ref{propopre} guarantees the presymplecticity of the method.
\begin{proposition}
The OCP geometric integrator in Definition~\ref{method} preserves the presymplectic 2-form
$\Omega_{f, d}=(j_f^d)^*\Omega_{12}$.
\end{proposition}

\subsection{Example}
As an academic example we consider the singular optimal control problem on $\mathbb{R}^2$ given by the control equations
$
\dot{x}=f(x)+u_1, \; \dot{y}=y\,
$,
where $u_1$, $u_2 \in \mathbb{R}$,
and the cost functional
\[
\int_{t_0}^{t_f} \left(
\frac{1}{2}x^2+\frac{1}{2}y^2+xu_1+yu_2+\frac{1}{2}u_1^2\right)\, dt\; .
\]
Then Pontryagin's Hamiltonian is
\[
H(x,y,p_x,p_y,u_1,u_2)=p_x (f(x)+u_1)+p_y y-\frac{1}{2}x^2-\frac{1}{2}y^2-xu_1-yu_2-\frac{1}{2}u_1^2\, .
\]
The Lagrangian submanifold ${\mathcal L}_H$ of $(T^*(T^*\mathbb{R}^2),\omega_{T^*\mathbb{R}^2})$ is
\begin{eqnarray*}
{\mathcal L}_H&=&
\{(x,y, p_x, p_y, P_{x}, P_y, P_{p_x}, P_{p_y})\ |\;
P_{x}=p_x\, f'(x)-x-u_1, P_y=p_y-y-u_2,\\
&&P_{p_x}=f(x)+u_1, P_{p_y}=y, \; p_x-x-u_1=0, \; y=0\}\, .
\end{eqnarray*}
Therefore, we obtain the following Lagrangian submanifold of $(TT^*\mathbb{R}^2,{\rm d}_T\omega_{\mathbb{R}^2})$:
\begin{eqnarray*}
{S}_H&=&
\{(x,y, p_x, p_y, \dot{x}, \dot{y}, \dot{p}_x, \dot{p}_y)\ |\;
\dot{x}=f(x)+u_1, \dot{y}=y,\dot{p}_x=-p_x\, f'(x)+x+u_1,\\
&& \dot{p}_y=-p_y+y+u_2, \; p_x-x-u_1=0, \; y=0\}\,.
\end{eqnarray*}
Applying the constraint algorithm we immediately deduce that
\[
P^f_H=\{(x,0, p_x, p_y)\in {\mathbb R}^4\}\equiv {\mathbb R}^3
\]
which is a presymplectic manifold with the 2-form $\omega_f=dx\wedge dp_x$
whose $\ker \omega_f=\hbox{span}\{ {\partial}/{\partial p_y}\}$. The corresponding Lagrangian submanifold of the presymplectic manifold
$(TP^f_H, d_T\omega_f)$
is:
\begin{eqnarray*}
S_H^f&=&\{(x, p_x, p_y, \dot{x}, \dot{p}_x, \dot{p}_y)\ |\;
\dot{x}=f(x)+u_1, \dot{p}_x=-p_x f'(x)+x+u_1,\\
&& \dot{p}_y=-p_y+u_2, \; p_x-x-u_1=0 \}\, .
\end{eqnarray*}
Applying now the discretization map derived in Example \ref{example3} we obtain that
\[
P^f_{H,d}=\{(x_k, 0, (p_x)_k, (p_y)_k, x_{k+1}, 0, (p_x)_{k+1}, (p_y)_{k+1})\in {\mathbb R}^8\}\equiv {\mathbb R}^3\times {\mathbb R}^3
\]
and the resulting presymplectic integrator is described by the equations:
\begin{eqnarray*}
\frac{x_{k+1}-x_k}{h}&=&f\left(\frac{x_k+x_{k+1}}{2}\right)+(u_1)_k,\\ \frac{(p_x)_{k+1}-(p_x)_k}{h}&=&-\frac{(p_x)_{k}+(p_x)_{k+1}}{2} f'\left(\frac{x_k+x_{k+1}}{2}\right)+
\frac{x_k+x_{k+1}}{2}+(u_1)_k,\\
 \frac{(p_y)_{k+1}-(p_y)_k}{h}&=&-\frac{(p_y)_{k}+(p_y)_{k+1}}{2}+(u_2)_k, \\
  0&=&\frac{(p_x)_{k}+(p_x)_{k+1}}{2}-\frac{x_k+x_{k+1}}{2}-(u_1)_k\, .
\end{eqnarray*}
The discrete equations for $p_y$ and $u_2$ are decoupled from the rest that can be written in a more compact way for $(x_k,(p_x)_k,x_{k+1},(p_x)_{k+1})$:
\begin{eqnarray*}
\frac{x_{k+1}-x_k}{h}&=&f\left(\frac{x_k+x_{k+1}}{2}\right)+\frac{(p_x)_{k}+(p_x)_{k+1}}{2}-\frac{x_k+x_{k+1}}{2},\\ \frac{(p_x)_{k+1}-(p_x)_k}{h}&=&-\frac{(p_x)_{k}+(p_x)_{k+1}}{2}f'\left(\frac{x_k+x_{k+1}}{2}\right)+
\frac{x_k+x_{k+1}}{2}\\
&&+\frac{(p_x)_{k}+(p_x)_{k+1}}{2}-\frac{x_k+x_{k+1}}{2}\,.
\end{eqnarray*}

\section{Conclusions and future
work}

In this paper we have only studied normal solutions, both regular and singular, for optimal control problems. The same technique can be applied for abnormal solutions where Pontryagin's Hamiltonian does not depend on the cost function. Moreover, we plan to construct geometric integrators for Dirac systems~\cite{Melvin}.

\section*{Appendix: Presymplectic geometry and Lagrangian submanifolds}

As introduced in~\cite{GuzMarr}, a preymplectic structure on a finite dimensional manifold $M$ is a closed 2-form $\omega$ on $M$. We say that $(M, \omega)$ is a presymplectic manifold. The kernel of the presymplectic structure at a point $x$ in $M$ is a vector subspace of the tangent space of $M$ at $x$ that it is  not necessarily zero as in the symplectic case. Remember that
$
\ker \omega_x=\{v\in T_x M\; |\; i_v \omega_x=0\}
$.
If $\dim \ker\omega_x=0$ for all $x\in M$, the presymplectic structure is non degenerate. Hence, it is a symplectic structure~\cite{90GuiStern,LiMarle}.

\begin{definition}
A submanifold ${\mathcal L}$ of dimension $r$ of the presymplectic manifold $(M, \omega)$ with canonical inclusion $i: {\mathcal L}\hookrightarrow M$ is said to be Lagrangian if the pullback of $\omega$ by the inclusion vanishes, that is, $i^*\omega=0$, and
\begin{equation}\label{eq:dim_preL}
r=\frac{ \dim M- \dim \left(\ker \omega_x\right)}{2} +\dim (T_x{\mathcal L}\cap \ker \omega_x), \quad \hbox{for all} \;  x\in M\,.
\end{equation}
\end{definition}
When $\omega$ is a symplectic structure, Equation~\eqref{eq:dim_preL} implies that the dimension of ${\mathcal L}$ is half of the dimension of $M$ and we recover the classical definition of Lagrangian submanifold in symplectic geometry~\cite{LiMarle,Weinstein}.

A smooth map $f: M\rightarrow N$ between two presymplectic manifolds $(M, \omega_M)$ and $(N, \omega_N)$ is a presymplectic map if $f$ preserves the presymplectic structures, that is,  $f^*\omega_N=\omega_M$. From that notion, it is possible to construct the following Lagrangian submanifolds.  The proofs come from the above definitions.
\begin{proposition} \label{Prop:preGraphf} Let $(M, \omega_M)$ and $(N, \omega_N)$ be presymplectic manifolds.
If $f: M\rightarrow N$  is a presymplectic diffeomorphism, then
\[
\hbox{\rm Graph} f=\{(x, f(x)) \, | \, x\in M\}
\]
is a $m$-dimensional  Lagrangian submanifold of $(M\times N, \Omega_{M\times N}=\omega_M-\omega_N)$, where $m$ is the dimension of $M$.
\end{proposition}

\begin{proposition}\label{propopre}
Let $(N, \omega_N)$ be a presymplectic manifold and $f: M\rightarrow N$ a diffeomorphism. Then $(M,\omega_M)$ is a presymplectic manifold with the presymplectic structure $\omega_M= f^*\omega_N$. Moreover, if ${\mathcal L}_N$  is a Lagrangian submanifold of $(N, \omega_N)$, then ${\mathcal L}_M=f^{-1}({\mathcal L}_N)$ is a Lagrangian submanifold of $(M,\omega_M)$.
\end{proposition}