# Distributionally Safe Reinforcement Learning under Model Uncertainty: A
Single-Level Approach by Differentiable Convex Programming

Alaa Eddine Chriat1 and Chuangchuang Sun1 1The authors are with the Aerospace
Engineering Department, Mississippi State University, Starkville, MS 39759,
USA. Emails<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Safety assurance is uncompromisable for safety-critical environments with the
presence of drastic model uncertainties (e.g., distributional shift),
especially with humans in the loop. However, incorporating uncertainty in safe
learning will naturally lead to a bi-level problem, where at the lower level
the (worst-case) safety constraint is evaluated within the uncertainty
ambiguity set. In this paper, we present a tractable distributionally safe
reinforcement learning framework to enforce safety under a distributional
shift measured by a Wasserstein metric. To improve the tractability, we first
use duality theory to transform the lower-level optimization from infinite-
dimensional probability space where distributional shift is measured, to a
finite-dimensional parametric space. Moreover, by differentiable convex
programming, the bi-level safe learning problem is further reduced to a
single-level one with two sequential computationally efficient modules: a
convex quadratic program to guarantee safety followed by a projected gradient
ascent to simultaneously find the worst-case uncertainty. This end-to-end
differentiable framework with safety constraints, to the best of our
knowledge, is the first tractable single-level solution to address
distributional safety. We test our approach on first and second-order systems
with varying complexities and compare our results with the uncertainty-
agnostic policies, where our approach demonstrates a significant improvement
on safety guarantees.

##  I INTRODUCTION

In many real-world applications, there often can be unmodelled dynamics and
uncertainties, both internally (e.g., systems failure/ dysfunctionality,
perceptional noise) and externally (e.g., gust, tough terrain). Those
uncertainties make control problems and training processes more challenging
and susceptible to errors, potentially leading to undesirable/ disastrous
outcomes. Moreover, there are often more restrictive settings where
environmental configurations change more drastically. For example, there can
always be “unknown unknowns” during the deployment of autonomous systems, such
as off-road autonomy, advanced air mobility, etc. Also, one of the fundamental
challenges towards (super) human-level intelligence for robots is the ability
of trained policies to generalize beyond the specific environments they
initially ever encountered during training. However, with the presence of
those pervasive uncertainties, safety is not compromisable in many safety-
critical scenarios with human health/ lives at stake. For example, in self-
driving cars, decision-making in real-time in complex and uncertain
environments can be very challenging, considering various surroundings, moving
pedestrians, etc.

Learning-based control has attracted lots of attention to combine the
advantages of machine learning and modern control theory. On one hand, data-
driven machine learning methods have arisen to learn control policies from
interactions with the environments without an accurate prior model. On the
other hand, while control theory has been extensively investigated and in
general has rigorous formal guarantees of performance, it often needs
relatively accurate models, which are expensive, if not impossible to obtain.
To achieve robust and safe decision-making and control with the presence of
uncertainties, there has been a large volume of learning-based approaches in
the literature. In general, robust learning and control treats uncertainty
with a soft criterion and aims to get the performant policy under the (worst-
case) uncertainty, yielding a minimax optimization problem [1, 2, 3, 4, 5, 6,
7]. While they can provide some degrees of safety, it is not guaranteed. On
the other hand, in safe learning, safety is imposed as a hard constraint [8,
9, 10, 11, 12, 13, 14, 15]. However, incorporating uncertainty in safe
learning will naturally lead to a bi-level problem, where at the lower level
the (worst-case) safety constraint is evaluated within the uncertainty
ambiguity set.

In this paper, we consider a distributinoally safe reinforcement learning
(DSRL) problem with a distributional shift measured by a Wasserstein metric.
Such distributional shift quantifies drastic uncertainties and it leads to a
distributional optimization at the lower level of the overall safe learning
problem to evaluate the worst-case safety constraint. Both the distributional
optimization and the bi-level nature make the distributional safe learning
problem intractable. We first use duality theory to transform the lower-level
optimization problem into the finite-dimensional parametric space instead of
infinite-dimensional probability space. Moreover, we further reduce the bi-
level safe learning problem as single-level by differentiable convex
programming [16] to sequentially find the worst case uncertainty and
simultaneously guarantee that the safety constraints are satisfied. Both
modules are computationally efficient, with the former as a projected gradient
ascent and the latter as a convex quadratic programming. Moreover, this
constrained learning pipeline is end-to-end differentiable. To the best of our
knowledge, this is the first attempt to address distributional safety with a
tractable single-level solution.

The rest of the paper is organized as follows. Section II reviews the
preliminaries needed to complete this work. In section III-A, we formulate our
distributionally safe approach using distributional robust optimization. In
section III-B we reduce the problem into a single-level learning-based
approach. Section IV contains numerical simulations and results comparisons,
and section V has concluding notes of our work.

###  I-A Related works on safe/robust learning

There are multiple ways to mitigate quantified uncertainty in the context of
safe/ robust learning, such as [17, 18], those focusing on the area of
robotics [19, 20], and comprehensive surveys [21, 22]. Specifically, the
uncertainty variable can be treated as a context variable representing
different tasks and can be subsequently solved as multi-task or meta-learning
problems [23, 24]. Moreover, given optimization theories, robust learning
algorithms have also been developed based on interior point methods [8, 9],
successive convexification [10] and (augmented) Lagrangian methods [11, 12,
13, 14, 15]. In learning-based control, Lyapunov theory, model predictive
control, and control barrier functions are also employed to develop robust
learning algorithms [25, 26, 27, 28, 29, 30]. Additionally, if the uncertainty
is considered in the worst-case scenario, minimax policy optimization [1, 2,
3] or its generalization Stackelberg games [4, 5, 6, 7] are often the
frameworks to promote resilience. Other works include meta-adaptive nonlinear
control integrating learning modules for fast adaptation in unpredictive
settings [31, 32].

In terms of robust learning and control under distributional shift, model-
based approaches [33] such as approximate dynamic programming [34, 35] and
model predictive control [36, 37, 38] have been proposed, either under a
Wasserstein metric or chance-constrained criterion. In the model-free regime,
one line of work is to generate environments/ tasks with distributional shifts
for policy training to achieve robustness [39, 40, 41]. Moreover, to balance
the worst-case (robustness) and average performance, [42] trains policies over
task groups by adding regularization to the worst possible outcomes. In
offline RL [43], [44] proposes a distributionally robust formulation with
tabular Markov decision processes with an uncertainty set specified by the
Kullback-Leibler (KL) divergence. Those approaches are often based on duality
in the distributionally robust optimization theories [45, 46, 47].

##  II PRELIMINARIES

###  II-A High-order CBF

In control theory, control barrier functions play a crucial role in ensuring
that a dynamic system can accomplish target objectives while ensuring that
safety is not compromised, a control barrier function essentially evaluates
the system’s safety and returns a scalar quantity. Consequently, our objective
is to determine a control input that keeps the system within the safety
boundaries as determined by the control barrier function. Mathematically,
consider the nonlinear control-affine system:

$\dot{x}(t)=f(x(t))+g(x(t))u(t)$ (1)

where $f$ and $g$ are globally Lipschitz, $x\in\mathbb{R}^{n}$ and
$u\in\mathbb{R}^{m}$ are the states and control inputs, respectively,
constrained in closed sets, with initial condition $x(t_{0})=x_{0}$. The time
dependency on $t$ can be omitted for notational simplicity.

###### Definition 1

[28]$h:\mathbb{R}^{n}\rightarrow\mathbb{R}$ is a barrier function for the set
$C=\left\\{x\in\mathbb{R}^{n}:h(x)\geqslant 0\right\\}$ if $\exists$ an
extended class-$\mathcal{K}$ function $\alpha(\bullet)$ such that:

$\begin{gathered}\sup_{u\in U}[L_{f}h(x)+L_{g}h(x)u+\alpha(h(x))]\geqslant
0\\\ \inf_{\text{int}(C)}[\alpha(h(x))]\geqslant 0\text{ \quad and
\quad}\lim_{\partial C}\alpha(h(x))=0\end{gathered}$ (2)

Because not all systems are first-order in inputs, we can use higher-order
control barrier functions to constrain higher-order systems.

###### Definition 2

[48]For the nonlinear system (1) with the $m^{th}$ differentiable function
$h(x)$ as a constraint, we define a sequence of functions $\psi_{i}$ with
$i\in\\{1,2,...,m\\}$, starting from $\psi_{0}=h(x)$:
$\psi_{i}(x,t)=\dot{\psi}_{i-1}(x,t)+\alpha_{i}\left(\psi_{i-1}(x,t)\right)$.
and define $C_{i}(t)$ sequence of safe sets associated with each $\psi_{i}$:
$C_{i}(t)=\left\\{x\in\mathbb{R}^{n}:\psi_{i-1}(x,t)\geqslant 0\right\\}$. The
function $h(x)$ is a high order control barrier function if there exist
extended class-$\mathcal{K}$ functions $\alpha_{i}(\bullet)$ such that
$\psi_{m}(x,t)\geqslant 0$.

Control barrier functions offer significant promise in the development of
secure dynamic systems, with applications in various robotic and autonomous
systems.

###  II-B Reinforcement learning: DDPG

Reinforcement Learning agents learn to make sequential decisions by
interacting with an environment to achieve a certain goal. One popular RL
approach is Deep Deterministic Policy Gradient (DDPG), an off-policy algorithm
that addresses continuous action space systems. In DDPG, the agent employs two
neural network architectures, the actor and the critic. The actor-network
represents the agent’s policy, mapping the observed state directly to a
specific action in the continuous action space
$\mu\left(s\mid\theta^{\mu}\right)$ and aiming to maximize the cumulative
reward function as $J(\theta)=\sum_{k=1}^{T}\gamma^{k}{R}(s_{k},a_{k})$, where
$\gamma$ and $R$ are the discount factor and rewards, respectively. The
critic-network, on the other hand, evaluates the quality of the actions chosen
by the actor $Q\left(s,a\mid\theta^{Q}\right)$. It estimates the expected
cumulative reward, known as the Q-value, by taking both the current state and
the action as input
$\mathcal{L}(\theta)=\mathbb{E}_{s,a,r,s^{\prime}}\left((y-Q(s,a|\theta^{Q})\right)^{2}$
and $y=R+\gamma
Q^{\prime}\left(s,\mu^{\prime}(s\mid\theta^{\mu^{\prime}})\mid\theta^{Q^{\prime}}\right)$.
These networks work together to guide the agent into achieving the desired
goal. DDPG employs a replay buffer, which stores a tuple of past experiences
$\langle S,A,\mathcal{R},S^{\prime}\rangle$, where $S$ is a set of agent
states in the environment, $A$ is a set of agent actions, $\mathcal{R}$ is the
reward function, and $S^{\prime}$ is a set of agent next states. During
training, the agent samples random batches of experiences from this buffer to
update the actor and critic networks. DDPG’s efficiency allows training agents
to perform accurate and complex tasks in a wide range of real-world scenarios.

###  II-C Differentiable convex programming

Differentiable convex programming is a powerful technique that allows
computation of the gradients of an optimization problem objective function
with respect to the parameters of the problem, by taking matrix
differentiation of the Karush-Kuhn-Tucker (KKT) conditions. A notable example
of a differentiable optimization method is CVXlayers [49], which incorporates
differentiable optimization problems in a disciplined fashion. In a broader
sense, this methodology can be applied to differentiate through disciplined
convex programs by initially mapping them into cone programs [50], computing
the gradients, and subsequently mapping back to the original problem. It is
worth noting that the convex quadratic program (QP) can be differentiated
through the KKT conditions [16], which serve as equivalent conditions for
global optimality. According to the KKT conditions, at the optimal solution,
the gradient of the Lagrangian function with respect to the program’s input
and parameters must be zero. Consequently, by taking the partial derivative of
the Lagrangian function with respect to the input and extending it through the
chain rule to the program’s parameters, their gradients can be obtained. For a
generalized QP:

$\displaystyle\min_{x}\ \ \frac{1}{2}x^{T}Qx+q^{T}x,\quad\text{ s.t.}\ \
Ax=b\quad Gx\leq h,$ (3)

we can write the Lagrangian of the problem as:

$L(z,\nu,\lambda)=\frac{1}{2}z^{T}Qz+q^{T}z+\nu^{T}(Az-b)+\lambda^{T}(Gz-h)$
(4)

where $\nu$ are the dual variables on the equality constraints and
$\lambda\geq 0$ are the dual variables on the inequality constraints. The KKT
conditions for stationarity, primal feasibility, and complementary slackness
are:

$\displaystyle Qz^{\star}+q+A^{T}\nu^{\star}+G^{T}\lambda^{\star}$
$\displaystyle=0$ (5) $\displaystyle Az^{\star}-b$ $\displaystyle=0$
$\displaystyle D\left(\lambda^{\star}\right)\left(Gz^{\star}-h\right)$
$\displaystyle=0$

By differentiating these conditions, we can shape the Jacobian of the problem
as follows.

$\left[\begin{array}[]{l}d_{z}\\\ d_{\lambda}\\\
d_{\nu}\end{array}\right]=-\left[\begin{array}[]{ccc}Q&G^{T}D\left(\lambda^{\star}\right)&A^{T}\\\
G&D\left(Gz^{\star}-h\right)&0\\\
A&0&0\end{array}\right]^{-1}\left[\begin{array}[]{c}\left(\frac{\partial\ell}{\partial
z^{\star}}\right)^{T}\\\ 0\\\ 0\end{array}\right]$ (6)

Using the chain rule, we can get the derivatives of any loss function with
respect to any of the parameters in the QP.

##  III Distributionally Safe Reinforcement Learning under Model Uncertainty

Stochastic optimization solves problems with stochastic variables under a
prior distribution. However, in many real-world engineering applications, such
a prior distribution is often inaccurate or can change in the deployment. As a
result, Distributionally Robust Optimization (DRO) addresses optimization
problems with the presence of uncertainty, whose prior probability
distribution might shift. Compared to robust optimization admitting parameters
within certain sets or stochastic optimization admitting a distribution of
parameters, DRO takes a more conservative formulation to solve a more
challenging problem, where the distribution of the parameters shifts from the
prior. Then following the minimax formulation, DRO aims to minimize the
objective function with a worst-case distributional shift in the ambiguity
set. This general strategy provides an adaptable/ resilient solution that
exhibits strong performance in different/ unseen testing environments.

###  III-A Distributional safety under model uncertainty using the Wasserstein
metric: a bi-level problem

Control barrier functions have found application in maintaining the safety
constraints of control systems. In the realm of reinforcement learning, CBFs
can guarantee that the actions taken by an agent comply with safety
constraints, while maximizing a cumulative reward function. There are several
methods to incorporate the CBF in reinforcement learning. For example, CBF can
be incorporated in the reward function of the agent. Specifically, rewarding
safe actions and penalizing the agent for violating the safety constraints
leads to an agent that is biased towards safe actions. Another way of using
CBF in RL is to incorporate the CBF module as a safety shield for the RL
actions. The CBF will then monitor the actions and rectify any action that
violates the safety constraint to the closest safe action [51]. In the
presence of model uncertainty or noise, the conventional CBF approach does not
hold the same efficacy. Some work addressed robust CBFs by estimating the
dynamics noise and incorporating it directly into the CBF [52], while some
other works have used the conditional value at risk of the constraints to
build a safe CBF [53].

In this section, we aim to integrate distributional shift measured by a
Wasserstein metric into CBF to guarantee the safety of the system under the
model uncertainty. Consider the nonlinear system with additive uncertainty:

$\dot{x}(t)=f(x(t))+g(x(t))u(t)+\omega$ (7)

where $x\in\mathbb{R}^{n}$ and $u\in\mathbb{R}^{m}$ and $\omega\sim
p_{0}(\omega)$ are the states, control inputs, and the disturbance acting on
the system, respectively. Here $p_{0}(\omega)$ is the prior probability
distribution of the model uncertainty. The metric to measure the distance
between two probability distributions, i.e., distributional shift, is the
Wasserstein metric ($W_{d}(\bullet,\bullet)$), which is often the case in many
works addressing distributional robustness [54]. Comparing the Wasserstein
metric with the KL divergence [55], the most obvious advantage of the former
is symmetry (i.e., $W_{d}(p_{1},p_{2})=W_{d}(p_{2},p_{1})$), which however
does not hold for the latter in general. Then the ambiguity set of the
perturbed distribution $p(\omega)$ from the nominal distribution
$p_{0}(\omega)$ under a Wasserstein metric [54] is expressed as
$\mathcal{P}=\\{p(\omega)\in
p(\mathbb{W})|W_{d}\big{(}p(\omega),p_{0}(\omega)\big{)}\leq\rho\\}$, with
$\rho$ as the threshold of such a shift. Moreover, $\mathbb{W}$ is the support
of $p(M)$ and is assumed to be convex and closed, which is a common assumption
in related works.

The objective of safe reinforcement learning is to generate a control input
$u_{r}$ to achieve certain goals characterized by the reward function in the
MDP while satisfying safety constraints. The typical way for goal-reaching
robotic navigation is to drive a potential function $V(x)$ to be zero, with
$V(x)=\left\|x-x_{f}\right\|_{2}^{2}$. The RL policy will generate an action
without safety guarantee first as
$u_{\text{RL}}(t)=\mu\left(x_{t}\mid\theta^{\mu}\right)+\mathcal{N}_{t}$,
where $\mu(\bullet\mid\theta^{\mu})$ is a policy parameterized by deep neural
networks $\theta^{\mu}$ and $\mathcal{N}$ is a random process for promoting
exploration. Then the barrier function method [28] ensures that the action
$u_{\text{RL}}$ complies with safety constraints111The state $x$ and $s$, the
control/ action $u$ and $a$, terminologies in control theory and reinforcement
learning, are used interchangeably here.. With the presence of the model
uncertainty $\omega$ under distributional shift, the CBF of the safety
constraint $h(x)\geq 0$, defined in definition(1), can be rewritten using the
chain rule as follows:

$\displaystyle\sup_{u\in U}\inf_{\begin{subarray}{c}\omega\sim p(\omega)\\\
p(\omega)\in\mathcal{P}\end{subarray}}\bigg{[}\frac{\partial h(x)}{\partial
x}(\dot{x}(t))+\alpha(h(x))\bigg{]}\geqslant 0$ (8) $\displaystyle=$
$\displaystyle\sup_{u\in U}\inf_{\begin{subarray}{c}\omega\sim p(\omega)\\\
p(\omega)\in\mathcal{P}\end{subarray}}\bigg{[}\underbrace{\frac{\partial
h(x)}{\partial
x}(f(x)+g(x)u+\omega)+\alpha(h(x))}_{:=H(x,u,\omega)}\bigg{]}\geqslant 0$

where $\alpha(\bullet)$ is an extended class-$\mathcal{K}$ function. For
simplicity, we use a linear function $\kappa(\bullet)$ as the
class-$\mathcal{K}$ function. The infimum tries to evaluate the worst-case
safety constraint while the supremum to find a feasible control input such
that the safety constraint is still satisfied even with the worst
perturbation. While the worst-case criterion is adopted here (possibly with
over-conservatism), other ones such as the chance-constrained criterion are
compatible with this framework as well.

Using (8) as the safety shield to solve for the rectified action
$u_{\text{R}}$ under the worst-case distribution of $\omega$ leading to the
following formulation

$\displaystyle\min_{u_{r}\in[\underline{u},\bar{u}]}||u_{r}-u_{\text{RL}}||^{2}$
(9) $\displaystyle\text{ s.t.
}\inf_{p(\omega)\in\mathcal{P}}\mathbb{E}_{\omega\sim
p(\omega)}\Big{\\{}H(x,u,\omega)\big{|}W_{d}\big{(}p(\omega),p_{0}(\omega)\big{)}\leq\rho\Big{\\}}\geqslant
0$

which is a bi-level programming. The supreme in (8) disappears as part of the
feasibility of the high-level minimization in (9). However, addressing this
bi-level DRO problem can be challenging due to the high computational
complexity of its bi-level nature on top of an already challenging
distributionally constrained optimization problem at the lower level.

For the low-level DRO, obtaining the infimum over the expectation of the
distributions $p(\omega)$ has been proved difficult. Various techniques, such
as convex relaxations, scenario approximations, or sample-based methods, are
used to handle the computational challenges associated with DRO. [56, 46, 34].
That is because the low-level infimum problem requires a search for the worst-
case safety violation $H(x,u,\omega)$ within the infinite-dimensional
probability space $\mathcal{P}$, making it intractable to solve [57, 58]. As a
result, we will demonstrate how we can solve the low-level safety estimation
problem in an efficient way based on distributionally robust optimization
[47]. Using the Kantorovich duality [59, 60], the infimum in (9) can be
further transformed into a more tractable problem

$\displaystyle\min_{u_{r}\in[\underline{u},\bar{u}]}||u_{r}-u_{\text{RL}}||^{2}$
(10) $\displaystyle\text{ s.t. }\mathbb{E}_{\omega_{0}\sim
p_{0}(\omega)}\inf_{\omega}\Big{\\{}H(x,u,\omega)\big{|}d\big{(}\omega,\omega_{0}\big{)}\leq\rho_{d}\Big{\\}}\geqslant
0$

where $d(\omega,\omega_{0})=\|\omega-\omega_{0}\|_{p}^{2},p\geq 1$ denotes the
the “cost” for an adversary to perturb $\omega_{0}$ to $\omega$ [45]. This
equivalent dual reformulation [45, 60] allows for solving the infimum over the
noisy state $\omega\in\mathbb{R}^{n}$, a parametric finite-dimensional space
instead of over $p(w)$ in the infinite-dimensional probability space
$\mathcal{P}$ in (9). The radius $\rho_{d}$ can be evaluated with Monte-Carlo
simulation based on the original ambiguity set measured by the Wasserstein
metric. However, even with the complexity of DRO alleviated, the problem is
still bi-level, and hence in our case, we propose a single-level reformulation
to improve tractability without compromising safety.

###  III-B Reduce the bi-level learning to single-level: a differentiable
convex programming-based approach

To address the complexity of the bi-level optimization problem (10), we will
take advantage of differentiable programming to reduce the problem into a
single-level optimization problem. The basic idea is that for general
optimization problems, there is a trade-off between optimizing objectives and
satisfying constraints. Specifically here, when the dynamics perturbation
$\omega$ makes the constraint satisfaction more challenging, then the control
input will lean more on constraint satisfaction rather than optimizing the
objective. Eventually, the effect of the model uncertainty will be quantified
by the objective function values. For example, for a variable value not
satisfying the constraints, the objective value tends to go to infinity. As a
result, we move the infimum of the constraints to the maximization of the
objective function. For safe reinforcement learning, the loss function is the
negative of the discounted cumulative reward with the safety-proofed actions.
We decompose the bi-level problem in (10) into two parts: the safety
calibration to get $u_{r}$ followed by the update of $\omega$ to get the worst
case uncertainty. Specifically, for the former, we can start by solving for
the optimal action subject to the sample noise under a prior distribution,
i.e., $\omega_{0}\sim p_{0}(w)$ via the simple CBF-based quadratic program:

$\displaystyle\min_{u_{r}\in[\underline{u},\bar{u}]}$
$\displaystyle||u_{r}-u_{\text{RL}}||^{2}$ (11) s.t.
$\displaystyle\mathbb{E}_{\omega_{0}\sim
p_{0}(\omega)}H(x,u,\omega_{0})\geqslant 0$

Note that in one-step propagation, multiple samples $\omega_{0}$’s will be
sampled as parameters to evaluate the mean in (11). Then we can propagate our
plant using the rectified action $u_{r}$, and calculate the episodic loss
function as

$\mathcal{L}=-\sum_{k=1}^{T}\gamma^{k}{R}(x_{k},u_{r,k}).$ (12)

By the end of each episode, we can extract a sample of the noise $\omega_{0}$
and the gradient of the loss function with respect to the rectified action
$\frac{\partial\mathcal{L}}{\partial u_{r}}$. As discussed before, the model
uncertainty will try to maximize the loss function, which will be used to find
the worst-case perturbation $\omega$ by gradient-based algorithms. Therefore,
we can use the chain rule to obtain the gradient of the loss function with
respect to the noise $\omega$ as:

$\frac{\partial\mathcal{L}}{\partial\omega}=\frac{\partial\mathcal{L}}{\partial
u_{r}}\frac{\partial u_{r}}{\partial\omega}$ (13)

To obtain the second term, differentiable convex programming will be leveraged
for solving (11) and extract the gradient of the rectified action ($u_{r}$,
the variable) with respect to the noise ($\omega$, the parameter) as
$\frac{\partial u_{r}}{\partial\omega}$, evaluated at the samples
$\omega_{0}$’s for trajectory rollout. Note that the gradients flow through
the QP (11) so that the whole pipeline is end-to-end differentiable. Since we
aim to find the worst-case uncertainty dynamics within the ambiguity set
$\mathcal{B}=\\{\omega|\mathbb{E}_{\omega_{0}\sim
p_{0}(w)}d(\omega-\omega_{0})\leqslant\rho_{d}\\}$ by maximizing the loss
function, projected gradient ascent will be employed to update the $\omega$ as

$\omega\leftarrow{\mbox{\text{Proj}}}_{\mathcal{B}}\bigg{[}\omega+\alpha\frac{\partial\mathcal{L}}{\partial\omega}\bigg{]}.$
(14)

The projection of $\bar{\omega}$ onto a ball can be done analytically by
performing a change of basis to the $l_{2}$ ball towards the center
$\mathbb{E}(\omega_{0})$. Then projecting $\bar{\omega}$ onto the closest
point within the ball leads to

${\mbox{\text{Proj}}}_{\mathcal{B}}(\bar{\omega})=\mathbb{E}(\omega_{0})+\rho_{d}\frac{\bar{\omega}-\mathbb{E}(\omega_{0})}{\max(\rho_{d},||\bar{\omega}-\mathbb{E}(\omega_{0})||)}$
(15)

Fig.1 illustrates the workflow of the proposed approach. Algorithm1 summarizes
the overall distributional safe RL framework with DDPG [61] and the learnable
worst-case noise $\omega$ under a distributional shift. In summary, we
transformed bi-level programming for distributional safety into a sequential
problem: a convex quadratic program followed by a projected gradient ascent
learning module for the worst-case noise. Both parts are computationally cheap
and differentiable convex programming is employed for end-to-end
differentiability to update both policy parameters with constraints and worst-
case uncertainty.

Figure 1: Overview of the single-level, end-to-end differentiable convex
programming learning-based approach for distributionally safe reinforcement
learning using a Wasserstein metric. Black (dashed) lines represent data flow
and trajectory rollout, while colored solid lines represent gradient back-
propagation for updating both policies and worst-case noise. Algorithm 1
Distributionally Safe Reinforcement Learning

1:Require: Environment setting, learning rates $\alpha,\beta$, discount factor
$\gamma$, target network update rate $\tau$, and samples $\omega_{0}\sim
p_{0}(\omega)$.

2:Initialize critic network $Q\left(s,a\mid\theta^{Q}\right)$, actor
$\mu\left(s\mid\theta^{\mu}\right)$ with weights $\theta^{Q}$ and
$\theta^{\mu}$.

3:Initialize target network $Q^{\prime}$ and $\mu^{\prime}$ with weights
$\theta^{Q^{\prime}}\leftarrow\theta^{Q},\theta^{\mu^{\prime}}\leftarrow\theta^{\mu}$
.

4:Initialize replay buffer $\mathcal{D}$.

5:Initialize a random $\omega$.

6:for episode $=1,\ldots,M$ do

7: Initialize a random noise $\omega_{0}\sim p_{0}(\omega)$ from samples .

8: Initialize a random process $\mathcal{N}$ for action exploration.

9: Receive initial observation state $s_{1}$ .

10: for ${t}=1,\ldots,T$ do

11: Select action
$a_{t}=\mu\left(s_{t}\mid\theta^{\mu}\right)+\mathcal{N}_{t}$ according to the
current policy and exploration noise.

12: Get rectified action $a_{t_{R}}$ via (11) and samples $\omega_{0}$.

13: Execute action $a_{t_{R}}$ and observe reward $R_{t}$ and new state
$s_{t+1}$.

14: Store transition
$\left(s_{t},a_{t},a_{t_{R}},R_{t},s_{t+1},\omega_{0}\right)$ in
$\mathcal{D}$.

15: Sample a random mini-batch of ${N}$ transitions
$\left(s_{t},a_{t},a_{t_{R}},R_{t},s_{t+1},\omega_{0}\right)$ from
$\mathcal{D}$.

16: Update critic using learning rate $\beta$.

17: Update the actor $\theta^{\mu}$ using the gradient ascent with the sampled
gradient of the return in (12).

18: $\theta^{\mu}\leftarrow\theta^{\mu}+\alpha\nabla_{\theta^{\mu}}J(\theta).$

19: Update the target networks with rate $\tau$ .

20: $\theta^{\prime}\leftarrow\tau\theta+(1-\tau)\theta^{\prime}$ .

21: end for

22: Sample noise $\omega_{0}$ from all transitions
$\left(s_{t},a_{t},a_{t_{R}},R_{t},s_{t+1},\omega_{0}\right)$ in $\mathcal{D}$
and evaluate the loss $\mathcal{L}$ and its gradient
$\frac{\partial\mathcal{L}}{\partial\omega}$.

23: Update $\omega$ using the projected gradient ascent
$\omega\leftarrow{\mbox{\text{Proj}}}_{\mathcal{B}}\big{[}\omega+\alpha\frac{\partial\mathcal{L}}{\partial\omega}\big{]}$
with the gradeient evaluation in (13) and the projection in (15).

24:end for

25:Return: $\theta^{\mu},\theta^{Q},\omega$.

##  IV SIMULATIONS AND RESULTS

In this section, we assess the performance of the suggested safe reinforcement
learning approach using two cases of Dubin’s car, a first-order and second-
order system, and a simplified 3D quadcopter. Our goal is to compare the
performance of the deterministically learned policy in the presence of
dynamics uncertainty, with the distributionally safe learned policy. The
comparison is carried out under the condition that all other settings and
parameters are kept identical to ensure a consistent assessment. By carrying
out this comparison, we aim to show the effectiveness of the distributionally
safe approach in dealing with uncertain environments.

####  IV-1 First-Order Dubins Car

The first simulation is carried out using the first-order Dubins car
environment with the following kinematics(16).

$\left(\begin{array}[]{c}\dot{x}\\\ \dot{y}\\\
\dot{\theta}\end{array}\right)=\left[\begin{array}[]{ccc}\cos\theta&-\sin\theta&0\\\
\sin\theta&\cos\theta&0\\\
0&0&1\end{array}\right]\left(\begin{array}[]{c}u_{x}\\\ u_{y}\\\
u_{\theta}\end{array}\right),$ (16)

where $u_{x},u_{y},u_{\theta}$ are the input velocity along the $x$ axis, the
sideways input velocity, and the angular input velocity, respectively. To
reach its final destination $x_{f}$ from an initial state $x_{o}$, we use a
reward that penalizes the squared distance between the car and the goal state
multiplied by a coefficient as $d\left\|x-x_{f}\right\|_{2}^{2}$, and
penalizes every time step by a constant $s$ for minimum time goal-reaching.
Hence, the reward is defined as $R=-d\left\|x-x_{f}\right\|_{2}^{2}-s$, with
$d>0$ and $s\geq 0$.

Fig. 2 presents a comparison between trajectories generated by a
deterministically learned policy, and distributionally safe policies before
and after training converges.

Figure 2: First-order Dubin’s car trajectories with the presence of model
uncertainty for deterministic policy, the DSRL policies before and after
convergence.

We can see that although the deterministic policy is trying to avoid the
obstacle by swerving in the same shape as the obstacle, it still violated the
safety constraint because it could not take into account the shift/ noise
present in the dynamics of the system. The DSRL policy adapted to the noise
after several episodes of training and started avoiding the obstacle while
keeping a very safe margin.

Figure 3: DSRL returns over the learning period for 3 random seeds for the
Dubin’s car: Left: first-order; Right: second order.

In Fig. 3 (left), we can see that although the learning of the DSRL policy is
fluctuating between episodes, the general trend of the cumulative reward
increases over the total number of episodes, demonstrating a promising
learning curve.

####  IV-2 Second-Order Dubins Car

To demonstrate our approach on a higher-order system, we use the second-order
Dubin’s car with the following dynamics:

$\left(\begin{array}[]{c}\ddot{x}\\\ \ddot{y}\\\
\ddot{\theta}\end{array}\right)=\left[\begin{array}[]{ccc}\cos\theta&-\sin\theta&0\\\
\sin\theta&\cos\theta&0\\\
0&0&1\end{array}\right]\left(\begin{array}[]{c}u_{x}\\\ u_{y}\\\
\tau_{c}\end{array}\right),$ (17)

We also have the following adjusted reward function to penalize the velocities
at the final destination for learning to brake as well:
$R=-d\left\|x-x_{f}\right\|_{2}^{2}-b\left\|v-v_{f}\right\|_{2}^{2}-s$. Fig. 4
and Fig.3(right) illustrate the trajectories of the benchmark and two episodes
from the learning episodes, and the cumulative rewards. In the second-order
case, we can see that the benchmark misses the finish line by a slight margin
due to the influence of noise. All other behaviors are similar to the first-
order case.

Figure 4: Second-order Dubin’s car trajectories with the presence of model
uncertainty for deterministic policy, the DSRL policies before and after
convergence.

####  IV-3 Quadcopter

To inspect the performance of the DSRL on a system with more complex dynamics
and safety constraints, we design high-level controllers for a heading-locked
quadcopter environment with the following dynamics

$\left(\begin{array}[]{c}\dot{x}\\\ \dot{y}\\\
\dot{z}\end{array}\right)=\left(\begin{array}[]{c}T\sin{\theta}\\\
T\cos{\theta}\sin{\phi}\\\ T\cos{\theta}\cos{\phi}-g\end{array}\right)$ (18)

where $x,y,z$ are the inertial displacements of the quadcopter and
$\theta,\phi,T$ are the pitch, roll, and total thrust, respectively. While
training, the quadcopter is penalized by the distance between the initial and
terminal positions: $R=-\left\|x-x_{f}\right\|_{2}^{2}$. In this section, we
simulate a landing scenario of the quadcopter within a prescribed glide slope
defined by a cone with a half angle $\delta$. The safety constraint can be
written concisely as:

$r_{I}^{T}M_{gs}r_{I}\geqslant 0$ (19)

where
$M_{gs}=[[-\cot(\delta_{gs})^{2},0,0]^{T},[0,-\cot(\delta_{gs})^{2},0]^{T},[0,0,1]^{T}]$
and $r_{I}=[x,y,z]$ is the inertial position of the quadcopter.

Figure 5: Quadcopter trajectories with the presence of model uncertainty for
deterministic policy and the DSRL policy after convergence.

Fig.5 presents the trajectories generated by the DSRL at the start and end of
the learning, and the one generated by the deterministic policy subject to
dynamic noise. The learned $\omega$ outperforms the deterministic policy, as
it always keeps a safe margin from the constraint, hence, maintaining safety
throughout the whole trajectory.

##  V CONCLUSIONS

This work addresses the challenges of tractable safe reinforcement learning
under distributional shift. We employ differentiable convex programming along
with distributinoally robust optimization to enable the safe learning over
unpredictable probability distributions of the model’s uncertainty in a
tractable manner. We evaluate our approach on first and second-order Dubin’s
car and a simplified quad-copter model and compare our results with the
deterministic policies’ performance, where our approach outperforms the
benchmark policies regarding safety guarantee under model uncertainty.

## References

  * [1] S. Li, Y. Wu, X. Cui, H. Dong, F. Fang, and S. Russell, “Robust multi-agent reinforcement learning via minimax deep deterministic policy gradient,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, no. 01, 2019, pp. 4213–4220.
  * [2] K. Zhang, Z. Yang, and T. Basar, “Policy optimization provably converges to nash equilibria in zero-sum linear quadratic games,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [3] C. Sun, D.-K. Kim, and J. P. How, “Romax: Certifiably robust deep multiagent reinforcement learning via convex relaxation,” _arXiv preprint arXiv:2109.06795_ , 2021.
  * [4] B. Yang, L. Zheng, L. J. Ratliff, B. Boots, and J. R. Smith, “Stackelberg maddpg: Learning emergent behaviors via information asymmetry in competitive games,” 2022.
  * [5] Z. Zhou and H. Xu, “Decentralized adaptive optimal tracking control for massive autonomous vehicle systems with heterogeneous dynamics: A stackelberg game,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 32, no. 12, pp. 5654–5663, 2021.
  * [6] N. Lauffer, M. Ghasemi, A. Hashemi, Y. Savas, and U. Topcu, “No-regret learning in dynamic stackelberg games,” _arXiv preprint arXiv:2202.04786_ , 2022.
  * [7] Y. Bai, C. Jin, H. Wang, and C. Xiong, “Sample-efficient learning of stackelberg equilibria in general-sum games,” _Advances in Neural Information Processing Systems_ , vol. 34, 2021.
  * [8] W. Jin, S. Mou, and G. Pappas, “Safe pontryagin differentiable programming,” _Advances in Neural Information Processing Systems_ , vol. 34, 2021.
  * [9] Y. Liu, J. Ding, and X. Liu, “Ipo: Interior-point policy optimization under constraints,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 04, 2020, pp. 4940–4947.
  * [10] J. Achiam, D. Held, A. Tamar, and P. Abbeel, “Constrained policy optimization,” in _Proceedings of the 34th International Conference on Machine Learning-Volume 70_. JMLR. org, 2017, pp. 22–31.
  * [11] D. Bertsekas and J. Tsitsiklis, _Parallel and distributed computation: numerical methods_. Athena Scientific, 2015.
  * [12] P. Geibel and F. Wysotzki, “Risk-sensitive reinforcement learning applied to control under constraints,” _Journal of Artificial Intelligence Research_ , vol. 24, pp. 81–108, 2005.
  * [13] Y. Chow, M. Ghavamzadeh, L. Janson, and M. Pavone, “Risk-constrained reinforcement learning with percentile risk criteria,” _The Journal of Machine Learning Research_ , vol. 18, no. 1, pp. 6070–6120, 2017.
  * [14] Y. Chow, O. Nachum, E. Duenez-Guzman, and M. Ghavamzadeh, “A lyapunov-based approach to safe reinforcement learning,” in _Advances in neural information processing systems_ , 2018, pp. 8092–8101.
  * [15] A. Stooke, J. Achiam, and P. Abbeel, “Responsive safety in reinforcement learning by pid lagrangian methods,” _arXiv preprint arXiv:2007.03964_ , 2020.
  * [16] B. Amos and J. Z. Kolter, “Optnet: Differentiable optimization as a layer in neural networks,” _arXiv preprint arXiv:1703.00443_ , 2017.
  * [17] H. Zhang, H. Chen, C. Xiao, B. Li, M. Liu, D. Boning, and C.-J. Hsieh, “Robust deep reinforcement learning against adversarial perturbations on state observations,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 21 024–21 037, 2020.
  * [18] K. Zhang, T. Sun, Y. Tao, S. Genc, S. Mallya, and T. Basar, “Robust multi-agent reinforcement learning with model uncertainty,” _Advances in Neural Information Processing Systems_ , vol. 33, pp. 10 571–10 583, 2020.
  * [19] L. Brunke, M. Greeff, A. W. Hall, c. Yuan, S. Zhou, J. Panerati, and A. P. Schoellig, “Safe learning in robotics: From learning-based control to safe reinforcement learning,” _Annual Review of Control, Robotics, and Autonomous Systems_ , vol. 5, 2021.
  * [20] B. Singh, R. Kumar, and V. P. Singh, “Reinforcement learning in robotic applications: a comprehensive survey,” _Artificial Intelligence Review_ , pp. 1–46, 2021.
  * [21] J. Garcıa and F. Fernández, “A comprehensive survey on safe reinforcement learning,” _Journal of Machine Learning Research_ , vol. 16, no. 1, pp. 1437–1480, 2015.
  * [22] J. Moos, K. Hansel, H. Abdulsamad, S. Stark, D. Clever, and J. Peters, “Robust reinforcement learning: A review of foundations and recent advances,” _Machine Learning and Knowledge Extraction_ , vol. 4, no. 1, pp. 276–315, 2022.
  * [23] H. Eghbal-zadeh, F. Henkel, and G. Widmer, “Learning to infer unseen contexts in causal contextual reinforcement learning,” in _Self-Supervision for Reinforcement Learning Workshop-ICLR 2021_ , 2021.
  * [24] K. Rakelly, A. Zhou, C. Finn, S. Levine, and D. Quillen, “Efficient off-policy meta-reinforcement learning via probabilistic context variables,” in _International conference on machine learning_. PMLR, 2019, pp. 5331–5340.
  * [25] J. Choi, F. Castañeda, C. J. Tomlin, and K. Sreenath, “Reinforcement learning for safety-critical control under model uncertainty, using control lyapunov functions and control barrier functions,” _arXiv preprint arXiv:2004.07584_ , 2020.
  * [26] L. Zheng, Y. Shi, L. J. Ratliff, and B. Zhang, “Safe reinforcement learning of control-affine systems with vertex networks,” in _Learning for Dynamics and Control_. PMLR, 2021, pp. 336–347.
  * [27] R. Cheng, G. Orosz, R. M. Murray, and J. W. Burdick, “End-to-end safe reinforcement learning through barrier functions for safety-critical continuous control tasks,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 33, 2019, pp. 3387–3395.
  * [28] A. D. Ames, X. Xu, J. W. Grizzle, and P. Tabuada, “Control barrier function based quadratic programs for safety critical systems,” _IEEE Transactions on Automatic Control_ , vol. 62, no. 8, pp. 3861–3876, 2016.
  * [29] F. Berkenkamp, M. Turchetta, A. Schoellig, and A. Krause, “Safe model-based reinforcement learning with stability guarantees,” in _Advances in neural information processing systems_ , 2017, pp. 908–918.
  * [30] C. Sun, D.-K. Kim, and J. P. How, “Fisar: Forward invariant safe reinforcement learning with a deep neural network-based optimizer,” in _2021 IEEE International Conference on Robotics and Automation (ICRA)_. IEEE, 2021, pp. 10 617–10 624.
  * [31] G. Shi, K. Azizzadenesheli, M. O’Connell, S.-J. Chung, and Y. Yue, “Meta-adaptive nonlinear control: Theory and algorithms,” _Advances in Neural Information Processing Systems_ , vol. 34, pp. 10 013–10 025, 2021.
  * [32] M. O’Connell, G. Shi, X. Shi, K. Azizzadenesheli, A. Anandkumar, Y. Yue, and S.-J. Chung, “Neural-fly enables rapid learning for agile flight in strong winds,” _Science Robotics_ , vol. 7, no. 66, p. eabm6597, 2022.
  * [33] B. P. Van Parys, D. Kuhn, P. J. Goulart, and M. Morari, “Distributionally robust control of constrained stochastic systems,” _IEEE Transactions on Automatic Control_ , vol. 61, no. 2, pp. 430–442, 2015.
  * [34] I. Yang, “Wasserstein distributionally robust stochastic control: A data-driven approach,” _IEEE Transactions on Automatic Control_ , vol. 66, no. 8, pp. 3863–3870, 2020.
  * [35] A. Hakobyan and I. Yang, “Wasserstein distributionally robust motion control for collision avoidance using conditional value-at-risk,” _IEEE Transactions on Robotics_ , vol. 38, no. 2, pp. 939–957, 2021.
  * [36] A. Bahari Kordabad, R. Wisniewski, and S. Gros, “Safe reinforcement learning using wasserstein distributionally robust mpc and chance constraint,” 2022.
  * [37] J. Coulson, J. Lygeros, and F. Dörfler, “Distributionally robust chance constrained data-enabled predictive control,” _IEEE Transactions on Automatic Control_ , vol. 67, no. 7, pp. 3289–3304, 2021.
  * [38] P. Coppens and P. Patrinos, “Data-driven distributionally robust mpc for constrained stochastic systems,” _IEEE Control Systems Letters_ , vol. 6, pp. 1274–1279, 2021.
  * [39] A. Z. Ren and A. Majumdar, “Distributionally robust policy learning via adversarial environment generation,” _IEEE Robotics and Automation Letters_ , vol. 7, no. 2, pp. 1379–1386, 2022.
  * [40] D. Morrison, P. Corke, and J. Leitner, “Egad! an evolved grasping analysis dataset for diversity and reproducibility in robotic manipulation,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 3, pp. 4368–4375, 2020.
  * [41] D. Wang, D. Tseng, P. Li, Y. Jiang, M. Guo, M. Danielczuk, J. Mahler, J. Ichnowski, and K. Goldberg, “Adversarial grasp objects,” in _2019 IEEE 15th International Conference on Automation Science and Engineering (CASE)_. IEEE, 2019, pp. 241–248.
  * [42] M. Xu, P. Huang, Y. Niu, V. Kumar, J. Qiu, C. Fang, K.-H. Lee, X. Qi, H. Lam, B. Li, _et al._ , “Group distributionally robust reinforcement learning with hierarchical latent variables,” _arXiv preprint arXiv:2210.12262_ , 2022.
  * [43] N. Kallus, X. Mao, K. Wang, and Z. Zhou, “Doubly robust distributionally robust off-policy evaluation and learning,” in _International Conference on Machine Learning_. PMLR, 2022, pp. 10 598–10 632.
  * [44] L. Shi and Y. Chi, “Distributionally robust model-based offline reinforcement learning with near-optimal sample complexity,” _arXiv preprint arXiv:2208.05767_ , 2022.
  * [45] A. Sinha, H. Namkoong, and J. Duchi, “Certifiable distributional robustness with principled adversarial training,” _arXiv preprint arXiv:1710.10571_ , vol. 2, 2017.
  * [46] J. Blanchet and K. Murthy, “Quantifying distributional model risk via optimal transport,” _Mathematics of Operations Research_ , vol. 44, no. 2, pp. 565–600, 2019.
  * [47] H. Rahimian and S. Mehrotra, “Distributionally robust optimization: A review,” _arXiv preprint arXiv:1908.05659_ , 2019.
  * [48] W. Xiao and C. Belta, “High-order control barrier functions,” _IEEE Transactions on Automatic Control_ , vol. 67, no. 7, pp. 3655–3662, 2021.
  * [49] A. Agrawal, B. Amos, S. Barratt, S. Boyd, S. Diamond, and J. Z. Kolter, “Differentiable convex optimization layers,” _Advances in neural information processing systems_ , vol. 32, 2019.
  * [50] A. Agrawal, S. Barratt, S. Boyd, E. Busseti, and W. M. Moursi, “Differentiating through a cone program,” _arXiv preprint arXiv:1904.09043_ , 2019.
  * [51] A. E. Chriat and C. Sun, “On the optimality, stability, and feasibility of control barrier functions: An adaptive learning-based approach,” _arXiv preprint arXiv:2305.03608_ , 2023.
  * [52] E. Daş, S. X. Wei, and J. W. Burdick, “Robust control barrier functions with uncertainty estimation,” _arXiv preprint arXiv:2304.08538_ , 2023.
  * [53] A. E. Chriat and C. Sun, “Wasserstein distributionally robust control barrier function using conditional value-at-risk with differentiable convex programming,” _Accepted to AIAA SciTech 2024_ , 2024.
  * [54] C. R. Givens and R. M. Shortt, “A class of wasserstein metrics for probability distributions.” _Michigan Mathematical Journal_ , vol. 31, no. 2, pp. 231–240, 1984.
  * [55] J. M. Joyce, “Kullback-leibler divergence,” in _International encyclopedia of statistical science_. Springer, 2011, pp. 720–722.
  * [56] R. T. Rockafellar, S. Uryasev, _et al._ , “Optimization of conditional value-at-risk,” _Journal of risk_ , vol. 2, pp. 21–42, 2000.
  * [57] Y. Wang, G. Zhang, and J. Ba, “On solving minimax optimization locally: A follow-the-ridge approach,” _arXiv preprint arXiv:1910.07512_ , 2019.
  * [58] K. K. Thekumparampil, P. Jain, P. Netrapalli, and S. Oh, “Efficient algorithms for smooth minimax optimization,” _Advances in Neural Information Processing Systems_ , vol. 32, 2019.
  * [59] M. Berger, B. Eckmann, P. Harpe, F. Hirzebruch, N. Hitchin, L. Hörmander, A. Kupiainen, G. Lebeau, M. Ratner, D. Serre, _et al._ , _Optimal transport: old and new_. Springer, 2009.
  * [60] Z. Wang, X. Wang, L. Shen, Q. Suo, K. Song, D. Yu, Y. Shen, and M. Gao, “Meta-learning without data via wasserstein distributionally-robust model fusion,” in _Uncertainty in Artificial Intelligence_. PMLR, 2022, pp. 2045–2055.
  * [61] T. P. Lillicrap, J. J. Hunt, A. Pritzel, N. Heess, T. Erez, Y. Tassa, D. Silver, and D. Wierstra, “Continuous control with deep reinforcement learning,” _arXiv preprint arXiv:1509.02971_ , 2015.


# Ionic liquid gating of SrTiO3 lamellas fabricated with a focused ion beam

Evgeny Mikheev  Department of Physics, Stanford University, Stanford, CA,
94305, USA  Stanford Institute for Materials and Energy Sciences, SLAC
National Accelerator Laboratory, Menlo Park, CA 94025, USA  Tino Zimmerling
Max-Planck-Institute for Chemical Physics of Solids, 01187 Dresden, Germany
Amelia Estry  Laboratory of Quantum Materials (QMAT), Institute of Materials
(IMX), École Polytechnique Fédérale de Lausanne (EPFL), Lausanne, Switzerland.
Max-Planck-Institute for Chemical Physics of Solids, 01187 Dresden, Germany
Philip J. W. Moll  Laboratory of Quantum Materials (QMAT), Institute of
Materials (IMX), École Polytechnique Fédérale de Lausanne (EPFL), Lausanne,
Switzerland.  Max-Planck-Institute for Chemical Physics of Solids, 01187
Dresden, Germany  David Goldhaber-Gordon  Department of Physics, Stanford
University, Stanford, CA, 94305, USA  Stanford Institute for Materials and
Energy Sciences, SLAC National Accelerator Laboratory, Menlo Park, CA 94025,
USA

###### Abstract

In this work, we combine two previously-incompatible techniques for defining
electronic devices: shaping three-dimensional crystals by focused ion beam
(FIB), and two-dimensional electrostatic accumulation of charge carriers. The
principal challenge for this integration is nanometer-scale surface damage
inherent to any FIB-based fabrication. We address this by using a sacrificial
protective layer to preserve a selected pristine surface. The test case
presented here is accumulation of 2D carriers by ionic liquid gating at the
surface of a micron-scale SrTiO3 lamella. Preservation of surface quality is
reflected in superconductivity of the accumulated carriers. This technique
opens new avenues for realizing electrostatic charge tuning in materials that
are not available as large or exfoliatable single crystals, and for patterning
the geometry of the accumulated carriers.

††preprint: APS/123-QED

## Introduction

Top-down fabrication using focused ion beam (FIB) milling has recently emerged
as a powerful technique to explore physical properties of quantum materials
[1, 2]. Selective material removal by FIB enables definition of intricate
device structures with features as small as hundreds of nanometers [3, 4],
with explicit control of three-dimensional shape [5, 6] and orientation of
conduction channel with respect to crystal axes [7, 8]. Successful
implementations of this technique have so far been in metallic compounds:
high-$T_{c}$ superconductors [9, 4], high-conductivity delafossites [3], heavy
fermion compounds [7], and arsenide semimetals [5, 10, 11]. The starting
materials are bulk-grown single crystals, where top-down patterning by
conventional lithography techniques is challenging due to chemical sensitivity
and/or unavailability of high quality thin films.

At present, a major constraint on FIB-based patterning techniques is from
surface damage [12, 13, 1]. Even very selective removal of material damages
the entire exposed surface: typically a nanometer-scale surface layer is
rendered partially amorphous and/or contains a high density of point defects.
Though the thickness of the damage layer can be reduced by selection of ion
species and acceleration voltages, it is rarely possible to reduce it below 2
nm [14]. For micron-scale structures with metallic bulk conduction such
surface disorder is negligible, but for semiconductors it could screen the
effect of electrostatic gates and dramatically-degrade mobility.

In this paper, we demonstrate voltage-driven accumulation of a superconducting
2D electron gas (2DEG) in a device defined by FIB. This combination is enabled
by protecting a pristine surface from FIB damage with a sacrificial conductive
layer that is removable in water. Electrons are subsequently accumulated at
the preserved channel surface by ionic liquid gating. Here, this method is
implemented on SrTiO3, a material where electron doping and ionic liquid
gating are both well known to drive a transition from an insulator to a two-
dimensional superconducting metal [15, 16, 17]. Observation of robust surface
charge transport and superconductivity in our device demonstrates feasibility
of integrating the FIB patterning toolbox with electrostatically tunable
devices.

Figure 1:  Schematic illustration of the fabrication procedure steps, see text
for details.

## Methods

The fabrication flow described in this section aims to demonstrate that a
first FIB-based patterning stage can be combined with subsequent conventional
lithography techniques, ionic liquid gating, and cryogenic transport
measurements. The account of each step of fabrication is accompanied by a
schematic illustration (Fig. 1.) Fig. 2 presents device images following
selected fabrication steps.

  * •

Step 1: A (001) Ti-terminated surface is prepared by soaking SrTiO3 single-
crystal substrates (MTI) for 20 minutes in de-ionized water heated to
approximately 65 °C, and annealing for 2 hours at 1000 °C in Ar/O2.

  * •

Step 2: A sacrificial layer of 200 nm Ge is deposited by e-beam evaporation.

  * •

Step 3: The crystal is cleaved using a diamond scriber to create a smooth side
surface.

  * •

Step 4: A rectangular lamella is defined by a Xenon FIB. The first cut is from
the cleaved side of the crystal, defining the bottom surface of the lamella
(“FIB 1” in Fig. 1).

  * •

Step 5: The second FIB cut is from the top, Ge-coated surface, defining the
side walls of the lamella (“FIB 2” in Fig. 1). Bridges on the sides are left
to hold the sample in place as seen in Fig. 2a. Next, the FIB is used to
polish the bottom surface of the lamella (same geometry as in the previous
step, “FIB 1” in Fig. 1), creating a more uniform rectangular sample and
allowing the surface of the Ge-coated STO to be flush with the SU8. For the
final FIB process, one of the bridges is cut and the other thinned down such
that the sample can be easily removed.

  * •

Step 6: The lamella is broken out of the bulk crystal using an ex-situ
micromanipulator. SU8 photoresist (GM1060, Gersteltec) is spun for 40 seconds
at 4000 rpm on a 5x5 mm oxide-coated silicon chip; silicon is easily
available, but other substrates would likely work similarly. With SU8 still in
the viscous liquid state, the lamella is slowly transferred onto its surface.

  * •

Step 7: The chip is left in air for 10 hours, heated in air to 65 °C for 5
minutes, heated to 95 °C for 10 minutes, then cooled to room temperature. This
process decreases the risk of the SU8 reflowing around the sample. Finally, to
solidify the SU8, the chip is baked at 200 °C for 1 hour. All ramps are done
at a rate of 2 °C per minute.

  * •

Step 8: Removal of the sacrificial Ge layer. 5 minute exposure to an O2 plasma
for residue descum, followed by a 4 hour soak in de-ionized water heated to
approximately 65 °C.

  * •

Step 9: First e-beam lithography (EBL) step. Deposition of alignment marks
around the device via liftoff of 15/60 nm of Ti/Au. This and subsequent
patterning steps are performed using spun-on PMMA resist, a 30 kV EBL system,
and a chilled 3:1 water:isopropyl alcohol mixture as a developer. Lift-off
steps are performed with acetone followed by an isopropyl alcohol rinse.

  * •

Step 10: Second EBL step. Ohmic contact definition by ion milling (900 V beam
voltage, 5 minutes) and lift-off of e-beam evaporated 10/40 nm of Ti/Au using
the same patterned resist layer. A Hall bar-like contact geometry is defined,
with two large source and drain contacts and six voltage probes.

  * •

Step 11: Third EBL step. Contact line, bond pads and side-gate deposition via
liftoff of 10/200 nm of Ti/Au.

  * •

Step 12: The lamella device and the side gate are covered by an ionic liquid
DEME-TFSI, and loaded into a dilution refrigerator system for charge transport
characterization.

Figure 2: (a) Scanning electron microscope image taken during FIB machining of
the lamella. (b-f) Optical images of the devices after fabrication steps 7, 8,
10, 11, 12 (see text). (g) Height image taken with a 3D laser scanning
confocal microscope (LSCM). (h) Comparison of height profile from LSCM and
AFM. Both are taken near the middle of the lamella, as indicated by the dashed
line in (g). (i) AFM image of the STO surface before sacrificial layer
deposition, and (j) after its removal, on the surface of the lamella. (k)
Perspective view of an LSCM height profile of the finished device.

The overarching idea behind this process is to use hardbaked SU8 resist as a
new planarized substrate for a micron-scale FIB-defined object that would
otherwise be too small for conventional lithography techniques. After
hardbaking, crosslinked SU8 is chemically robust to solvents, developers, and
typical e-beam doses required for PMMA patterning. It is also sufficiently
robust mechanically to withstand basic handling during device fabrication, and
subsequent repeated thermal cycling to cryogenic temperatures.

One key enabler for this process is the resulting device surface topography,
which is approximately flush between the FIB-defined lamella and the
surrounding SU8. Remarkably, micron scale objects placed on the surface of
liquid SU8 typically do not drown or float on their bottom surface, but
instead partially submerge down to the exposed top surface. This is
illustrated in Fig. 1 in the device cross section for step 6 onwards. Images
of actual device surface topography after hardbaking are shown in figures 2g
and 2h, using standard atomic force microscopy (AFM) and 3D laser scanning
confocal microscopy (LSCM, Keyence VK-X250). The resulting boundary between
SU8 and SrTiO3 can easily be crossed with patterned 210 nm thick contact metal
lines. In absence of the SU8 layer, metal lines would need to climb over a
micron-scale cliff between the lamella top surface and the underlying solid
substrate, which is challenging to realize robustly with reasonable metal
thicknesses. Additionally, side wall contacts for some materials are likely to
get shorted by conductive amorphous surface layers (see e.g. [10]).

A second key enabler is the use of a sacrificial Ge layer to preserve the
SrTiO3 surface from damage during the FIB milling steps. Ge satisfies two
necessary criteria: 1) it is sufficiently robust to withstand secondary
exposure to FIB (photoresist, in contrast, would be milled away too rapidly),
2) it can be selectively removed by simply soaking the device in water [18].
The removal progress can be tracked optically (Fig. 2b,c) as the lamella
surface color transitions from metallic (Ge) to transparent (SrTiO3). The
optical contrast on the device channel in Fig. 2c-e is from the lamella back
side topography. We found that for devices exposed to FIB, heating the water
above room temperature was necessary for complete removal, likely due to
implantation of resputtered material by FIB into Ge. The AFM images in Fig.
2i,j present a first indication that the sacrificial layer is successful in
shielding the intended device surface from FIB damage. The originally prepared
SrTiO3 surface has a terrace structure with unit cell steps resulting from
substrate miscut angle. This morphology is successfully preserved throughout
Ge deposition, FIB milling and Ge removal. The difference seen in step spacing
is a consequence of pre-existing miscut angle inhomogeneity across the 5x5 mm
substrate.

For transport measurements, the electrical contact between gold bond pads and
gold wire was made by cured conductive silver epoxy. The ionic liquid
Diethylmethyl(2-methoxyethyl)ammonium bis(trifluoromethylsulfonyl)imide (DEME-
TFSI) was deposited to cover the lamella and the adjacent side gate. The
device was then loaded into a dilution refrigerator and vacuum pumped
overnight to reduce ambient water contamination in the ionic liquid.

## Results

In the following, we present a detailed analysis of the superconducting
transition in the FIB-prepared slabs. This is important for two reasons.
First, the observed superconducting and normal state transport behavior
compare well to prior measurements on ionic liquid gated SrTiO3 surfaces in
devices fabricated without FIB [15, 19, 20, 21, 22, 23, 24, 25] and to the
closely-related 2DEG heterostructure systems SrTiO3/LaAlO3 and SrTiO3/LaTiO3
[26, 27, 28, 29, 30, 31, 32, 33, 34, 35]. These similarities leave no doubt
that the main goal of pristine surface preservation was accomplished in this
fabrication approach. Second, the SrTiO3 sample dimensions
(3$\times$15$\times$54 $\mu$m) are far smaller than in any previous SrTiO3
device used for electrical transport studies. In this new regime, local
strains fields can play an important role [36]. Moreover, these dimensions are
comparable to the typical size of tetragonal domains, which spontaneously form
in SrTiO3 below a cubic-to-tetragonal transition near 100 K [37] and are known
to affect normal state [38] and superconducting transport [39, 40].

Figure 3:  Gate tuning of the charge transport in the lamella device channel.
(a) Two-probe measurement near room temperature of DC current through the
channel ($I_{\text{DC}}$) and leakage current to the gate electrode through
the ionic liquid ($I_{\text{GIL}}$), as a function of gate voltage
($V_{\text{GIL}}$). (b) Temperature-dependent four probe measurements of
device channel resistance, taken at different fixed values of
$V_{\text{GIL}}$. For $V_{\text{GIL}}=$ 2.5 V, a two-probe measurement in the
same configuration as (a) is shown.

The gate voltage $V_{\text{GIL}}$ is applied on the side gate near room
temperature, polarizing the ionic liquid and accumulating a 2DEG on the
surface of the lamella, as illustrated in step 12 of Fig. 1. Fig. 3a shows the
2-terminal measurement of DC current $I_{\text{DC}}$ through the device
channel, measured while ramping $V_{\text{GIL}}$ up. A constant DC source bias
$V_{\text{SRC}}=$ 100 mV is applied to one of the ohmic contacts, while all
other contacts are grounded. Above $V_{\text{GIL}}=$ 1.2 V, carrier density
accumulation in the channel results in a gradual transition from an insulator
to a robust conductor. The device is then cooled down at constant
$V_{\text{GIL}}$. The ionic liquid freezes near 220 K, fixing the 2DEG state
until it is thermally cycled back to near room temperature for further
adjustment of $V_{\text{GIL}}$. The temperature dependence of 4-terminal
channel resistance was measured upon heating of the cryostat (Fig. 3b). For
$V_{\text{GIL}}$ above 3 V, the channel becomes metallic. At $V_{\text{GIL}}=$
3.1 V, a Kondo-like resistance upturn saturating at low temperature is seen,
similarly to [19, 21].

At $V_{\text{GIL}}=$ 3.5 and 4 V, robust metallicity is observed. The residual
resistance ratio between 200 K and base temperature is 20, and $\approx 45$ if
we extrapolate $R(T)$ to room temperature. These ratios imply low-temperature
carrier mobility within the 100-1000 cm2/Vs range typically observed in ionic
liquid gated SrTiO3 [15, 19, 25]. The cusp feature near 80 K is typical for
SrTiO3-based 2DEGs measured upon heating, and is usually attributed to charge
trapping and de-trapping at domain walls and/or point defects [30, 31, 35].

A superconducting transition is observed in the mK range (Fig. 3b,c). Taking
the midpoint of the resistance decrease as a measure of transition temperature
$T_{c}$, it is 400 and 380 mK for $V_{\text{GIL}}=$ 3.5 and 4 V, respectively.
400 mK is the upper bound of previously reported $T_{c}$ in ionic liquid gated
SrTiO3 [15, 20, 22, 23, 24, 25]. It is also notably higher than any such
devices previously fabricated in the authors’ laboratory (without FIB), which
typically have $T_{c}<$ 350 mK [22, 23, 24, 25]. Superconductivity in SrTiO3
is sensitive to strain, likely by tuning proximity to ferroelectric quantum
criticality [41, 42, 43]. The unexpectedly high $T_{c}$ in a micron-scale
lamella device may thus reflect local strain fields from differential thermal
expansion with the surrounding SU8 and the underlying silicon substrate [36].

Figure 4: Suppression of superconductivity with DC current $I_{\text{DC}}$,
temperature $T$, and perpendicular field $B$. (a) 4 probe DC measurement of
current-voltage non-linearity in the device channel. Purple dashed lines are
power law fits, with the red solid line indicating $a=3$ exponent criterion
for a BKT transition. (b) $T$ dependence of the exponent $a$ from power law
fits in (a). Solid line is the $T$ dependence of channel resistance (from Fig.
3b), scaled by a factor of 10 for comparison with $a(T)$. Dotted lines show
the transition point criteria $a=3$ and $R/R_{N}=0.5$. Dashed red line is a
fit to equation 1. (c) $T$ dependence of the critical current and a fit to
equation 2. (d) Same measurement as in (a), but as a function of $B$ at 32 mK.
(e) Critical current in (d) and a fit to equation 3. (f) 4 probe AC resistance
measurement as a function of $B$ and $T$. (g) $T$ dependence of the critical
field, normalized to base temperature value. Black dotted line is the 3D BCS
lineshape. Red dotted and solid lines are WHH model lineshapes [44], with
$\alpha=$ 0 and 1.8. $V_{\text{GIL}}=$ 4 V in (a-e), 3.5 V in (f, g).

Fig. 4 illustrates the critical scales for suppression of superconductivity
with DC current $I_{\text{DC}}$ and magnetic field $B$. Fig. 4a shows a strong
non-linearity in 4-terminal DC measurements. $I_{\text{DC}}$ is the DC current
excitation sourced through the large source and drain contacts, and
$V_{\text{DC}}$ the DC bias measured by adjacent voltage probes on the
channel. This non-linearity follows a power law
$I_{\text{DC}}\sim(V_{\text{DC}})^{a}$ with a temperature-dependent exponent
$a(T)$ shown in Fig. 4b. This is reminiscent of the BKT
Berezinskii–Kosterlitz–Thouless (BKT) transition framework, based on breaking
of vortex-antivortex pairs by DC current in a 2D superconductor. The BKT
transition point signaled by $a=3$ is at 365 mK, only slightly below $T_{c}=$
380 mK. Very similar transitions have been reported in LaAlO3/SrTiO3 2DEGs
[26, 32, 45, 34]. As pointed out in [45, 34] a failure point of the BKT
description is that it entails a discontinuous jump in $a$, rather than a
gradual transition. The rate of increase in $a$ with $T$ in Fig. 4b is similar
to that in [45, 34].

An alternate picture of this transition is in terms of mesoscopically
inhomogeneous superconductivity. A common approach for SrTiO3-based 2DEGs is
to describe the system by effective medium theory [34, 28, 29]: normal state
regions coexist with superconducting puddles with a probabilistic distribution
of $T_{c}$. For the thermally-driven resistive transition this predicts

$R(T)=R_{N}\left(1-w+w\,\mathrm{erf}\left(\frac{T-\overline{T}_{c}}{\sqrt{2}\gamma}\right)\right),$
(1)

Where $w$ is the fraction of superconducting puddles, $\overline{T}_{c}$ and
$\gamma$ are the average value and width of the Gaussian distribution for
$T_{c}$, $R_{N}$ is the normal state resistance. A fit line with $w=0.5$,
$\overline{T}_{c}=$ 380 K, $\gamma=$ 35 K is shown in Fig. 4b. Similarly to
[29, 33, 34], this simple model successfully accounts for the large width of
the resistance drop, but fails to describe the prolonged tail above $T_{c}$.
The latter shortcoming may reflect a more complex statistical distribution of
local $T_{c}$, for example a bimodal distribution including a second set of
puddles with $T_{c}\approx$ 400-700 mK [28].

Further support for the inhomogeneous superconductivity picture comes from
analysis of the superconducting critical current $I_{c}$ (Fig. 4c-e). $I_{c}$
was defined as $I_{\text{DC}}$ at which $V_{\text{DC}}$ reaches 5 $\mu$V.
Similarly to [33, 34], its temperature dependence is well described by the
Ambegaokar-Baratoff model of a dirty superconducting weak link:

$I_{c}(T)=\frac{\pi\Delta(T)}{2eR_{N}^{*}}\tanh\frac{\Delta(T)}{2k_{B}T},$ (2)

Where $\Delta(T)=\Delta(0)\tanh{(1.74\sqrt{(T_{c}-T)/T})}$ is the BCS
superconducting gap, and $R_{N}^{*}$ is the normal state weak link resistance.
For the inhmogeneous case, this model corresponds to an array of Josephson
junctions between superconducting puddles with the same functional form of
$I_{c}(T)$ and a total normal state resistance $R_{N}^{*}$ [33]. The fit in
Fig. 4c is to a normalized lineshape with $I_{c}(0)=$ 430 nA and $T_{c}=$ 380
mK.

At base temperature, the suppression of $I_{c}$ with a magnetic field $B$
orthogonal to the 2DEG plane (Fig. 4d) is well described by the random
junction array model [33, 46, 47]:

$I_{c}(B)=I_{c}(0)\cdot\frac{1}{1+(B/B_{0})^{\beta}},$ (3)

Where $B_{0}=\Phi_{0}/\pi A_{0}$, $A_{0}$ is the junction area, $\beta$ is an
array geometry-dependent exponent, typically 1-2 [33, 46, 47]. Fitting in Fig.
4e gives $\beta=2.3$, $\sqrt{A_{0}}=80$ nm, consistent with a mesoscopic array
of Josephson junctions.

Figures 4f,g show how superconductivity is suppressed with $B$ and $T$ in a
4-probe AC measurement at zero DC bias (for $V_{\text{GIL}}=$ 3.5 V). The
temperature dependence of the critical field $B_{c}$ (defined as the midpoint
of the resistance drop in $B$) presents a very rapid decrease close to
$T_{c}$. The concave curvature is stronger than both the 2D BCS dependence
($B_{c}/B_{c}(0)=1-T/T_{c}$) and 3D BCS dependence
($B_{c}/B_{c}(0)=\sqrt{1-T/T_{c}}$). A similarly strong curvature in
$B_{c}(T)$ can appear in the Werthamer-Helfand-Hohenberg (WHH) model [44, 48,
49] in the limit of strong orbital depairing ($\alpha$) and low spin-orbit
scattering strength ($\lambda_{\text{SO}}$) [49]. The data in Fig. 4g are well
matched by the normalized WHH curve with $\alpha=1.6$ and
$\lambda_{\text{SO}}=0$. Spin-orbit scattering values in STO span a wide range
as a function of carrier density and between different experiments [17]. So
negligible $\lambda_{\text{SO}}$ is not too surprising, but it is a departure
from the common pattern of high spin-orbit coupling approximately coinciding
with optimal superconducting $T_{c}$ [27, 17].

Doping-induced insulator-to-metal transition and inhomogeneous
superconductivity are key features of the carrier density-temperature phase
diagram in SrTiO3-based 2DEGs. Both are preserved in our device despite its
exposure to FIB milling. These are transport-based indications that our
approach can preserve a chosen channel surface and the relevant host material
physics.

## Conclusions

We have demonstrated a unique combination of top-down patterning with FIB,
EBL, and voltage-driven accumulation of carrier density. The key enabler is
the use of sacrificial Ge layer to protect a pristine surface from FIB damage.
This technique can be used to fabricate gated devices from materials other
than SrTiO3, for which device fabrication by standard lithography is not
feasible because large single crystal or thin films are not available.

Separately, combining ionic liquid gating and/or conventional electrostatic
gating with three dimensional shape control on the micron scale by FIB
presents exciting opportunities. For example, 2DEG functionality could be
controlled with engineered local strain [36] or piezoelectric actuation of
stretching, compression, or bending of FIB-defined lamellas [6]. Provided the
lamella material is piezoelectric or has a large electrostrictive coefficient
(as is the case for SrTiO3 [50]) bulk or surface acoustic waves could be
generated and acoustoelectrically coupled to the induced 2DEG [51]. Here, FIB
would both define key dimensions and thus frequencies of resonant acoustic
waves, and enable releasing micron-scale structures from a larger substrate.

## References

  * Moll [2018] P. J. W. Moll, Focused Ion Beam Microstructuring of Quantum Matter, Annual Review of Condensed Matter Physics 9, 147 (2018).
  * Bruchhaus  _et al._ [2017] L. Bruchhaus, P. Mazarov, L. Bischoff, J. Gierak, A. D. Wieck, and H. Hövel, Comparison of technologies for nano device prototyping with a special focus on ion beams: A review, Applied Physics Reviews 4, 011302 (2017).
  * Moll  _et al._ [2016a] P. J. Moll, P. Kushwaha, N. Nandi, B. Schmidt, and A. P. Mackenzie, Evidence for hydrodynamic electron flow in PdCoO2, Science 351, 1061 (2016a).
  * Yasui  _et al._ [2020] Y. Yasui, K. Lahabi, V. F. Becerra, R. Fermin, M. S. Anwar, S. Yonezawa, T. Terashima, M. V. Milošević, J. Aarts, and Y. Maeno, Spontaneous emergence of Josephson junctions in homogeneous rings of single-crystal Sr2RuO4, NPJ Quantum Materials 5, 21 (2020).
  * Moll  _et al._ [2016b] P. J. Moll, N. L. Nair, T. Helm, A. C. Potter, I. Kimchi, A. Vishwanath, and J. G. Analytis, Transport evidence for Fermi-arc-mediated chirality transfer in the Dirac semimetal Cd3As2, Nature 535, 266 (2016b).
  * Kim  _et al._ [2018] H.-H. Kim, S. M. Souliou, M. E. Barber, E. Lefrançois, M. Minola, M. Tortora, R. Heid, N. Nandi, R. A. Borzi, G. Garbarino, A. Bosak, J. Porras, T. Loew, M. König, P. J. W. Moll, A. P. Mackenzie, B. Keimer, C. W. Hicks, and M. Le Tacon, Uniaxial pressure control of competing orders in a high-temperature superconductor, Science 362, 1040 (2018).
  * Ronning  _et al._ [2017] F. Ronning, T. Helm, K. Shirer, M. Bachmann, L. Balicas, M. K. Chan, B. Ramshaw, R. D. Mcdonald, F. F. Balakirev, M. Jaime, E. D. Bauer, and P. J. W. Moll, Electronic in-plane symmetry breaking at field-tuned quantum criticality in CeRhIn5, Nature 548, 313 (2017).
  * Maniv  _et al._ [2021] E. Maniv, N. L. Nair, S. C. Haley, S. Doyle, C. John, S. Cabrini, A. Maniv, S. K. Ramakrishna, Y.-L. Tang, P. Ercius, R. Ramesh, Y. Tserkovnyak, A. P. Reyes, and J. G. Analytis, Antiferromagnetic switching driven by the collective dynamics of a coexisting spin glass, Science Advances 7, eabd8452 (2021).
  * Ooi  _et al._ [2002] S. Ooi, T. Mochiku, and K. Hirata, Periodic Oscillations of Josephson-Vortex Flow Resistance in Bi2Sr2CaCu2O8+y, Physical Review Letters 89, 247002 (2002).
  * Bachmann  _et al._ [2017] M. D. Bachmann, N. Nair, F. Flicker, R. Ilan, T. Meng, N. J. Ghimire, E. D. Bauer, F. Ronning, J. G. Analytis, and P. J. W. Moll, Inducing superconductivity in Weyl semimetal microstructures by selective ion sputtering, Science Advances 3, e1602983 (2017).
  * Ramshaw  _et al._ [2018] B. Ramshaw, K. A. Modic, A. Shekhter, Y. Zhang, E.-A. Kim, P. J. Moll, M. D. Bachmann, M. Chan, J. Betts, F. Balakirev, A. Migliori, N. Ghimire, E. Bauer, F. Ronning, and R. McDonald, Quantum limit transport and destruction of the Weyl nodes in TaAs, Nature Communications 9, 2217 (2018).
  * Kato  _et al._ [1999] N. I. Kato, Y. Kohno, and H. Saka, Side-wall damage in a transmission electron microscopy specimen of crystalline Si prepared by focused ion beam etching, Journal of Vacuum Science & Technology A: Vacuum, Surfaces, and Films 17, 1201 (1999).
  * Giannuzzi and Stevie [1999] L. A. Giannuzzi and F. A. Stevie, A review of focused ion beam milling techniques for TEM specimen preparation, Micron 30, 197 (1999).
  * Kelley  _et al._ [2013] R. D. Kelley, K. Song, B. Van Leer, D. Wall, and L. Kwakman, Xe+ FIB milling and measurement of amorphous silicon damage, Microscopy and Microanalysis 19, 862 (2013).
  * Ueno  _et al._ [2008] K. Ueno, S. Nakamura, H. Shimotani, A. Ohtomo, N. Kimura, T. Nojima, H. Aoki, Y. Iwasa, and M. Kawasaki, Electric-field-induced superconductivity in an insulator, Nature Materials 7, 855 (2008).
  * Sulpizio  _et al._ [2014] J. A. Sulpizio, S. Ilani, P. Irvin, and J. Levy, Nanoscale phenomena in oxide heterostructures, Annual Review of Materials Research 44, 117 (2014).
  * Pai  _et al._ [2018] Y.-Y. Pai, A. Tylan-Tyler, P. Irvin, and J. Levy, Physics of SrTiO3-based heterostructures and nanostructures: a review, Reports on Progress in Physics 81, 036503 (2018).
  * Melosh  _et al._ [2020] N. A. Melosh, M. R. Angle, M.-E. S. Hanna, and Y. Kong, Method for lithograghic patterning of sensitive materials (2020), US Patent 10,581,003.
  * Lee  _et al._ [2011a] M. Lee, J. R. Williams, S. Zhang, C. D. Frisbie, and D. Goldhaber-Gordon, Electrolyte gate-controlled kondo effect in srtio3, Physical Review Letters 107, 256601 (2011a).
  * Lee  _et al._ [2011b] Y. Lee, C. Clement, J. Hellerstedt, J. Kinney, L. Kinnischtzke, X. Leng, S. D. Snyder, and A. M. Goldman, Phase diagram of electrostatically doped SrTiO3, Physical Review Letters 106, 136809 (2011b).
  * Li  _et al._ [2012] M. Li, T. Graf, T. D. Schladt, X. Jiang, and S. S. P. Parkin, Role of percolation in the conductance of electrolyte-gated SrTiO3, Physical Review Letters 109, 196803 (2012).
  * Stanwyck  _et al._ [2013] S. W. Stanwyck, P. Gallagher, J. R. Williams, and D. Goldhaber-Gordon, Universal conductance fluctuations in electrolyte-gated SrTiO3 nanostructures, Applied Physics Letters 103, 213504 (2013).
  * Gallagher  _et al._ [2014] P. Gallagher, M. Lee, J. R. Williams, and D. Goldhaber-Gordon, Gate-tunable superconducting weak link and quantum point contact spectroscopy on a strontium titanate surface, Nature Physics 10, 748 (2014).
  * Gallagher  _et al._ [2015] P. Gallagher, M. Lee, T. A. Petach, S. W. Stanwyck, J. R. Williams, K. Watanabe, T. Taniguchi, and D. Goldhaber-Gordon, A high-mobility electronic system at an electrolyte-gated oxide surface, Nature Communications 6, 6437 (2015).
  * Mikheev  _et al._ [2021] E. Mikheev, I. T. Rosen, and D. Goldhaber-Gordon, Quantized critical supercurrent in SrTiO3-based quantum point contacts, Science Advances 7, eabi6520 (2021).
  * Reyren  _et al._ [2007] N. Reyren, S. Thiel, A. D. Caviglia, L. F. Kourkoutis, G. Hammerl, C. Richter, C. W. Schneider, T. Kopp, A.-S. Rüetschi, D. Jaccard, M. Gabay, D. A. Muller, J.-M. Triscone, and J. Mannhart, Superconducting interfaces between insulating oxides, Science 317, 1196 (2007).
  * Caviglia  _et al._ [2010] A. Caviglia, S. Gariglio, C. Cancellieri, B. Sacépé, A. Fete, N. Reyren, M. Gabay, A. F. Morpurgo, and J.-M. Triscone, Two-dimensional quantum oscillations of the conductance at LaAlO3/SrTiO3 interfaces, Physical Review Letters 105, 236802 (2010).
  * Caprara  _et al._ [2011] S. Caprara, M. Grilli, L. Benfatto, and C. Castellani, Effective medium theory for superconducting layers: A systematic analysis including space correlation effects, Physical Review B 84, 014514 (2011).
  * Caprara, S and Biscaras, J and Bergeal, N and Bucheli, D and Hurand, S and Feuillet-Palma, C and Rastogi, A and Budhani, RC and Lesueur, J and Grilli, M [2013] Caprara, S and Biscaras, J and Bergeal, N and Bucheli, D and Hurand, S and Feuillet-Palma, C and Rastogi, A and Budhani, RC and Lesueur, J and Grilli, M, Multiband superconductivity and nanoscale inhomogeneity at oxide interfaces, Physical Review B 88, 020504 (2013).
  * Seri  _et al._ [2013] S. Seri, M. Schultz, and L. Klein, Thermally activated recovery of electrical conductivity in LaAlO3/SrTiO3, Physical Review B 87, 125110 (2013).
  * Biscaras  _et al._ [2014] J. Biscaras, S. Hurand, C. Feuillet-Palma, A. Rastogi, R. C. Budhani, N. Reyren, E. Lesne, J. Lesueur, and N. Bergeal, Limit of the electrostatic doping in two-dimensional electron gases of La$X$O3($X$ = Al, Ti)/SrTiO3, Scientific Reports 4, 6788 (2014).
  * Han  _et al._ [2014] Y.-L. Han, S.-C. Shen, J. You, H.-O. Li, Z.-Z. Luo, C.-J. Li, G.-L. Qu, C.-M. Xiong, R.-F. Dou, L. He, D. Naugle, G.-P. Guo, and J.-C. Nie, Two-dimensional superconductivity at (110) LaAlO3/SrTiO3 interfaces, Applied Physics Letters 105, 192603 (2014).
  * Prawiroatmodjo  _et al._ [2016] G. E. D. K. Prawiroatmodjo, F. Trier, D. V. Christensen, Y. Chen, N. Pryds, and T. S. Jespersen, Evidence of weak superconductivity at the room-temperature grown LaAlO3/SrTiO3 interface, Physical Review B 93, 184504 (2016).
  * Venditti  _et al._ [2019] G. Venditti, J. Biscaras, S. Hurand, N. Bergeal, J. Lesueur, A. Dogra, R. Budhani, M. Mondal, J. Jesudasan, P. Raychaudhuri, S. Caprara, and L. Benfatto, Nonlinear $I\text{-}V$ characteristics of two-dimensional superconductors: Berezinskii-Kosterlitz-Thouless physics versus inhomogeneity, Physical Review B 100, 064506 (2019).
  * Yin  _et al._ [2020] C. Yin, A. E. M. Smink, I. Leermakers, L. M. K. Tang, N. Lebedev, U. Zeitler, W. G. Van Der Wiel, H. Hilgenkamp, and J. Aarts, Electron trapping mechanism in LaAlO3/SrTiO3 heterostructures, Physical Review Letters 124, 017702 (2020).
  * Bachmann  _et al._ [2019] M. D. Bachmann, G. M. Ferguson, F. Theuss, T. Meng, C. Putzke, T. Helm, K. Shirer, Y.-S. Li, K. A. Modic, M. Nicklas, M. Konig, D. Low, S. Ghosh, A. P. Mackenzie, F. Arnold, E. Hassinger, R. D. McDonald, L. E. Winter, E. D. Bauer, F. Ronning, B. J. Ramshaw, K. C. Nowack, and P. J. W. Moll, Spatial control of heavy-fermion superconductivity in CeIrIn5, Science 366, 221 (2019).
  * Gray Jr  _et al._ [2016] D. J. Gray Jr, T. A. Merz, Y. Hikita, H. Y. Hwang, and H. Mabuchi, Orientation-resolved domain mapping in tetragonal SrTiO3 using polarized Raman spectroscopy, Physical Review B 94, 214107 (2016).
  * Kalisky  _et al._ [2013] B. Kalisky, E. M. Spanton, H. Noad, J. R. Kirtley, K. C. Nowack, C. Bell, H. K. Sato, M. Hosoda, Y. Xie, Y. Hikita, C. Woltmann, G. Pfanzelt, R. Jany, C. Richter, H. Y. Hwang, J. Mannhart, and K. A. Moler, Locally enhanced conductivity due to the tetragonal domain structure in LaAlO3/SrTiO3 heterointerfaces, Nature Materials 12, 1091 (2013).
  * Noad  _et al._ [2016] H. Noad, E. M. Spanton, K. C. Nowack, H. Inoue, M. Kim, T. A. Merz, C. Bell, Y. Hikita, R. Xu, W. Liu, A. Vailionis, H. Y. Hwang, and K. A. Moler, Variation in superconducting transition temperature due to tetragonal domains in two-dimensionally doped SrTiO3, Physical Review B 94, 174516 (2016).
  * Noad  _et al._ [2019] H. Noad, P. Wittlich, J. Mannhart, and K. A. Moler, Modulation of Superconducting Transition Temperature in LaAlO3/SrTiO3 by SrTiO3 Structural Domains, Journal of Superconductivity and Novel Magnetism 32, 821 (2019).
  * Ahadi  _et al._ [2019] K. Ahadi, L. Galletti, Y. Li, S. Salmani-Rezaie, W. Wu, and S. Stemmer, Enhancing superconductivity in SrTiO3 films with strain, Science Advances 5, eaaw0120 (2019).
  * Herrera  _et al._ [2019] C. Herrera, J. Cerbin, A. Jayakody, K. Dunnett, A. V. Balatsky, and I. Sochnikov, Strain-engineered interaction of quantum polar and superconducting phases, Physical Review Materials 3, 124801 (2019).
  * Enderlein  _et al._ [2020] C. Enderlein, J. F. de Oliveira, D. A. Tompsett, E. B. Saitovitch, S. S. Saxena, G. G. Lonzarich, and S. E. Rowley, Superconductivity mediated by polar modes in ferroelectric metals, Nature Communications 11, 4852 (2020).
  * Werthamer  _et al._ [1966] N. R. Werthamer, E. F. Helfand, and P. C. Hohenberg, Temperature and purity dependence of the superconducting critical field, $H_{c2}$. III. Electron spin and spin-orbit effects, Physical Review 147, 295 (1966).
  * Monteiro  _et al._ [2017] A. M. R. V. L. Monteiro, D. J. Groenendijk, I. Groen, J. de Bruijckere, R. Gaudenzi, H. S. J. Van Der Zant, and A. D. Caviglia, Two-dimensional superconductivity at the (111) LaAlO3/SrTiO3 interface, Physical Review B 96, 020504 (2017).
  * Müller  _et al._ [1992] K.-H. Müller, D. N. Matthews, and R. Driver, Critical current density of ceramic high-temperature superconductors in a low magnetic field, Physica C: Superconductivity 191, 339 (1992).
  * Van der Laan  _et al._ [2001] D. C. Van der Laan, H. J. N. Van Eck, B. ten Haken, J. Schwartz, and H. H. J. ten Kate, Temperature and magnetic field dependence of the critical current of Bi2Sr2Ca2Cu3Ox tape conductors, IEEE transactions on applied superconductivity 11, 3345 (2001).
  * Kim  _et al._ [2012] M. Kim, Y. Kozuka, C. Bell, Y. Hikita, and H. Y. Hwang, Intrinsic spin-orbit coupling in superconducting $\delta$-doped SrTiO3 heterostructures, Physical Review B 86, 085121 (2012).
  * Solenov  _et al._ [2017] D. Solenov, M. Nikolo, J. Singleton, J. Jiang, J. Weiss, and E. Hellstrom, The impact of Maki parameter and spin orbit scattering constant in the WHH model on upper critical magnetic fields in Ni-and Co-doped pnictide bulk superconductors, in  _AIP Conference Proceedings_ , Vol. 1895 (2017) p. 060004.
  * Grupp and Goldman [1997] D. E. Grupp and A. M. Goldman, Giant piezoelectric effect in strontium titanate at cryogenic temperatures, Science 276, 392 (1997).
  * Uzun  _et al._ [2020] Y. Uzun, A. E. M. Smink, M. P. De Jong, H. Hilgenkamp, and W. G. Van Der Wiel, Acoustoelectric charge transport at the LaAlO3/SrTiO3 interface, Applied Physics Letters 116, 011601 (2020).

## Acknowledgments

We thank Carsten Putzke for help with focused ion beam work; Trevor Petach,
Xiao Zhang, and Joe Finney for help with the development of lithographic
fabrication processes; and Ilan Rosen and Eli Fox for help with cryogenic
measurements. Our use of Germanium as a water-soluble protective layer was
inspired by discussions with Nicholas Melosh, prior to publication of his work
on Germanium as a protective layer for lithography on organic materials.

Funding: Experimental work (fabrication and measurement) by E.M. was primarily
supported by the Air Force Office of Scientific Research through grant no.
FA9550-16-1-0126. E. M. was also supported by the Nano- and Quantum Science
and Engineering Postdoctoral Fellowship at Stanford University and (during
paper-writing and analysis) by the Department of Energy, Office of Science,
Basic Energy Sciences, Materials Sciences and Engineering Division, under
Contract DE-AC02-76SF00515. Measurement infrastructure was funded in part by
the Gordon and Betty Moore Foundation’s EPiQS Initiative through grant
GBMF3429 and grant GBMF9460. Part of this work was performed at the Stanford
Nano Shared Facilities (SNSF), supported by the National Science Foundation
under award ECCS-1542152.

Author contributions: E.M., P.J.W.M., T.Z. and D.G.-G. designed the
experiment. T.Z., A.E. and P.J.W.M. carried out FIB-based device fabrication.
E.M. carried out EBL-based device fabrication and cryogenic measurements. All
authors discussed the results. E.M. wrote the manuscript with input from other
authors.


# Differentiability of the $n$-Variable Function Deduced by the
Differentiability of the $n-1$-Variable Function

Zhenglin Ye, Qianqiao Guo  Zhenglin Ye, School of Mathematics and Statistics,
Northwestern Polytechnical University, Xi’an, 710129, Shaanxi, China
<EMAIL_ADDRESS>Qianqiao Guo, School of Mathematics and Statistics,
Northwestern Polytechnical University, Xi’an, 710129, Shaanxi, China
<EMAIL_ADDRESS>
###### Abstract.

In this paper, some sufficient conditions for the differentiability of the
$n$-variable real-valued function are obtained, which are given based on the
differentiability of the $n-1$-variable real-valued function and are weaker
than classical conditions.

###### Key words and phrases:

Sufficient Condition, Differentiability, $n$-Variable Function

###### 2010 Mathematics Subject Classification:

26B05

##  1\. Introduction

Consider a real-valued function $f:U(a)\to\mathbb{R}$ defined in a
neighbourhood $U(a)\subset\mathbb{R}^{n},n\geq 2$ of the point
$a:=(a_{1},a_{2},\cdots,a_{n})$.

A classical sufficient condition for the differentiability of $f$ at the point
$a$ states as (see for example [1]):

######  Proposition 1.1.

If the function $f$ has all partial derivatives
$f^{\prime}_{x_{1}},\cdots,f^{\prime}_{x_{n}}$ at each point of the
neighbourhood $U(a)$ and they are continuous at $a$, then $f$ is
differentiable at $a$.

In this paper, we show a different sufficient condition for the
differentiability of $f(x_{1},x_{2},\cdots,x_{n})$ at the point $a$, which are
weaker than that of Proposition 1.1.

The main result is stated below.

######  Theorem 1.2.

Assume the function $f$ has all partial derivatives at the point $a$, and has
partial derivatives
$f^{\prime}_{x_{i_{0}}}(x_{1},\cdots,x_{i_{0}-1},a_{i_{0}},x_{i_{0}+1},\cdots,x_{n})$
at some neighbourhood of the point
$\widehat{a}_{i_{0}}:=(a_{1},\cdots,a_{i_{0}-1},\wedge,a_{i_{0}+1},\cdots,a_{n})\in\mathbb{R}^{n-1}$,
for fixed $1\leq i_{0}\leq n$. Then the function $f$ is differentiable at $a$
if and only if

(i) the $n-1$-variable function
$f(x_{1},\cdots,x_{i_{0}-1},a_{i_{0}},x_{{i_{0}}+1},\cdots,x_{n})$ is
differentiable at the point $\widehat{a}_{i_{0}}$;

(ii) the $n-1$-variable function
$f^{\prime}_{x_{i_{0}}}(x_{1},\cdots,x_{{i_{0}}-1},a_{i_{0}},x_{{i_{0}}+1},\cdots,x_{n})$
is continuous at the point $\widehat{a}_{i_{0}}$.

######  Remark 1.3.

(1) The interesting point of Theorem 1.2 is that the differentiability of
$f(x_{1},x_{2},\cdots,x_{n})$ at the point $a$ can be deduced by the
differentiability of the $n-1$-variable real-valued function
$f(x_{1},x_{2},\cdots,x_{i_{0}-1},a_{i_{0}},x_{i_{0}+1},\cdots,x_{n})(1\leq
i_{0}\leq n)$ with some additional assumptions.

(2) Notice that if $f^{\prime}_{x_{k}}(x_{1},\cdots,x_{n})$ is continuous at
$a$, then the conditions (i) and (ii) hold. But we cannot deduce the
continuity of $f^{\prime}_{x_{k}}(x_{1},\cdots,x_{n})$ at $a$ for any $1\leq
k\leq n$ by (i) and (ii). Thus, the conditions of Theorem 1.2 are weaker than
that of Proposition 1.1.

By inducement, it is easy to see that the following holds.

######  Theorem 1.4.

Assume the function $f$ has all partial derivatives at the point $a$, and has
partial derivatives
$f^{\prime}_{x_{k}}(a_{1},\cdots,a_{k},x_{k+1},\cdots,x_{n})$ at some
neighbourhood of the point
$\widehat{a}^{n-k}:=(a_{k+1},\cdots,a_{n})\in\mathbb{R}^{n-k}(k=1,\cdots,n-1)$.
Then the function $f$ is differentiable at $a$ if and only if

(1) the $n-1$-variable function $f^{\prime}_{x_{1}}(a_{1},x_{2},\cdots,x_{n})$
is continuous at the point $\widehat{a}^{n-1}$;

(2) the $n-2$-variable function
$f^{\prime}_{x_{2}}(a_{1},a_{2},x_{3},\cdots,x_{n})$ is continuous at the
point $\widehat{a}^{n-2}$;

$\cdots$

(n-1) the $1$-variable function
$f^{\prime}_{x_{n-1}}(a_{1},\cdots,a_{n-1},x_{n})$ is continuous at the point
$\widehat{a}^{n-(n-1)}=a_{n}$.

As an example, the case $n=3$ is stated as follows.

######  Corollary 1.5.

Assume the function $f(x_{1},x_{2},x_{3})$ has all partial derivatives at the
point $a$, and has partial derivative $f^{\prime}_{x_{1}}(a_{1},x_{2},x_{3})$
at some neighbourhood of the point $(a_{2},a_{3})$, and has partial derivative
$f^{\prime}_{x_{2}}(a_{1},a_{2},x_{3})$ at some neighbourhood of the point
$a_{3}$. Then the function $f$ is differentiable at $a$ if and only if
$f^{\prime}_{x_{1}}(a_{1},x_{2},x_{3})$ is continuous at the point
$(a_{2},a_{3})$, and $f^{\prime}_{x_{2}}(a_{1},a_{2},x_{3})$ is continuous at
the point $a_{3}$.

######  Remark 1.6.

The functions $f^{\prime}_{x_{1}}(a_{1},x_{2},x_{3})$ and
$f^{\prime}_{x_{2}}(a_{1},a_{2},x_{3})$ above can also be replaced by, for
example, $f^{\prime}_{x_{2}}(x_{1},a_{2},x_{3})$ and
$f^{\prime}_{x_{3}}(x_{1},a_{2},a_{3})$, or
$f^{\prime}_{x_{3}}(x_{1},x_{2},a_{3})$ and
$f^{\prime}_{x_{1}}(a_{1},x_{2},a_{3})$.

##  2\. Proofs of the main results

###### Proof of Theorem 1.2.

(i@) The “if” part.

Denote

$\displaystyle\Delta R_{n}$ $\displaystyle:=$ $\displaystyle f(a_{1}+\Delta
x_{1},a_{2}+\Delta x_{2},\cdots,a_{n}+\Delta
x_{n})-f(a_{1},a_{2},\cdots,a_{n})-\sum\limits_{k=1}^{n}f^{\prime}_{x_{k}}(a_{1},a_{2},\cdots,a_{n})\Delta
x_{k},$ $\displaystyle\Delta R_{n-1}^{i_{0}}$ $\displaystyle:=$ $\displaystyle
f(a_{1}+\Delta x_{1},\cdots,a_{{i_{0}}-1}+\Delta
x_{{i_{0}}-1},a_{i_{0}},a_{{i_{0}}+1}+\Delta x_{{i_{0}}+1},\cdots,a_{n}+\Delta
x_{n})-f(a_{1},a_{2},\cdots,a_{n})$
$\displaystyle-\sum\limits_{k=1,k\neq{i_{0}}}^{n}f^{\prime}_{x_{k}}(a_{1},a_{2},\cdots,a_{n})\Delta
x_{k}$

and

$\rho_{n}:=\Bigg{(}\sum\limits_{k=1}^{n}(\Delta
x_{k})^{2}\Bigg{)}^{\frac{1}{2}},\ \
\rho_{n-1}^{i_{0}}:=\Bigg{(}\sum\limits_{k=1,k\neq{i_{0}}}^{n}(\Delta
x_{k})^{2}\Bigg{)}^{\frac{1}{2}}.$

Since $f$ is differentiable at $a$, we have

$\lim\limits_{\rho_{n}\to 0}\frac{\Delta R_{n}}{\rho_{n}}=0.$

Let $\Delta x_{i_{0}}=0$. Then it is easy to see that

$\lim\limits_{\rho_{n-1}^{i_{0}}\to 0}\frac{\Delta
R_{n-1}^{i_{0}}}{\rho_{n-1}^{i_{0}}}=0,$ (2.3)

which implies that the $n-1$-variable function
$f(x_{1},\cdots,x_{{i_{0}}-1},a_{i_{0}},x_{{i_{0}}+1},\cdots,x_{n})$ is
differentiable at the point $\widehat{a}_{i_{0}}$.

On the other hand, denote

$\displaystyle\Delta
f^{\prime}_{x_{i_{0}}}(a_{1},\cdots,a_{{i_{0}}-1},\wedge,a_{{i_{0}}+1},\cdots,a_{n})$
$\displaystyle:=$ $\displaystyle f^{\prime}_{x_{i_{0}}}(a_{1}+\Delta
x_{1},\cdots,a_{{i_{0}}-1}+\Delta x_{{i_{0}}-1},a_{i_{0}},a_{{i_{0}}+1}+\Delta
x_{{i_{0}}+1}\cdots,a_{n}+\Delta
x_{n})-f^{\prime}_{x_{i_{0}}}(a_{1},\cdots,a_{n}).$

Since $f$ has all derivatives at $a$, by using (2) and (2), we have

$\displaystyle\Delta R_{n}$ (2.4) $\displaystyle=$ $\displaystyle
f(a_{1}+\Delta x_{1},\cdots,a_{{i_{0}}-1}+\Delta
x_{{i_{0}}-1},a_{i_{0}},a_{{i_{0}}+1}+\Delta x_{{i_{0}}+1},\cdots,a_{n}+\Delta
x_{n})$ $\displaystyle+f^{\prime}_{x_{i_{0}}}(a_{1}+\Delta
x_{1},\cdots,a_{{i_{0}}-1}+\Delta x_{{i_{0}}-1},a_{i_{0}},a_{{i_{0}}+1}+\Delta
x_{{i_{0}}+1},\cdots,a_{n}+\Delta x_{n})\Delta x_{i_{0}}+\alpha_{i_{0}}\Delta
x_{i_{0}}$
$\displaystyle-f(a_{1},a_{2},\cdots,a_{n})-\sum\limits_{k=1}^{n}f^{\prime}_{x_{k}}(a_{1},a_{2},\cdots,a_{n})\Delta
x_{k},$ $\displaystyle=$ $\displaystyle\Delta R_{n-1}^{i_{0}}+\Delta
f^{\prime}_{x_{i_{0}}}(a_{1},\cdots,a_{{i_{0}}-1},\wedge,a_{{i_{0}}+1},\cdots,a_{n})\Delta
x_{i_{0}}+\alpha_{i_{0}}\Delta x_{i_{0}},$

where $\alpha_{i_{0}}\to 0$ as $\Delta x_{i_{0}}\to 0$. Take $\Delta
x_{i_{0}}=\rho_{n-1}^{i_{0}}$. Then
$\rho_{n}=2^{\frac{1}{2}}\rho_{n-1}^{i_{0}}$. Noticing that $f$ is
differentiable at $a$, by using (2.3) and (2.4) we have

$\displaystyle 0$ $\displaystyle=$ $\displaystyle\lim\limits_{\rho_{n}\to
0}\frac{\Delta R_{n}}{\rho_{n}}$ $\displaystyle=$
$\displaystyle\lim\limits_{\rho_{n}\to 0}\Bigg{(}\frac{\Delta
R_{n-1}^{i_{0}}}{\rho_{n}}+\Delta
f^{\prime}_{x_{i_{0}}}(a_{1},\cdots,a_{{i_{0}}-1},\wedge,a_{{i_{0}}+1},\cdots,a_{n})\frac{\Delta
x_{i_{0}}}{\rho_{n}}+\alpha_{i_{0}}\frac{\Delta x_{i_{0}}}{\rho_{n}}\Bigg{)}$
$\displaystyle=$ $\displaystyle
2^{-\frac{1}{2}}\lim\limits_{\rho_{n-1}^{i_{0}}\to 0}\Delta
f^{\prime}_{x_{i_{0}}}(a_{1},\cdots,a_{{i_{0}}-1},\wedge,a_{{i_{0}}+1},\cdots,a_{n}),$

which implies that the partial derivative
$f^{\prime}_{x_{i_{0}}}(x_{1},\cdots,x_{{i_{0}}-1},a_{i_{0}},x_{{i_{0}}+1},\cdots,x_{n})$
is continuous at the point $\widehat{a}_{i_{0}}$.

(ii@) The “only if” part.

Again by using (2.4), we have

$\displaystyle\lim\limits_{\rho_{n}\to 0}\frac{\Delta R_{n}}{\rho_{n}}$
$\displaystyle=$ $\displaystyle\lim\limits_{\rho_{n}\to 0}\Bigg{(}\frac{\Delta
R_{n-1}^{i_{0}}}{\rho_{n}}+\Delta
f^{\prime}_{x_{i_{0}}}(a_{1},\cdots,a_{{i_{0}}-1},\wedge,a_{{i_{0}}+1},\cdots,a_{n})\frac{\Delta
x_{i_{0}}}{\rho_{n}}+\alpha_{i_{0}}\frac{\Delta x_{i_{0}}}{\rho_{n}}\Bigg{)},$

where $\alpha_{i_{0}}\to 0$ as $\Delta x_{i_{0}}\to 0$. Since the
$n-1$-variable function
$f(x_{1},\cdots,x_{{i_{0}}-1},a_{i_{0}},x_{{i_{0}}+1},\cdots,x_{n})$ is
differentiable at the point $\widehat{a}_{i_{0}}$, then (2.3) holds. Thus

$\lim\limits_{\rho_{n}\to 0}\frac{\Delta
R_{n-1}^{i_{0}}}{\rho_{n}}=\lim\limits_{\rho_{n-1}^{i_{0}}\to 0}\frac{\Delta
R_{n-1}^{i_{0}}}{\rho_{n-1}^{i_{0}}}\cdot\frac{\rho_{n-1}^{i_{0}}}{\rho_{n}}=0.$
(2.6)

On the other hand, the partial derivative
$f^{\prime}_{x_{i_{0}}}(x_{1},\cdots,x_{{i_{0}}-1},a_{i_{0}},x_{{i_{0}}+1},\cdots,x_{n})$
is continuous at the point $\widehat{a}_{i_{0}}$. Then

$\lim\limits_{\rho_{n-1}^{i_{0}}\to 0}\Delta
f^{\prime}_{x_{i_{0}}}(a_{1},\cdots,a_{{i_{0}}-1},\wedge,a_{{i_{0}}+1},\cdots,a_{n})=0.$
(2.7)

Therefore by using (2), (2.6) and (2.7), we conclude that

$\lim\limits_{\rho_{n}\to 0}\frac{\Delta R_{n}}{\rho_{n}}=0.$

That is, the function $f$ is differentiable at $a$. ∎

###### Proof of Theorem 1.4.

By using Theorem 1.2, it is easy to end the proof. We omit it here. ∎

## References

  * [1] Vladimir A. Zorich, Mathematical Analysis I, Springer-Verlag, Berlin, 2004, 219-232.