# Some properties of the torsion function with Robin boundary conditions

Rossano Sannipoli  Dipartimento di Matematica e Applicazioni “R. Caccioppoli”,
Università degli studi di Napoli Federico II
Via Cintia, Complesso Universitario Monte S. Angelo, 80126 Napoli, Italy.
<EMAIL_ADDRESS>
###### Abstract

In this paper we study some properties of the torsion function with Robin
boundary conditions. Here we write the shape derivative of the $L^{\infty}$
and $L^{p}$ norms, for $p\geq 1$, of the torsion function, seen as a
functional on a bounded simply connected open set
$\Omega\subset\mathbb{R}^{n}$, and prove that the balls are critical shapes
for these functionals, when the volume of $\Omega$ is preserved.
MSC 2010: 35J05 - 35J15 - 35J20 - 35J25
Keywords: Torsion problem - Robin boundary conditions - Shape derivative

AbstractAbstract Partial Differential Equations

##  1 Introduction

Let $\Omega\subset\mathbb{R}^{n}$, $n\geq 2$, be a bounded $C^{2,\alpha}$ and
simply connected open set. Let us consider the following torsion problem with
Robin boundary conditions

$\begin{cases}-\Delta u=1&\mbox{in}\ \Omega\vspace{0.2cm}\\\ \dfrac{\partial
u}{\partial\nu}+\beta u=0&\mbox{on}\ \partial\Omega\vspace{0.2cm},\end{cases}$
(1.1)

where $\nu$ stands for the outer unit normal to $\partial\Omega$ and $\beta>0$
is a positive real number, known as Robin boundary parameter.
A weak solution to (1.1) is a function $u\in H^{1}(\Omega)$ which satisfies

$\int_{\Omega}\nabla
u\nabla\phi\,dx+\beta\int_{\partial\Omega}u\phi\,d\mathcal{H}^{n-1}=\int_{\Omega}\phi,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\forall\phi\in
H^{1}(\Omega).$

When $\beta\to\infty$ we have the well known torsion problem of elasticity or
Saint Venant problem, that has been studied by many authors: different
estimates and qualitative properties have been studied, for example, by Pólya
and Szegö in [PS], or Payne in [P].
Other problems with non-linear operators and a different geoemtry of the
domain can be found. One example is given in [PPT]: if $\Omega$ is a convex
set with a hole, the authors prove that the spherical shell is a minimizer for
the torsional rigidity of the $p$-Laplacian, with outer Robin and internal
Neumann boundary conditions, when the volume and the outer perimeter are
fixed.
It is well known that the solution to problem (1.1) is unique and positive
whenever $\partial\Omega$ is sufficiently smooth.
A comparison result à la Talenti has been proved in [ANT]. They consider the
symmetrized problem, that is

$\begin{cases}-\Delta v=1&\mbox{in}\ \Omega^{\sharp}\vspace{0.2cm}\\\
\dfrac{\partial v}{\partial\nu}+\beta v=0&\mbox{on}\
\partial\Omega^{\sharp}\vspace{0.2cm},\end{cases}$ (1.2)

where $\Omega^{\sharp}$ is the ball centered in the origin having the same
measure as $\Omega$. In dimension $2$, they proved that

$u^{\sharp}(x)\leq v(x)\;\;\;\;\;\;\;x\in\Omega^{\sharp},$

where $v$ is the solution to the problem (1.2) and $u^{\sharp}$ is the Schwarz
symmetrization of $u$ (for more details see [K]). Moreover they proved, for
$n\geq 3$, that

$\|u\|_{L^{p,1}(\Omega)}\leq\|v\|_{L^{p,1}(\Omega^{\sharp})}$ (1.3)

and

$\|u\|_{L^{2p,2}(\Omega)}\leq\|v\|_{L^{2p,2}(\Omega^{\sharp})},$

for all $0<p\leq\frac{n}{n-2}$. Here $L^{p,q}(\Omega)$ is the Lorentz space
(for more detail see [L]). It follows that in dimensions greater than 3

$\|u\|_{L^{p}(\Omega)}\leq\|v\|_{L^{p}(\Omega^{\sharp})}$

when $p=1,2$.
In [BG] the authors proved, with different arguments, that
$\|u\|_{L^{1}(\Omega)}\leq\|v\|_{L^{1}(\Omega^{\sharp})}$, where $u$ is the
solution to (1.1) and $v$ is the solution to the symmetrized problem (1.2).

Problem (1.1) is a particular case of the following

$\begin{cases}-\Delta u=f&\mbox{in}\ \Omega\vspace{0.2cm}\\\ \dfrac{\partial
u}{\partial\nu}+\beta u=0&\mbox{on}\ \partial\Omega\vspace{0.2cm},\end{cases}$
(1.4)

where $f\in L^{2}(\Omega)$ and its symmetrized problem is

$\begin{cases}-\Delta v=f^{\sharp}&\mbox{in}\ \Omega^{\sharp}\vspace{0.2cm}\\\
\dfrac{\partial v}{\partial\nu}+\beta v=0&\mbox{on}\
\partial\Omega^{\sharp}\vspace{0.2cm},\end{cases}$ (1.5)

where $f^{\sharp}$ is the Schwarz symmetrization of $f$. When in (1.4) and
(1.5) we have Dirichlet boundary conditions, Talenti (see [T]), via the
rearrangement of a function and the Schwarz symmetrization, showed the
following

$u^{\sharp}(x)\leq v(x)\;\;\;\;\;x\in\Omega^{\sharp}$ (1.6)

Still in [ANT], a Talenti comparison result for the problem (1.4) has been
proved: when $n\geq 2$

$\|u\|_{L^{p,1}(\Omega)}\leq\|v\|_{L^{p,1}(\Omega^{\sharp})}\,\,\,\,\,\,\,\text{for
all}\,\,\,0<p\leq\frac{n}{2n-2},$

and

$\|u\|_{L^{2p,2}(\Omega)}\leq\|v\|_{L^{2p,2}(\Omega^{\sharp})}\,\,\,\,\,\,\,\text{for
all}\,\,\,0<p\leq\frac{n}{3n-4},$

where $v$ is solution to (1.5). So, in dimension 2, we have that

$\|u\|_{L^{p}(\Omega)}\leq\|v\|_{L^{p}(\Omega^{\sharp})},$ (1.7)

when $p=1,2$. One may ask if (1.7) is still true for larger values of $p$ in
dimension 2 or if it is valid in every dimension and value of $p$. The
authors, though, found counterexamples of the untruthfulness of these
questions when $n=2$ and $p=\infty$, and when $n=3$ and $p=2$. This led to the
following open problems:

  * •

$u^{\sharp}\leq v$ in $\Omega^{\sharp}$ for $n\geq 3$ and $f\equiv 1$ ;

  * •

$\|u\|_{L^{1}(\Omega)}\leq\|v\|_{L^{1}(\Omega^{\sharp})}$ for $n\geq 3$ and
$f\in L^{2}(\Omega)$.

In this paper we move the first steps in these directions.

In particular we set

$M(\Omega)=\|u\|_{L^{\infty}(\Omega)},$

and for every $p\in[1,+\infty)$ we denote the following functional

$F_{p}(\Omega)=\int_{\Omega}|u(x)|^{p}\,dx=\int_{\Omega}u^{p}(x)\,dx=\|u\|_{L^{p}(\Omega)}^{p},$

where $u$ is solution to (1.1). We are interested in computing the shape
derivative (see [HP]) of these two functionals and prove that the ball
centered at the origin is a critical shape for them.
Namely, if $\Omega\subset\mathbb{R}^{n}$ is a bounded $C^{2,\alpha}$ simply
connected open set, let us consider a first order perturbation

$\Omega_{t}=(\mathbbm{1}_{\mathbb{R}^{n}}+tV)(\Omega),$

with $\mathbbm{1}_{\mathbb{R}^{n}}$ being the identity function, $V$ a
$C^{2,\alpha}(\mathbb{R}^{n},\mathbb{R}^{n})$ vector field and $t$ a small
real number. We are interested in the study of the limits (if they exist)

$M^{\prime}(\Omega,V)=\lim_{t\to 0}\frac{M(\Omega_{t})-M(\Omega)}{t}$ (1.8)

and

$F_{p}^{\prime}(\Omega,V)=\lim_{t\to
0}\frac{F_{p}(\Omega_{t})-F_{p}(\Omega)}{t}.$ (1.9)

Let us denote by $B_{R}$ the ball centered at the origin in $\mathbb{R}^{n}$
with radius $R>0$. Then problem (1.1) becomes

$\begin{cases}-\Delta u=1&\mbox{in}\ B_{R}\vspace{0.2cm}\\\ \dfrac{\partial
u}{\partial\nu}+\beta u=0&\mbox{on}\ \partial B_{R}\vspace{0.2cm},\end{cases}$
(1.10)

where $\nu=\frac{x}{R}$ is the outer unit normal to the boundary.
We are now able to state the main theorem.

Main result. The ball $B_{R}$ is a critical shape for the functionals
$M(\Omega)$ and $F_{p}(\Omega)$, $p\geq 1$, i.e.

$M^{\prime}(B_{R},V)=F^{\prime}_{p}(B_{R},V)=0,$

where $V$ is a $C^{2,\alpha}(\mathbb{R}^{n},\mathbb{R}^{n})$ vector field
volume preserving of the first order and where $M^{\prime}(\cdot,v)$ and
$F_{p}^{\prime}(\cdot,v)$ are the shape derivatives of $M$ and $F_{p}$
respectively.

For the precise definition of vector field volume preserving of the first
order see section 2.

Next section will include some preliminary results in order to compute the
shape derivative of the two functionals defined before. In section $3$ there
will be the computation of the $L^{\infty}$-norm of the torsion function with
the consequent proof of the main result. In section $4$, the same thing for
the $L^{p}$-norm will be done.

##  2 Preliminaries

Throughout this paper we will denote by
$B_{R}=\\{x\in\mathbb{R}^{n}:\|x\|<R\\}$ the ball centered at the origin with
radius $R>0$, where $\|\cdot\|$ is the classical euclidean distance; by
$\Omega$ a bounded $C^{2,\alpha}$ and simply connected open set with finite
Lebeasgue measure, where $C^{2,\alpha}$ stands for the $\alpha$-Hölderian
space. We denote by $\mathcal{H}^{n-1}$ the $(n-1)$-dimensional Hausdorff
measure in $\mathbb{R}^{n}$ and by $|\cdot|$ the Lebeasgue measure in
$\mathbb{R}^{n}$. If $x,y\in\mathbb{R}^{n}$ we will indicate by $x\cdot y$ the
scalar product between the two vectors.

###  2.1 Shape derivative: some definitions and computations

Here we give some preliminary definitions and results that the reader can find
in [BW] and [HP]. We point out that in this and next subsection, we will use
the Einstein summation convention for the repeated indexes.
Let $\Omega\subset\mathbb{R}^{n}$ be a bounded and simply connected open set.
We consider a family of perturbations $\\{\Omega_{t}\\}_{t}$ of the form

$\Omega_{t}=\displaystyle\\{y=x+tV(x):x\in\Omega,t\,\,\text{small enough}\\},$
(2.1)

where $V$ is a $C^{2,\alpha}(\mathbb{R}^{n},\mathbb{R}^{n})$ vector field.

The Jacobian matrix of the transformation

$y:=y(t,\Omega)=x+tV(x),\,\,\,\,\,\,x\in\Omega,\,\,t\,\,\text{small}$ (2.2)

is

$D_{y}=I+tD_{V},$

where $I$ is the identity matrix and $(D_{V})_{ij}=\frac{\partial
V_{i}}{\partial x_{j}}$. By Jacobi’s formula, for small t, the Jacobian
determinant is given by

$J(t)=1+t\,\mathrm{div}\,V.$ (2.3)

It is clear that for t small enough, $J(t)\approx 1$, so $y(t,\Omega)$ is a
diffeomorphism and in this case we can consider its inverse transformation
$x(y)$.
In particular we can write the measure of $\Omega_{t}$ in terms of the
perturbations defined before

$|\Omega_{t}|=\int_{\Omega}J(t)\,dx=|\Omega|+t\int_{\Omega}\mathrm{div}\,V\,dx.$

######  Definition 2.1.

$y(t,\Omega)$ is said to be volume preserving of the first order if

$\int_{\Omega}\mathrm{div}\,V\,dx=0.$

Let $w\in H^{1}(\Omega_{t})$ and let us consider the following energy
functional

$\mathcal{E}(\Omega_{t},w)=\int_{\Omega_{t}}|\nabla_{y}w(y,t)|^{2}\,dy-2\int_{\Omega_{t}}w(y,t)\,dy+\beta\int_{\partial\Omega_{t}}w^{2}(y,t)\,d\sigma_{t},$
(2.4)

where with $\nabla_{y}$ we denoted the gradient operator with respect to $y$
and $d\sigma_{t}$ is the surface element of $\Omega_{t}$.
A critical point $u\in H^{1}(\Omega_{t})$ of (2.4) satisfies the Euler-
Lagrange equations

$\begin{cases}\Delta_{y}u(y,t)+1=0&\mbox{in}\ \Omega_{t}\vspace{0.2cm}\\\
\dfrac{\partial u}{\partial\nu_{t}}(y,t)+\beta u(y,t)=0&\mbox{on}\
\partial\Omega_{t}\vspace{0.2cm},\end{cases}$ (2.5)

where $\Delta_{y}$ is the Laplacian operator with respect to $y$ and $\nu_{t}$
is the outer normal to $\partial\Omega_{t}$.
We want to transform the integrals in (2.4) in integrals onto $\Omega$ and
$\partial\Omega$. Indeed by a change of variables, using the inverse function
$x(y)$ (which exists for small t), we get

$\begin{split}\mathcal{E}(\Omega,u)&=\int_{\Omega}\frac{\partial u}{\partial
x_{i}}(x+tV(x),t)\frac{\partial u}{\partial x_{j}}(x+tV(x),t)\frac{\partial
x_{i}}{\partial y_{k}}\frac{\partial x_{j}}{\partial y_{k}}J(t)\,dx\\\
&-2\int_{\Omega}u(x+tV(x),t)J(t)\,dx+\beta\int_{\partial\Omega}u^{2}(x+tV(x),t)m(t)\,d\sigma.\end{split}.$
(2.6)

Here $m(t)$ is the index of deformation when passing from $d\sigma$ to
$d\sigma_{t}$ (with $d\sigma$ being the surface element of $\Omega$). If we
define the tangential divergence of the vector field V as follows

$\mathrm{div}_{\partial\Omega}\,V:=\mathrm{div}\,V-\nu\cdot D_{V}\nu,$ (2.7)

then, up to first order terms, $m(t)$ can be approximated by (See [BW],
section $2.2.2$)

$m(t)=1+t\mathrm{div}_{\partial\Omega}\,V.$ (2.8)

If we denote by

$\tilde{u}(t):=u(x+tV(x),t)$ (2.9)

and

$A=(A_{ij}(t)):=\frac{\partial x_{i}}{\partial y_{k}}\frac{\partial
x_{j}}{\partial y_{k}}J(t),$ (2.10)

we can write (2.6) in a more concise form

$\mathcal{E}(t):=\int_{\Omega}\nabla\tilde{u}(t)A\nabla\tilde{u}(t)\,dx-2\int_{\Omega}\tilde{u}(t)J(t)\,dx+\beta\int_{\partial\Omega}\tilde{u}^{2}(t)m(t)\,d\sigma.$
(2.11)

If we simplify one more time the notations and indicate by

$L_{A}=\frac{\partial}{\partial x_{j}}(A_{ij}(t)\frac{\partial}{\partial
x_{i}})$ (2.12)

and

$\partial_{\nu_{A}}=\nu_{i}A_{ij}(t)\frac{\partial}{\partial x_{j}},$ (2.13)

then the transformed function $\tilde{u}(t)$ solves the Eulero-Lagrange
equations

$\begin{cases}L_{A}\tilde{u}(t)+J(t)=0&\mbox{in}\ \Omega\vspace{0.2cm}\\\
\partial_{\nu_{A}}\tilde{u}(t)+\beta m(t)\tilde{u}(t)=0&\mbox{on}\
\partial\Omega\vspace{0.2cm}.\end{cases}$ (2.14)

A key role will be played by (2.9). If we expand it in a neighbourhood of
$t=0$ we have

$\tilde{u}(t)=\tilde{u}(0)+t\dot{\tilde{u}}(0)+o(t).$

We remark that the dot notation stands for the derivative with respect to t.
In particular the first two coefficients will be

$\tilde{u}(0):=u(x)$ (2.15)

and

$\begin{split}\dot{\tilde{u}}(0)=\left[\frac{d}{dt}\tilde{u}(t)\right]_{t=0}&=\left[\frac{\partial\tilde{u}}{\partial
t}(t)+V(x)\cdot\nabla\tilde{u}(t)\right]_{t=0}\\\
&=\frac{\partial\tilde{u}}{\partial t}(0)+V\cdot\nabla u.\end{split}$ (2.16)

######  Definition 2.2.

We will call shape derivative of $\tilde{u}$, and it will be denoted by
$u^{\prime}$, the following function

$u^{\prime}(x):=\frac{\partial\tilde{u}}{\partial t}(0).$ (2.17)

Hence we can write (2.16) in this way

$\dot{\tilde{u}}(0)=u^{\prime}+V\cdot\nabla u.$ (2.18)

Besides the expansion of $\tilde{u}(t)$, it will be helpful to write the
taylor series of some other of the function seen until now. Next Lemma will
collect all the necessary coefficients of the expansions just mentioned

######  Lemma 2.3.

We have that

$\displaystyle 1)$ $\displaystyle J(0)=1,$ $\displaystyle 2)$
$\displaystyle\dot{J}(0)=\mathrm{div}\,V,$ (2.19) $\displaystyle 3)$
$\displaystyle m(0)=1,$ $\displaystyle 4)$
$\displaystyle\dot{m}(0)=\mathrm{div}_{\partial\Omega}\,V,$ $\displaystyle 5)$
$\displaystyle\tilde{u}(0)=u(x),$ $\displaystyle 6)$
$\displaystyle\dot{\tilde{u}}(0)=u^{\prime}+V\cdot\nabla u,$ $\displaystyle
7)$ $\displaystyle A_{ij}(0)=\delta_{ij},$ $\displaystyle 8)$
$\displaystyle\dot{A}_{ij}(0)=\mathrm{div}\,V\delta_{ij}-\frac{\partial
V_{i}}{\partial x_{j}}-\frac{\partial V_{j}}{\partial x_{i}}.$

###### Proof.

To compute $1),2),3),4)$ it is sufficient to differentiate (2.3) and (2.8) and
evaluate for $t=0$. $5)$ and $6)$ are given by (2.15) and (2.18) respectively.
Some more effort will be needed for the matrix $A$, defined in (2.10).
Remembering that the Jacobian matrix of the transformation $y(t,\Omega)$ is

$D_{y}=I+tD_{V},$

if $t$ is small enough, we have

$\begin{split}\frac{\partial x_{i}}{\partial
y_{k}}=(D_{y}^{-1})_{ik}&=(I+tD_{V})^{-1}_{ik}\\\
&=(I-tD_{V}+o(t))_{ik}=\delta_{ik}-t\frac{\partial V_{k}}{\partial
x_{i}}+o(t).\end{split}$

This allows us to obtain

$A_{ij}(0)=\delta_{ik}\delta_{jk}=\delta_{ij}$

and

$\begin{split}\dot{A}_{ij}(0)&=\bigg{[}\frac{\partial x_{i}}{\partial
y_{k}}\frac{\partial x_{j}}{\partial
y_{k}}\dot{J}(t)+\bigg{(}\frac{d}{dt}\frac{\partial x_{i}}{\partial
y_{k}}\bigg{)}\frac{\partial x_{j}}{\partial y_{k}}J(t)+\frac{\partial
x_{i}}{\partial y_{k}}\bigg{(}\frac{d}{dt}\frac{\partial x_{j}}{\partial
y_{k}}\bigg{)}J(t)\bigg{]}_{t=0}\\\ &=\mathrm{div}V\delta_{ij}-\frac{\partial
V_{i}}{\partial x_{j}}-\frac{\partial V_{j}}{\partial x_{i}}.\end{split}$

∎

We want to find the equations that are solved by $u^{\prime}$ in $B_{R}$ and
on its boundary.

###  2.2 An equation for $u^{\prime}$ in $B_{R}$

Let us consider problem (1.10). It is well known that it admits a unique and
positive solution, given by

$u(x)=\frac{R}{\beta n}+\frac{1}{2n}(R^{2}-|x|^{2}),$ (2.20)

which is a radial and strictly concave function, whose maximum and minimum are
achieved in $0$ and on $\partial B_{R}$ respectevely. More precisely

$u_{\max}=u(0)=\frac{R}{\beta
n}+\frac{R^{2}}{2n},\,\,\,\,\,\,\,u_{\min}=u(R)=\frac{R}{\beta n}.$ (2.21)

In order to prove next proposition, it will be useful to keep in mind the
gradient and the Hessian matrix of (2.20). The gradient is

$\nabla u(x)=-\frac{x}{n}.$ (2.22)

In particular, if $x\in\partial B_{R}$, being $\nu=\frac{x}{R}$ the outer unit
normal to the boundary of $B_{R}$, then

$\nabla u(x)=-\frac{R}{n}\nu,\,\,\,\,\,\,\,\,\,\,\frac{\partial
u}{\partial\nu}=-\frac{R}{n}.$ (2.23)

The Hessian matrix is clearly negative definite and it is given by

$\text{Hess}_{u}(x)=-\frac{I}{n},$ (2.24)

where I is the identity matrix.

######  Proposition 2.4.

Let $V$ be a $C^{2,\alpha}(\mathbb{R}^{n},\mathbb{R}^{n})$ vector field and
$R>0$ a positive real number. The function $u^{\prime}$ solves the following
boundary value problem in the ball with radius $R$

$\begin{cases}\Delta u^{\prime}=0&\mbox{in}\ B_{R}\vspace{0.2cm}\\\
\displaystyle\frac{\partial u^{\prime}}{\partial\nu}+\beta
u^{\prime}=\displaystyle\left(\frac{1+\beta R}{n}\right)(V\cdot\nu)&\mbox{on}\
\partial B_{R}\vspace{0.2cm},\end{cases}$ (2.25)

where $\nu$ is the outer unit normal to $\partial B_{R}$.

###### Proof.

If we differentiate the first equation in (2.14) with respect to $t$ and
evaluate for $t=0$, we obtain

$L_{A(0)}\dot{\tilde{u}}(0)+L_{\dot{A}(0)}\tilde{u}(0)+\dot{J}(0)=0.$

It will be helpful to write explicitly $L_{A}$ and $L_{\dot{A}}$. By applying
Lemma 2.3 we have

$\begin{split}\frac{\partial}{\partial
x_{j}}\bigg{(}\delta_{ij}\frac{\partial}{\partial
x_{i}}\bigg{)}&\bigg{(}u^{\prime}+V_{k}\frac{\partial u}{\partial
x_{k}}\bigg{)}\\\ &+\frac{\partial}{\partial x_{j}}\bigg{(}\frac{\partial
V_{k}}{\partial x_{k}}\delta_{ij}\frac{\partial}{\partial
x_{i}}-\frac{\partial V_{j}}{\partial x_{i}}\frac{\partial}{\partial
x_{i}}-\frac{\partial V_{i}}{\partial x_{j}}\frac{\partial}{\partial
x_{i}}\bigg{)}u+\frac{\partial V_{k}}{\partial x_{k}}=0.\end{split}$

Hence

$\begin{split}\frac{\partial}{\partial x_{j}}\bigg{(}\frac{\partial
u^{\prime}}{\partial x_{j}}+\frac{\partial}{\partial
x_{j}}&\bigg{(}V_{k}\frac{\partial u}{\partial x_{k}}\bigg{)}\bigg{)}\\\
&+\frac{\partial}{\partial x_{j}}\bigg{(}\frac{\partial V_{k}}{\partial
x_{k}}\frac{\partial u}{\partial x_{j}}-\frac{\partial V_{j}}{\partial
x_{i}}\frac{\partial u}{\partial x_{i}}-\frac{\partial V_{i}}{\partial
x_{j}}\frac{\partial u}{\partial x_{i}}\bigg{)}+\frac{\partial V_{k}}{\partial
x_{k}}=0.\end{split}$

Renaming the indexes $k$

$\begin{split}\frac{\partial^{2}u^{\prime}}{\partial
x_{j}^{2}}&+\frac{\partial^{2}V_{i}}{\partial x_{j}^{2}}\frac{\partial
u}{\partial x_{i}}+2\frac{\partial V_{i}}{\partial
x_{j}}\frac{\partial^{2}u}{\partial x_{j}\partial
x_{i}}+V_{i}\frac{\partial^{3}u}{\partial^{2}x_{j}\partial x_{i}}\\\
&+\frac{\partial^{2}V_{i}}{\partial x_{j}\partial x_{i}}\frac{\partial
u}{\partial x_{j}}+\frac{\partial V_{i}}{\partial
x_{i}}\frac{\partial^{2}u}{\partial
x_{j}^{2}}-\frac{\partial^{2}V_{j}}{\partial x_{j}\partial
x_{i}}\frac{\partial u}{\partial x_{i}}\\\ &-\frac{\partial V_{j}}{\partial
x_{i}}\frac{\partial^{2}u}{\partial x_{j}\partial
x_{i}}-\frac{\partial^{2}V_{i}}{\partial x_{j}^{2}}\frac{\partial u}{\partial
x_{i}}-\frac{\partial V_{i}}{\partial x_{j}}\frac{\partial^{2}u}{\partial
x_{j}\partial x_{i}}+\frac{\partial V_{i}}{\partial x_{i}}=0.\end{split}$

Considering that $(\mathrm{Hess}_{u}(x))_{ij}=\frac{\partial^{2}u}{\partial
x_{j}\partial x_{i}}=0$ whenever $i\neq j$ and the fact that
$\frac{\partial^{2}u}{\partial x_{j}^{2}}=\Delta u=-1$, we have

$\frac{\partial u^{\prime}}{\partial
x_{j}^{2}}+\frac{\partial^{2}V_{i}}{\partial x_{j}\partial
x_{i}}\frac{\partial u}{\partial x_{j}}-\frac{\partial^{2}V_{j}}{\partial
x_{j}\partial x_{i}}\frac{\partial u}{\partial x_{i}}=0.$

Eventually, by Schwarz Theorem, we obtain $\Delta u^{\prime}=0$ in $B_{R}$.
If we now differentiate the boundary conditions in (2.14) and evaluate them
for $t=0$, then

$\partial_{\nu_{A(0)}}\dot{\tilde{u}}(0)+\partial_{\nu_{\dot{A}(0)}}\tilde{u}(0)+\beta\dot{m}(0)\tilde{u}(0)+\beta
m(0)\dot{\tilde{u}}(0)=0.$ (2.26)

Let us compute every term in the previous equation. Considering the boundary
conditions satisfied by $u$, Lemma (2.3) and (2.23) we get

$\partial_{\nu_{A(0)}}\dot{\tilde{u}}(0)=\frac{\partial
u^{\prime}}{\partial\nu}+\frac{\partial}{\partial\nu}(V\cdot\nabla
u)=\frac{\partial
u^{\prime}}{\partial\nu}-\frac{R}{n}\nabla(V\cdot\nu)\cdot\nu,$
$\begin{split}\partial_{\nu_{\dot{A}(0)}}\tilde{u}(0)&=\frac{\partial
u}{\partial\nu}\,\mathrm{div}\,V-\nu\cdot D_{V}\nabla u-\nabla u\cdot
D_{V}\nu\\\ &=-\frac{R}{n}\mathrm{div}\,V+\frac{2R}{n}\nu\cdot
D_{V}\nu,\end{split}$ $\beta\dot{m}(0)\tilde{u}(0)=\beta
u\,\mathrm{div}_{\partial B_{R}}V=\frac{R}{n}\mathrm{div}V-\frac{R}{n}\nu\cdot
D_{V}\nu,$ $\beta m(0)\dot{\tilde{u}}(0)=\beta u^{\prime}+\beta V\cdot\nabla
u=\beta u^{\prime}-\frac{R\beta}{n}V\cdot\nu.$

Substituting in (2.26) we have

$\frac{\partial u^{\prime}}{\partial\nu}+\beta
u^{\prime}=\frac{R}{n}\nabla(V\cdot\nu)\cdot\nu-\frac{R}{n}\nu\cdot
D_{V}\nu+\frac{\beta R}{n}V\cdot\nu.$

Now

$\frac{R}{n}\nabla(V\cdot\nu)=\frac{R}{n}\nu\cdot
D_{V}\nu+\frac{1}{n}V\cdot\nu.$

Hence

$\frac{\partial u^{\prime}}{\partial\nu}+\beta
u^{\prime}=\frac{1}{n}V\cdot\nu+\frac{\beta R}{n}V\cdot\nu=\left(\frac{1+\beta
R}{n}\right)(V\cdot\nu).$

∎

As a consequence of the previous proposition, we deduce that

######  Corollary 2.5.

If $\,V\in C^{2,\alpha}(\mathbb{R}^{n},\mathbb{R}^{n})$ is a vector field
volume preserving of the first order, then the solution $u^{\prime}$ to (2.25)
is a function with zero mean value in $B_{R}$, that is

$\int_{B_{R}}u^{\prime}\,dx=\int_{\partial B_{R}}u^{\prime}\,d\sigma=0.$

###### Proof.

Let us integrate the first equation in (1.10)

$0=\int_{B_{R}}\Delta u^{\prime}\,dx=\int_{B_{R}}\mathrm{div}(\nabla
u^{\prime})\,dx=\int_{\partial B_{R}}\frac{\partial
u^{\prime}}{\partial\nu}\,d\sigma.$ (2.27)

Since $V$ is volume preserving of the first order

$\int_{B_{R}}\mathrm{div}\,V\,dx=\int_{\partial B_{R}}V\cdot\nu\,d\sigma=0.$
(2.28)

So integrating the second equation on the boundary, by (2.27) and (2.28)

$\int_{\partial B_{R}}u^{\prime}\,d\sigma=\frac{1}{\beta}\left(\frac{1+\beta
R}{n}\right)\int_{\partial B_{R}}V\cdot\nu\,d\sigma=0.$

In conclusion, being $u^{\prime}$ a harmonic function, by the mean value
theorem

$\int_{B_{R}}u^{\prime}\,dx=\frac{R}{n}\int_{\partial
B_{R}}u^{\prime}\,d\sigma=0.$

∎

##  3 Main results

The shape derivative $u^{\prime}$, solution to problem (2.25) will play a
central role to prove the desired results.
Let us recall the definition of the $L^{\infty}$ and $L^{p}$ functionals. If
$\Omega$ is a bounded and simply connected open subset of $\mathbb{R}^{n}$ and
$u$ is the solution to problem (1.1), we will denote the
$L^{\infty}$-functional over $\Omega$ by

$M(\Omega)=\|u\|_{L^{\infty}(\Omega)},$

and, if $p\geq 1$, the $L^{p}$-functional

$F_{p}(\Omega)=\|u\|_{L^{p}(\Omega)}^{p}.$

###  3.1 Shape derivative of the $L^{\infty}$-norm

We prove the next result following the proof that can be found in [HLP].

######  Theorem 3.1.

Let $B_{R}$ be a ball centered at the origin with radius $R>0$. Then for every
$C^{2,\alpha}(\mathbb{R}^{n},\mathbb{R}^{n})$ vector field $V$, the shape
derivative of $M$ at $B_{R}$ in any direction $V$ exists and it is given by

$M^{\prime}(B_{R},V)=u^{\prime}(0),$

where $0$ is the maximum point of (2.20) and $u^{\prime}$ is the solution to
(2.25).

###### Proof.

Let us perturb the ball $B_{R}$ as we have seen in the previous section

$B_{R,t}=\\{y=x+tV(x):x\in B_{R},\,\,t\,\,\text{small enough}\\},$

and consider problem (2.5), with $B_{R,t}$ in place of $\Omega_{t}$.
Being $x=0$ the unique maximum point of (2.20), then $\nabla u(0,0)=0$. On the
other hand, by the strict concavity of the torsion function $u$ on $B_{R}$,
the matrix

$D_{y}\nabla u(0,0)=\mathrm{Hess}_{u}(0)$

is invertible, since $\mathrm{Hess}_{u}$ is negative definite. Hence by the
implicit function theorem, in a neighbourhood of the origin and for $t$ small
enough, there exists a unique $y_{t}$ such that $\nabla u(y_{t},t)=0$.
Moreover the function $t\longrightarrow y_{t}$ is differentiable and $y_{t}$
must be a maximum, so $M(B_{R,t})=u(y_{t},t)$.
We want to prove that

$\lim_{t\to 0}\frac{M(B_{R,t})-M(B_{R})}{t}=u^{\prime}(0),$

where

$\begin{split}\frac{M(B_{R,t})-M(B_{R})}{t}&=\frac{u(y_{t},t)-u(0,0)}{t}\\\
&=\frac{u(y_{t},t)-u(0,t)}{t}+\frac{u(0,t)-u(0,0)}{t}.\end{split}$

By the differentiability of the map $t\longrightarrow u(\cdot,t)$ and the the
fact that $\nabla u(0,0)=0$, we have

$\lim_{t\to
0}\frac{u(0,t)-u(0,0)}{t}=\frac{d}{dt}\left[u(0,t)\right]_{t=0}=u^{\prime}(0)+V(0)\cdot\nabla
u(0,0)=u^{\prime}(0).$

Furthermore, by the differentiability of $t\longrightarrow y_{t}$, by Lagrange
theorem on the segment $[0,y_{t}]$, the mean value property of $\nabla
u(\cdot,t)$ and the regularity of $u(\cdot,t)$, we get

$\begin{split}\frac{u(y_{t},t)-u(0,t)}{t}&=\nabla
u(\xi_{t},t)\frac{y_{t}}{t}\\\
&=\left(\frac{1}{|B_{r}(\xi_{t})|}\int_{B_{r}(\xi_{t})}\nabla
u(y,t)\,dy\right)\frac{y_{t}}{t},\end{split}$

with $\xi_{t}$ a suitable point in $[0,y_{t}]$. Hence

$\lim_{t\to 0}\frac{u(y_{t},t)-u(0,t)}{t}=\nabla
u(0,0)\left[\frac{dy_{t}}{dt}\right]_{t=0}=0.$

This conclude the proof. ∎

######  Corollary 3.2.

The ball is a critical shape for the functional $M$, for every $V\in
C^{2,\alpha}(\mathbb{R}^{n},\mathbb{R}^{n})$ which is volume preserving of the
first order, i.e.

$M^{\prime}(B_{R},V)=u^{\prime}(0)=0.$

###### Proof.

As a consequence of Lemma 2.5 and Theorem 3.1, applying the mean value
theorem, we have that

$M^{\prime}(B_{R},V)=u^{\prime}(0)=\frac{1}{n\omega_{n}R^{n-1}}\int_{\partial
B_{R}}u^{\prime}(x)\,d\sigma=\frac{1}{\omega_{n}R^{n}}\int_{B_{R}}u^{\prime}\,dx=0.$

∎

###  3.2 Shape derivative of the $L^{p}$-norm

Next theorem will be a straightforward computation of the shape derivative of
the functional $F_{p}(\Omega)$.

######  Theorem 3.3.

For every $C^{2,\alpha}(\mathbb{R}^{n},\mathbb{R}^{n})$ vector field $V$, the
shape derivative of $F_{p}$ at $\Omega$ in any direction $V$ exists and it is
given by

$F^{\prime}_{p}(\Omega,V)=p\int_{\Omega}u^{p-1}u^{\prime}\,dx+\int_{\partial\Omega}u^{p}(V\cdot\nu)\,d\sigma,$

where $u^{\prime}$ is the shape derivative of $u$, solution to problem (1.1),
and $\nu$ is the outer unit normal to $\partial\Omega$.

###### Proof.

Let $u$ be the solution to the perturbed problem

$\begin{cases}-\Delta_{y}u(y,t)=1&\mbox{in}\ \Omega_{t}\vspace{0.2cm}\\\
\dfrac{\partial u}{\partial\nu_{t}}(y,t)+\beta u(y,t)=0&\mbox{on}\
\partial\Omega_{t}\vspace{0.2cm},\end{cases}$

where $\Omega_{t}$ is the perturbed domain defined in (2.1). Then

$F_{p}(\Omega_{t})=\int_{\Omega_{t}}u^{p}(y,t)\,dy=\int_{\Omega}\tilde{u}^{p}(t)J(t)\,dx,$

with $\tilde{u}(t)=\tilde{u}(x+tV(x),t)$ and $J(t)$ the Jacobian determinant
as in (2.3).

Then it is possible to differentiate under the sign of integral and

$\begin{split}\frac{d}{dt}F_{p}(\Omega_{t})&=\frac{d}{dt}\int_{\Omega}\tilde{u}^{p}(t)J(t)\,dx=\int_{\Omega}\frac{d}{dt}[\tilde{u}^{p}(t)J(t)]\,dx\\\
&=p\int_{\Omega}\tilde{u}^{p-1}(t)\dot{\tilde{u}}(t)J(t)\,dx+\int_{\Omega}\tilde{u}^{p}(t)\dot{J}(t)\,dx.\end{split}$

Evaluating this derivative for $t=0$

$\left[\frac{d}{dt}F_{p}(\Omega_{t})\right]_{t=0}=p\int_{\Omega}\tilde{u}^{p-1}(0)\dot{\tilde{u}}(0)J(0)\,dx\\\
+\int_{\Omega}\tilde{u}^{p}(0)\dot{J}(0)\,dx.$

Applying Lemma (2.3), we get

$\begin{split}F^{\prime}_{p}(\Omega,V)&=p\int_{\Omega}u^{p-1}u^{\prime}\,dx+p\int_{\Omega}u^{p-1}V\cdot\nabla
u\,dx+\int_{\Omega}u^{p}\mathrm{div}\,V\,dx\\\
&=p\int_{\Omega}u^{p-1}u^{\prime}\,dx+p\int_{\Omega}u^{p-1}V\cdot\nabla
u\,dx+\int_{\partial\Omega}u^{p}(V\cdot\nu)\,d\sigma\\\
&-p\int_{\Omega}u^{p-1}V\cdot\nabla
u\,dx=p\int_{\Omega}u^{p-1}u^{\prime}\,dx+\int_{\partial\Omega}u^{p}(V\cdot\nu)\,d\sigma.\end{split}$

∎

When $\Omega=B_{R}$, we can use the symmetry properties of (2.20) and the
property of $u^{\prime}$ to be a zero mean function, to prove that

######  Corollary 3.4.

The ball $B_{R}$ centered at the origin with radius $R>0$ is a critical shape
for the functional $F_{p}$, for every $1\leq p<+\infty$ and every vector field
$V\in C^{2,\alpha}(\mathbb{R}^{n},\mathbb{R}^{n})$ which is volume preserving
of the first order, i.e.

$F^{\prime}_{p}(B_{R},V)=0.$

###### Proof.

By previous theorem we know that

$F^{\prime}_{p}(B_{R},V)=p\int_{B_{R}}u^{p-1}u^{\prime}\,dx+\int_{\partial
B_{R}}u^{p}(V\cdot\nu)\,d\sigma.$

Being $u$ constant on the boundary and $V$ a vector field volume preserving of
the first order

$\int_{\partial B_{R}}u^{p}(V\cdot\nu)\,d\sigma=u_{\min}^{p}\int_{\partial
B_{R}}(V\cdot\nu)\,d\sigma=u_{\min}^{p}\int_{B_{R}}\mathrm{div}\,V=0,$

where $u_{\min}=u(R)=\frac{R}{\beta n}$. By corollary 3.2, we know that
$u^{\prime}(0)=0$ and so by the mean value theorem, we have that

$\int_{\partial
B_{r}}u^{\prime}(x)\,d\mathcal{H}^{n-1}=u^{\prime}(0)n\omega_{n}r^{n-1}=0,$

for every $r\in[0;R]$. Eventually, applying the Coarea Formula

$\begin{split}&p\int_{B_{R}}u^{p-1}u^{\prime}\,dx=p\int_{B_{R}}\left(\frac{R}{\beta
n}-\frac{1}{2n}(R^{2}-|x|^{2})\right)^{p-1}u^{\prime}\,dx\\\
=&\,p\int_{0}^{R}\left(\frac{R}{\beta
n}-\frac{1}{2n}(R^{2}-r^{2})\right)^{p-1}\int_{\partial
B_{r}}u^{\prime}(x)\,d\sigma\,dr=0.\end{split}$

Hence

$F^{\prime}_{p}(B_{R},V)=0.$

∎

## References

  * [ANT] A. Alvino, C. Nitsch, C. Trombetti, A Talenti comparison result for solutions to elliptic problems with robin boundary conditions, ArXiv (2019): 1-15.
  * [BG] D. Bucur, A. Giacomini,  Faber-Krahn inequalities for the Robin-Laplacian: A free discontinuity approach, Arch. Rational Mech. Anal. 218, no. 2, (2015), 757-824.
  * [BW] C. Bandle, A. Wagner, Second domain variation for problems with Robin boundary conditions, ArXiv (2015): 1-42.
  * [HLP] A. Henrot, L. Lucardesi, G. Philippin, On two functionals involving the maximum of the torsion function, ArXiv (2017): 1-15.
  * [HP] A. Henrot, M. Pierre, Variation et optimisation de formes, Mathématique & Applications, 48, Springer, (2005), 167-178.
  * [K] S. Kesavan, Symmetrization & applications, Series in Analysis, 3, World Scientific Publishing Co. Pte. Ltd., (2006).
  * [L] G. G. Lorentz, Some new functional spaces, Ann. of Math. (2) 51, (1950), 37-55.
  * [P] L.E. Payne, Isoperimetric inequalities and their applications, SIAM Rev., 9, (1967), 453-488.
  * [PPT] G. Paoli, G. Piscitelli, L. Trani, Sharp estimates for the first p-Laplacian eigenvalue and for the p-torsional rigidity on convex sets with holes, ESAIM: COCV, (in press).
  * [PS] G. Pólya, G. Szegö, Isoperimetric Inequalities in Mathematical Physics, Annals of Mathematics Studies, no. 27, Princeton University Press, (1951).
  * [T] G. Talenti, Elliptic equations and rearrangements, Ann. Scuola Norm. Sup. Pisa Cl. Sci. (4) 3, no. 4, (1976), 697-718.


††thanks<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
#  Neutrinoless double beta decays tell nature of right-handed neutrinos

Takehiko Asaka  Department of Physics, Niigata University, Niigata 950-2181,
Japan  Hiroyuki Ishida  KEK Theory Center, IPNS, Tsukuba, Ibaraki 305-0801,
Japan  Kazuki Tanaka  Graduate School of Science and Technology, Niigata
University Niigata, 950-2181, Japan

###### Abstract

We consider the minimal seesaw model, the Standard Model extended by two
right-handed neutrinos, for explaining the neutrino masses and mixing angles
measured in oscillation experiments. When one of right-handed neutrinos is
lighter than the electroweak scale, it can give a sizable contribution to
neutrinoless double beta ($0\nu\beta\beta$) decay. We show that the detection
of the $0\nu\beta\beta$ decay by future experiments gives a significant
implication to the search for such light right-handed neutrino.

††preprint: KEK-TH-2292

The Standard Model (SM) of the particle physics preserve two accidental global
symmetries in the (classical) Lagrangian, namely the baryon and lepton number
symmetries. It is well known that these global symmetries are non-
perturbatively broken at the quantum level tHooft:1976rip ; tHooft:1976snw ,
especially at high temperature of the universe Dimopoulos:1978kv ;
Manton:1983nd ; Klinkhamer:1984di ; Kuzmin:1985mm . Even at the quantum level,
however, a baryon minus lepton symmetry, often called $U(1)_{\rm B\mathchar
45\relax L}$ #1#1#1We do not specifically consider the symmetry as the gauge
symmetry., has to be preserved in the SM.

The simplest way to break the $U(1)_{\rm B\mathchar 45\relax L}$ symmetry
without loss of the renormalizability is introducing right-handed neutrinos
(RH$\nu$s) into the SM. Since RH$\nu$s are singlet under the SM gauge
symmetries, we can write the mass term, called Majorana mass term, of it
without conflicting the gauge principle. The Majorana mass term breaks the
lepton number symmetry by two units. Therefore, the phenomena of the lepton
number violation can be a definite signal of the existence of RH$\nu$s.

The existence of RH$\nu$s is not only for the violation of the $U(1)_{\rm
B\mathchar 45\relax L}$ symmetry but also important to solve the origin of the
observed tiny neutrino masses. In the renormalizable Lagrangian with RH$\nu$s,
we can obtain two kind of the neutrino masses, one is called Dirac masses and
another is called Majorana masses. When enough hierarchy between these masses
is realized, we can simply explain the tiny neutrino masses by the seesaw
mechanism Minkowski:1977sc ; Yanagida:1979as ; Yanagida:1980xy ; Ramond:1979 ;
GellMann:1980vs ; Glashow:1979 ; Mohapatra:1979ia . In addition, the violation
of $U(1)_{\rm B\mathchar 45\relax L}$ can seeds the origin of the baryon
asymmetry of the universe #2#2#2There are a bunch of possibilities to provide
the baryon asymmetry through the lepton number violation. But the detail of
the mechanism is independent of the discussions below. .

One of the most promising signals of the $U(1)_{\rm B\mathchar 45\relax L}$
violation is the neutrinoless double beta decay, which breaks the lepton
number by two units while keeping the baryon number. (See, for example,
articles Doi:1985dx ; Pas:2015eia ; DellOro:2016tmg ; Dolinski:2019nrj .) The
rate of the decay is characterized by the effective mass defined by the
neutrino masses and mixing angles. When we simply add Majorana masses of three
(active or left-handed) neutrinos which are responsible for the neutrino
oscillation into the SM, the effective mass can be predicted depending on the
lightest active neutrino mass together with the unknown CP violating phases.

In view of the fundamental models for the origin of the neutrino masses, the
mass of the lightest active neutrino cannot be determined uniquely, leading to
different predictions on the effective mass. It should be noted that the
effective mass can vanish in the normal hierarchy (NH) case of the active
neutrinos in a certain parameter region. In such a case, the contribution from
new physics (other than active neutrinos) including RH$\nu$s would be more
important for the detection. So far, no neutrinoless double beta decay is
detected and the upper bounds on the effective mass have been imposed by
various experiments.#3#3#3In a recent analysis 1833580 , the differential rate
of two neutrinoless double beta decay is discussed to constrain mixing
elements of RH$\nu$s with masses at $\mathcal{O}(0.1\mathchar 45\relax
10)~{}{\rm MeV}$. The most stringent bound at present is $61$-$165$ meV by the
KamLAND-Zen experiment KamLAND-Zen:2016pfg . Since this limit is approaching
to the predicted range in the inverted hierarchy (IH) case, the experimental
results in near future can give us some implication on RH$\nu$s.

There are several interesting possibilities that the effective mass can be
significantly modified due to the destructive or constructive contribution
from RH$\nu$s. This additional contribution becomes important when the masses
of RH$\nu$s are smaller or comparable to the typical scale of Fermi momentum
in the decaying nucleus ($\sim\mathcal{O}(100)~{}{\rm MeV}$).

Recently, we have pointed out one interesting possibility is that RH$\nu$ may
hide one of the neutrinoless double beta decay processes Asaka:2020wfo ;
Asaka:2020lsx  (see also Refs. Halprin:1983ez ; Leung:1984vy ). This is due to
the destructive contribution of RH$\nu$ to the effective mass. Note that the
impact of RH$\nu$ does depend on the decaying nuclei. If this is the case, the
mixing elements of RH$\nu$ can be predicted in terms of its mass in a certain
range which is a good target of future search experiments.

In this paper, we project out the consequences of the opposite situation,
namely the case when the neutrinoless double beta decay is observed in some
nucleus, and discuss the impacts on the mixing elements of RH$\nu$s.

First of all, let us explain the framework of the present analysis, the
minimal seesaw model. It is the Standard Model extended by two right-handed
neutrinos $\nu_{RI}$ ($I=1,2)$, which Lagrangian is given by

$\displaystyle{\cal L}=$ $\displaystyle{\cal L}_{\rm
SM}+i\overline{\nu_{RI}}\gamma^{\mu}\partial_{\mu}\nu_{RI}$
$\displaystyle-\left(F_{\alpha
I}\overline{L_{\alpha}}\Phi\nu_{RI}+\frac{M_{I}}{2}\overline{\nu_{RI}^{c}}\nu_{RI}+h.c.\right)\,,$
(1)

where $L_{\alpha}=(\nu_{L\alpha},e_{L\alpha})^{T}$ ($\alpha=e,\mu,\tau)$ and
$\Phi$ are the weak doublets of left-handed lepton and Higgs, respectively.
The Yukawa coupling constants and the Majorana masses for neutrinos are
denoted by $F_{\alpha I}$ and $M_{I}$. By assuming that the Dirac masses
$F_{\alpha I}\langle\Phi\rangle$ are much smaller than the Majorana mass
$M_{I}$, the seesaw mechanism works, and the mass eigenstates of neutrinos are
three active neutrinos $\nu_{i}$ ($i=1,2,3)$ with masses $m_{i}$ and two heavy
neutral leptons (HNLs) $N_{I}$ with masses $M_{I}$.

The mass ordering of active neutrinos is not determined by the oscillation
data, and two possibilities, the normal hierarchy (NH) with
$m_{3}>m_{2}>m_{1}=0$ and the inverted hierarchy (IH) with
$m_{2}>m_{1}>m_{3}=0$, are allowed. Note that the lightest active neutrino is
massless in the considering situation. On the other hand, we can take the
masses of HNLs as $M_{2}\geq M_{1}$ without loss of generality. The left-
handed (flavor) neutrinos are then written as

$\displaystyle\nu_{L\alpha}=\sum_{i}U_{\alpha
i}\,\nu_{i}+\sum_{I}\Theta_{\alpha I}\,N_{I}^{c}\,,$ (2)

where $U_{\alpha i}$ is the mixing matrix of active neutrinos called as the
PMNS matrix while $\Theta_{\alpha I}$ is that of HNLs.

One of the most important consequences of the seesaw mechanism is that active
neutrinos and HNLs are both Majorana particles. In this case the lepton number
violating processes are induced by these particles, which is a clear signature
of physics beyond the SM. One promising example is the $0\nu\beta\beta$ decay,
and the quest for the decay is going on by various experiments.

The rate for the $0\nu\beta\beta$ decay mediated by active neutrinos and HNLs
is proportional $|m_{\rm eff}|^{2}$, where $m_{\rm eff}$ is the so-called
effective (neutrino) mass in the $0\nu\beta\beta$ decay. In the minimal seesaw
model it is given by

$\displaystyle m_{\rm eff}=m_{\rm eff}^{\nu}+m_{\rm eff}^{N}\,.$ (3)

Here the first term in the right-hand side represents the contributions from
the active neutrinos, which is given by

$\displaystyle m_{\rm eff}^{\nu}$ $\displaystyle=\sum_{i}U_{ei}^{2}\,m_{i}\,.$
(4)

On the other hand, the contributions from HNLs are expressed as

$\displaystyle m_{\rm eff}^{N}$
$\displaystyle=\sum_{I}\Theta_{eI}^{2}\,M_{I}\,f_{\beta}(M_{I})\,,$ (5)

where $f_{\beta}$ is the suppression factor compared to $m_{\rm eff}^{\nu}$
due to the heaviness of HNLs $M_{I}\gg m_{i}$. Here we apply the result in
Ref. Faessler:2014kka ; Barea:2015zfa  and assume the following form

$\displaystyle
f_{\beta}(M)=\frac{\Lambda_{\beta}^{2}}{\Lambda_{\beta}^{2}+M^{2}}\,,$ (6)

where $\Lambda_{\beta}={\cal O}(10^{2})$ MeV denotes the typical scale of the
Fermi momentum in the $0\nu\beta\beta$ decay. Hereafter we take
$\Lambda_{\beta}=200$ MeV as a representative value.

In this letter we consider the impacts of the detection of the
$0\nu\beta\beta$ decay by future experiments on the properties of HNLs. The
measurement of the decay rate gives the value of $|m_{\rm eff}|$. Note that
$m_{\rm eff}$ is a complex number. First, we consider the case when right-
handed neutrinos possess the hierarchical masses $M_{2}\gg M_{1}$. We then
find that the mixing element $|\Theta_{e1}|^{2}$ of the lighter HNL is given
by

$\displaystyle\Theta_{e1}^{2}=\frac{m_{\rm eff}-m_{\rm
eff}^{\nu}\left[1-f_{\beta}(M_{2})\right]}{M_{1}\left[f_{\beta}(M_{1})-f_{\beta}(M_{2})\right]}\,.$
(7)

Here we have used the intrinsic relation between mixing elements in the seesaw
mechanism

$\displaystyle 0$
$\displaystyle=\sum_{i}U_{ei}^{2}\,m_{i}+\sum_{I}\Theta_{eI}^{2}\,M_{I}\,.$
(8)

Importantly, the mixing element $|\Theta_{e1}|^{2}$ is given by $m_{\rm eff}$
and $m_{\rm eff}^{\nu}$ together with masses $M_{1}$ and $M_{2}$. This means
that, if $|m_{\rm eff}|$ is found by the detection of the $0\nu\beta\beta$
decay, the range of $|\Theta_{e1}|^{2}$ can be predicted. In practice both
upper and lower bounds on $|\Theta_{e1}|^{2}$ are obtained by varying the
unknown parameters in $m_{\rm eff}^{\nu}$ (i.e., the Majorana phase $\eta$ and
the mass ordering) and the phase of $m_{\rm eff}$.

Figure 1:  Upper and lower bounds on $|\Theta_{e1}|^{2}$ for the NH (red solid
lines) and IH (blue dashed lines) cases. Here $M_{1}=1$ GeV and $M_{2}=200$
GeV.

When $M_{1}=1$ GeV and $M_{2}=200$ GeV, these bounds are shown in Fig. 1 in
terms of the (would-be) observed value of $|m_{\rm eff}|$ denoted by $m_{\rm
eff}^{\rm obs}$. In the present analysis we take the central values of the
mass squared differences, the mixing angles and the Dirac phase in the PMNS
matrix given in Ref. Esteban:2020cvm  for the estimation of $|m_{\rm
eff}^{\nu}|$. We find that $|m_{\rm eff}^{\nu}|=1.45$–$3.68$ meV and 18.6–48.4
meV for the NH and IH cases, respectively. It is found from Eq. (7) that the
lower bound on $|\Theta_{e1}|^{2}$ vanishes when $m_{\rm eff}^{\rm
obs}=|m_{\rm eff}^{\nu}|(1-f_{\beta}(M_{2}))$.

Figure 2:  Upper and lower bounds on $|\Theta_{e1}|^{2}$ for the NH (left) and
IH (right) cases. We take $m_{\rm eff}^{\rm obs}$ =100 meV (red sold lines),
50 meV (blue dashed lines), and 10 meV (green dot-dashed lines). Here
$M_{2}=200$ GeV. The current (conservative) upper bound on $|\Theta_{e1}|^{2}$
from $|m_{\rm eff}|<165$ meV is shown by black solid line (and the light-gray
region is exluded). The dark-gray regions are excluded by the direct search
experiments. The dotted lines shows the sensitivities by the future
experiments. See the detail in the main text.

The predicted range of $|\Theta_{e1}|^{2}$ is shown in Fig. 2 where the
current upper bounds and the sensitivities on $|\Theta_{e1}|^{2}$ by future
search experiments are also shown PIENU:2011aa ; Aguilar-Arevalo:2017vlf ;
NA62:2020mcv ; Blondel:2014bra ; SHiP:2018xqw ; Krasnov:2019kdc ;
Alpigiani:2020tva . We take the (would-be) observed value of the effective
mass as $|m_{\rm eff}|$ = 100 meV, 50 meV, and 10 meV. Importantly, the most
of the predicted range can be tested by the future experiments.

We should note that the understanding of $f_{\beta}(M)$ is important for the
precise prediction of the mixing elements, since it contains the uncertainty
of the order unity. For this purpose the better understanding of the nuclear
matrix elements of the $0\nu\beta\beta$ decay mediated by HNL is crucial.

Next, let us consider the case when the masses of HNLs are degenerate

$\displaystyle M_{1}=M_{2}=M_{N}\,.$ (9)

In this case, the total effective mass is given by

$\displaystyle m_{\rm eff}=m_{\rm
eff}^{\nu}\left[1-f_{\beta}(M_{N})\right]\,,$ (10)

and hence the total value is always smaller than the that from active
neutrinos $|m_{\rm eff}|<|m_{\rm eff}^{\nu}|$ as long as HNLs participate the
$0\nu\beta\beta$ decay. Note that the arguments of $m_{\rm eff}$ and $m_{\rm
eff}^{\nu}$ are the same. In this case, we find the interesting consequences
if $|m_{\rm eff}|$ is measured: First, the mass of degenerate HNLs is
determined depending on the measured value of $|m_{\rm eff}|$ as

$\displaystyle M_{N}=\Lambda_{\beta}\sqrt{\frac{m_{\rm eff}^{\rm obs}}{|m_{\rm
eff}^{\nu}|-m_{\rm eff}^{\rm obs}}}\,.$ (11)

This shows that, once $m_{\rm eff}^{\rm obs}$ is fixed, the unknown Majorana
phase in $m_{\rm eff}^{\nu}$ determines $M_{N}$. Second, the sum of the mixing
elements is found to be

$\displaystyle\left|\Theta_{e1}^{2}+\Theta_{e2}^{2}\right|=\frac{|m_{\rm
eff}^{\nu}|}{\Lambda_{\beta}}\sqrt{\frac{|m_{\rm eff}^{\nu}|-m_{\rm eff}^{\rm
obs}}{m_{\rm eff}^{\rm obs}}}\,.$ (12)

Figure 3:  The degenerate mass $M_{N}$ and mixing element
$|\Theta_{e1}^{2}+\Theta_{e2}^{2}|$ in terms of the observed value $m_{\rm
eff}^{\rm obs}$ in the NH (red solid line) or IH (blue dashed line). We take
the Majorana phase $\eta=0$.

Figure 4:  Range of the mixing element $|\Theta_{e1}^{2}+\Theta_{e2}^{2}|$ in
terms of the degenerate mass $M_{N}$ by taking the Majorana phase
$\eta=0$–$\pi$ in the NH (red solid line) or IH (blue dashed line).

These results are shown in Fig. 3. Here we take the Majorana phase as
$\eta=0$, and $|m_{\rm eff}^{\nu}|$ = 3.54 meV and 48.4 meV for the NH and IH
cases, respectively. It is seen that the observed effective mass $m_{\rm
eff}^{\rm obs}$ of a few $10$ meV corresponds to the Majorana mass
$M_{N}\simeq{\cal O}(0.1-1)$ GeV and the mass ordering is the IH since HNL
contributions are always destructive to the active neutrino ones. The relation
between $M_{N}$ and $|\Theta_{e1}^{2}+\Theta_{e2}^{2}|$ is shown in Fig. 4. We
find that in order to test the degenerate case the improvement of the
sensitivity by future experiments is required especially for the NH case.

Before concluding the paper, we stress the impact of the difference among the
$0\nu\beta\beta$ decay nuclei Asaka:2020lsx . Throughout this paper, we have
assumed the approximated form of the suppression function $f_{\beta}$ to be
Eq. (6) and fixed the typical Fermi momentum as $\Lambda_{\beta}=200~{}{\rm
MeV}$. The important point is that the nuclear matrix elements including the
suppression factor due to HNLs are different depending on the decaying nuclei
used in the $0\nu\beta\beta$ experiments. This effect may be quantified by the
choice the typical Fermi momentum in this analysis.

Figure 5:  Upper and lower bounds of predicted effective mass with
$\tilde{\Lambda}_{\beta}=100~{}{\rm MeV}$ (left) and
$\tilde{\Lambda}_{\beta}=400~{}{\rm MeV}$ (right) in the NH case. We assume
that the effective mass observed in the nucleus with
$\Lambda_{\beta}=200~{}{\rm MeV}$ is $100~{}{\rm meV}$ (red, solid),
$50~{}{\rm meV}$ (blue, bashed), and $10~{}{\rm meV}$ (green, dot-dashed).
Here, we fix $M_{2}=200~{}{\rm GeV}$.

In Fig. 5, we plot the upper and lower values of the predicted effective mass
with different Fermi momentum from $200~{}{\rm MeV}$ while assuming the
$0\nu\beta\beta$ decay is observed at the experiment with
$\Lambda_{\beta}=200~{}{\rm MeV}$ in the NH case. We can obtain similar
behavior straightforwardly in the IH case as well. We take the observed value
of the effective mass to be $100~{}{\rm meV}$, $50~{}{\rm meV}$, or $10~{}{\rm
meV}$. Interestingly, the predicted effective mass can be significantly
enhanced when $\Lambda_{\beta}$ becomes larger enough than $200~{}{\rm MeV}$
and $M_{1}$ gets heavier. By inserting Eq. (7) into the expression of the
effective mass, we can obtain

$\displaystyle\tilde{m}_{\rm eff}$
$\displaystyle=\left[1-\tilde{f}_{\beta}(M_{2})\right]m_{\rm eff}^{\nu}$
$\displaystyle+\left[m_{\rm eff}-m_{\rm
eff}^{\nu}\left[1-f_{\beta}(M_{2})\right]\right]\frac{\tilde{f}_{\beta}(M_{1})-\tilde{f}_{\beta}(M_{2})}{f_{\beta}(M_{1})-f_{\beta}(M_{2})}\,,$
(13)

where $\Lambda_{\beta}=200$ MeV in $f_{\beta}$ but $\Lambda_{\beta}\neq 200$
MeV in $\tilde{f}_{\beta}$ which is denoted as $\tilde{\Lambda}$. Since the
last fraction in the RHS of Eq. (13) can be simplified as

$\displaystyle\frac{\tilde{f}_{\beta}(M_{1})-\tilde{f}_{\beta}(M_{2})}{f_{\beta}(M_{1})-f_{\beta}(M_{2})}=\frac{\tilde{\Lambda}_{\beta}^{2}}{\Lambda_{\beta}^{2}}\frac{\left(\Lambda_{\beta}^{2}+M_{1}^{2}\right)\left(\Lambda_{\beta}^{2}+M_{2}^{2}\right)}{\left(\tilde{\Lambda}_{\beta}^{2}+M_{1}^{2}\right)\left(\tilde{\Lambda}_{\beta}^{2}+M_{2}^{2}\right)}\,,$
(14)

Namely, the effective mass is enhanced as $M_{1}$ gets greater/suppressed than
the typical Fermi momentum in $\tilde{f}_{\beta}$ by the factor
$\tilde{\Lambda}_{\beta}^{2}/\Lambda_{\beta}^{2}$. As clearly seen, since
significant enhancement/suppression could happen depending on the values of
$\Lambda_{\beta}$ due to the contributions from HNLs. Thus, we can claim that
the multiple detection by the $0\nu\beta\beta$ experiments using different
nuclei is crucial to reveal the properties of HNLs.

In conclusions, we have considered the minimal seesaw model with two right-
handed neutrinos. It has been shown that, if the effective mass in the
$0\nu\beta\beta$ decay will be measured by future experiments, the possible
range of the mixing elements for the lighter heavy neutral lepton (right-
handed neutrino) is determined. Especially, when two heavy neutral leptons are
hierarchical and the lighter mass is below the electroweak scale, $N_{1}$ is a
good target of the direct search experiments.

It has also been shown that the predicted effective mass can depend on nucleus
of the experiment. Therefore, comprehensive studies on the neutrinoless double
beta decays in the seesaw mechanism is necessary to extract the concrete
information of the heavy neutral leptons.

## Acknowledgments

The work of T.A. was partially supported by JSPS KAKENHI Grants No. 17K05410,
No. 18H03708, No. 19H05097, and No. 20H01898. The work of H.I. was supported
by JSPS KAKENHI Grant No. 18H03708.

## References

  * (1) G. ’t Hooft, Phys. Rev. Lett. 37 (1976), 8-11 doi:10.1103/PhysRevLett.37.8
  * (2) G. ’t Hooft, Phys. Rev. D 14 (1976), 3432-3450 [erratum: Phys. Rev. D 18 (1978), 2199] doi:10.1103/PhysRevD.14.3432
  * (3) S. Dimopoulos and L. Susskind, Phys. Rev. D 18 (1978), 4500-4509 doi:10.1103/PhysRevD.18.4500
  * (4) N. S. Manton, Phys. Rev. D 28 (1983), 2019 doi:10.1103/PhysRevD.28.2019
  * (5) F. R. Klinkhamer and N. S. Manton, Phys. Rev. D 30 (1984), 2212 doi:10.1103/PhysRevD.30.2212
  * (6) V. A. Kuzmin, V. A. Rubakov and M. E. Shaposhnikov, Phys. Lett. B 155 (1985), 36 doi:10.1016/0370-2693(85)91028-7
  * (7) P. Minkowski, Phys. Lett. B 67, 421 (1977).
  * (8) T. Yanagida, in Proceedings of the Workshop on Unified Theory and Baryon Number of the Universe, edited by.O. Sawada and A. Sugamoto (KEK, Tsukuba, Ibaraki 305- 0801 Japan, 1979) p. 95.
  * (9) T. Yanagida, Prog. Theor. Phys. 64, 1103 (1980).
  * (10) P. Ramond, in Talk given at the Sanibel Symposium, Palm Coast, Fla., Feb. 25-Mar. 2, 1979, preprint CALT-68-709 (retroprinted as hep-ph/9809459).
  * (11) M. Gell-Mann, P. Ramond, and R. Slansky, in Supergravity, edited by.P. van Niewwenhuizen and D. Freedman (North Holland, Amsterdam, 1979) [arXiv:1306.4669 [hep-th]].
  * (12) S. L. Glashow, in Proc. of the Cargése Summer Institute on Quarks and Leptons, Cargése, July 9-29, 1979, eds. M. Lévy et. al, , (Plenum, 1980, New York), p707.
  * (13) R. N. Mohapatra and G. Senjanovic, Phys. Rev. Lett. 44 (1980) 912.
  * (14) M. Doi, T. Kotani and E. Takasugi, Prog. Theor. Phys. Suppl. 83 (1985), 1 doi:10.1143/PTPS.83.1
  * (15) H. Päs and W. Rodejohann, New J. Phys. 17 (2015) no.11, 115010 doi:10.1088/1367-2630/17/11/115010 [arXiv:1507.00170 [hep-ph]].
  * (16) S. Dell’Oro, S. Marcocci, M. Viel and F. Vissani, Adv. High Energy Phys. 2016 (2016), 2162659 doi:10.1155/2016/2162659 [arXiv:1601.07512 [hep-ph]].
  * (17) M. J. Dolinski, A. W. Poon and W. Rodejohann, Ann. Rev. Nucl. Part. Sci. 69 (2019), 219-251 doi:10.1146/annurev-nucl-101918-023407 [arXiv:1902.04097 [nucl-ex]].
  * (18) P. D. Bolton, F. F. Deppisch, L. Gráf and F. Šimkovic, [arXiv:2011.13387 [hep-ph]].
  * (19) A. Gando et al. [KamLAND-Zen], Phys. Rev. Lett. 117 (2016) no.8, 082503 doi:10.1103/PhysRevLett.117.082503 [arXiv:1605.02889 [hep-ex]].
  * (20) T. Asaka, H. Ishida and K. Tanaka, [arXiv:2012.12564 [hep-ph]].
  * (21) T. Asaka, H. Ishida and K. Tanaka, [arXiv:2012.13186 [hep-ph]].
  * (22) A. Halprin, S. T. Petcov and S. P. Rosen, Phys. Lett. B 125 (1983), 335-338 doi:10.1016/0370-2693(83)91296-0
  * (23) C. N. Leung and S. T. Petcov, Phys. Lett. B 145 (1984), 416-420 doi:10.1016/0370-2693(84)90071-6
  * (24) A. Faessler, M. González, S. Kovalenko and F. Šimkovic, Phys. Rev. D 90, no.9, 096010 (2014) doi:10.1103/PhysRevD.90.096010 [arXiv:1408.6077 [hep-ph]].
  * (25) J. Barea, J. Kotila and F. Iachello, Phys. Rev. D 92 (2015), 093001 doi:10.1103/PhysRevD.92.093001 [arXiv:1509.01925 [hep-ph]].
  * (26) I. Esteban, M. C. Gonzalez-Garcia, M. Maltoni, T. Schwetz and A. Zhou, JHEP 09 (2020), 178 doi:10.1007/JHEP09(2020)178 [arXiv:2007.14792 [hep-ph]]; NuFIT 5.0 (2020), www.nu-fit.org.
  * (27) M. Aoki et al. [PIENU], Phys. Rev. D 84 (2011), 052002 doi:10.1103/PhysRevD.84.052002 [arXiv:1106.4055 [hep-ex]].
  * (28) A. Aguilar-Arevalo et al. [PIENU], Phys. Rev. D 97 (2018) no.7, 072012 doi:10.1103/PhysRevD.97.072012 [arXiv:1712.03275 [hep-ex]].
  * (29) E. Cortina Gil et al. [NA62], Phys. Lett. B 807 (2020), 135599 doi:10.1016/j.physletb.2020.135599 [arXiv:2005.09575 [hep-ex]].
  * (30) A. Blondel et al. [FCC-ee study Team], Nucl. Part. Phys. Proc. 273-275 (2016), 1883-1890 doi:10.1016/j.nuclphysbps.2015.09.304 [arXiv:1411.5230 [hep-ex]].
  * (31) C. Ahdida et al. [SHiP], JHEP 04 (2019), 077 doi:10.1007/JHEP04(2019)077 [arXiv:1811.00930 [hep-ph]].
  * (32) I. Krasnov, Phys. Rev. D 100 (2019) no.7, 075023 doi:10.1103/PhysRevD.100.075023 [arXiv:1902.06099 [hep-ph]].
  * (33) C. Alpigiani et al. [MATHUSLA], [arXiv:2009.01693 [physics.ins-det]].


11institutetext: Universidad Politécnica de Madrid (UPM), José Gutiérrez
Abascal, 2, 28006 Madrid, Spain.11email<EMAIL_ADDRESS>22institutetext: School of Science and Technology, IE University, Spain.
22email<EMAIL_ADDRESS>
33institutetext: Department of Mathematics, University of Michigan, Ann Arbor,
MI 48109, USA.
33email<EMAIL_ADDRESS>Centre for Automation and Robotics
(CSIC-UPM), Ctra. M300 Campo Real, Km 0,200, Arganda del Rey - 28500 Madrid,
Spain.44email<EMAIL_ADDRESS>
# Virtual Affine Nonholonomic Constraints††thanks: E. Stratoglou and A.
Anahory have contributed equally. The authors acknowledge financial support
from Grant PID2019-106715GB-C21 funded by MCIN/AEI/ 10.13039/501100011033 and
the LINC Global project from CSIC ”Wildlife Monitoring Bots” INCGL20022. A.B.
was partially supported by NSF grants DMS-1613819 and DMS-2103026, and AFOSR
grant FA 9550-22-1-0215.

Efstratios Stratoglou 1*1* Alexandre Anahory Simoes 2*2* Anthony Bloch 33
Leonardo Colombo 44

###### Abstract

Virtual constraints are relations imposed in a control system that become
invariant via feedback, instead of real physical constraints acting on the
system. Nonholonomic systems are mechanical systems with non-integrable
constraints on the velocities. In this work, we introduce the notion of
virtual affine nonholonomic constraints in a geometric framework. More
precisely, it is a controlled invariant affine distribution associated with an
affine connection mechanical control system. We show the existence and
uniqueness of a control law defining a virtual affine nonholonomic constraint.

###### Keywords:

Nonholonomic Systems Virtual Constraints Nonlinear Control.

##  1 Introduction

Virtual constraints are relations on the configuration variables of a control
system which are imposed through feedback control and the action of actuators,
instead of through physical connections such as gears or contact conditions
with the environment. The class of virtual holonomic constraints became
popular in applications to biped locomotion where it was used to express a
desired walking gait, as well as for motion planning to search for periodic
orbits and its employment in the technique of transverse linearization to
stabilize such orbits.

Virtual nonholonomic constraints are a class of virtual constraints that
depend on velocities rather than only on the configurations of the system.
Those constraints were introduced in [4] to design a velocity-based swing foot
placement in bipedal robots. In particular, this class of virtual constraints
has been used in [5, 6, 7, 8] to encode velocity-dependent stable walking
gaits via momenta conjugate to the unacatuated degrees of freedom of legged
robots and prosthetic legs.

The recent work [9] intorduces an approach to defining rigorously virtual
nonholonomic constraints, but the nonlinear nature of the constraints makes
difficult a thorough mathematical analysis. In the work [1], we provide a
formal definition of linear virtual nonholonomic constraints, i.e.,
constraints that are linear on the velocities. Our definition is based on the
invariance property under the closed-loop system and coincides with that of
[9] in the linear case. In this paper we extend the results of [1] to the case
of affine constraints on the velocities.

In particular, a virtual affine nonholonomic constraint is described by an
affine non-integrable distribution on the configuration manifold of the system
for which there is a feedback control making it invariant under the flow of
the closed-loop system. We provide conditions for the existence and uniqueness
of such a feedback law defining the virtual affine nonholonomic constraint.

The remainder of the paper is structured as follows. Section 2 introduces
nonholonomic systems with affine constraints. We define virtual nonholonomic
constraints in Section 3, where we provide conditions for the existence and
uniqueness of a control law defining a virtual nonholonomic constraint. In
Section 4, we provide an example to illustrate the theoretical results.

##  2 Nonholonomic systems with affine constraints

We begin with a Lagrangian system on an $n$-dimensional configuration space
$Q$ and Lagrangian $L:TQ\to\mathbb{R}$. We assume the Lagrangian has the
mechanical form

$L=K-V\circ\pi,$ (1)

where $K$ is a function on $TQ$ describing the kinetic energy of the system,
that is, $K=\frac{1}{2}\mathscr{G}(v_{q},v_{q})$, where $\mathscr{G}$ is the
Riemannian metric on $Q$, $V:Q\to\mathbb{R}$ is a function on $Q$ representing
the potential energy, and $\pi:TQ\to Q$ is the tangent bundle projection,
locally given by $\pi(q,\dot{q})=q$ with $(q,\dot{q})$ denoting local
coordinates on $TQ$. In addition, we denote by $\mathfrak{X}(Q)$ the set of
vector fields on $Q$ and by $\Omega^{1}(Q)$ the set of $1$-forms on $Q$. If
$X,Y\in\mathfrak{X}(Q),$ then $[X,Y]$ denotes the standard Lie bracket of
vector fields.

Consider affine nonholonomic constraints, that is, for each $q\in Q$ the
velocities belong to an affine subspace $\mathscr{A}_{q}$ of the tangent space
$T_{q}Q.$ Thus, $\mathscr{A}_{q}$ can be written as a sum of a vector field
$X\in\mathfrak{X}(Q)$ and a nonintegrable distribution $\mathscr{D}$ on $Q$,
i.e. $\mathscr{A}_{q}=X(q)+\mathscr{D}_{q}$, where $\mathscr{D}$ is of
constant rank $r$, with $1<r<n$. In this case, we say that the affine space
$\mathscr{A}_{q}$ is modelled on the vector subspace $\mathscr{D}_{q}$. In
local coordinates $\mathscr{D}$ can be expressed as the null space of a
$q$-dependent matrix $S(q)$ of dimension $m\times n$ and
$\operatorname{rank}S(q)=m,$ with $m=n-r$ as $\mathscr{D}_{q}=\\{\dot{q}\in
T_{q}Q:S(q)\dot{q}=0\\}$. The rows of $S(q)$ can be represented by the
coordinate functions of $m$ independent 1-forms $\mu^{b}=\mu^{b}_{i}dq^{i}$,
$1\leq i\leq n,\;1\leq b\leq m$. The affine distribution is
$\mathscr{A}_{q}=\\{\dot{q}\in T_{q}Q:S(q)(\dot{q}-X(q))=0\\}$, hence
$\mathscr{A}=\\{(q,\dot{q})\in TQ:\phi(q,\dot{q})=0\\}$, with
$\phi(q,\dot{q})=S(q)\dot{q}+Z(q)$ and $Z(q)=-S(q)X(q)\in\mathbb{R}^{m}$ (see
[2], [3] for instance).

###### Definition 1

A mechanical system with affine nonholonomic constraints on a smooth manifold
$Q$ is given by the triple $(\mathscr{G},V,\mathscr{A})$, where $\mathscr{G}$
is a Riemannian metric on $Q,$ representing the kinetic energy of the system,
$V:Q\rightarrow\mathbb{R}$ is a smooth function representing the potential
energy, and $\mathscr{A}$ an affine distribution on $Q$ describing the affine
nonholonomic constraints.

In any Riemannian manifold, there is a unique connection
$\nabla^{\mathscr{G}}:\mathfrak{X}(Q)\times\mathfrak{X}(Q)\rightarrow\mathfrak{X}(Q)$
called the Levi-Civita connection satisfying the following two properties:

  1. 1.

$[X,Y]=\nabla_{X}^{\mathscr{G}}Y-\nabla_{Y}^{\mathscr{G}}X$ (symmetry)

  2. 2.

$X(\mathscr{G}(Y,Z))=\mathscr{G}(\nabla_{X}^{\mathscr{G}}(Y,Z)+\mathscr{G}(Y,\nabla_{X}^{\mathscr{G}}Z)$
(compatibillity of the metric).

The Levi-Civita connection helps us describe the trajectories of a mechanical
Lagrangian system. The trajectories $q:I\rightarrow Q$ of a mechanical
Lagrangian determined by a Lagrangian function as in (1) satisfy the following
equations

$\nabla_{\dot{q}}^{\mathscr{G}}\dot{q}+\text{grad}_{\mathscr{G}}V(q(t))=0,$
(2)

where the vector field $\text{grad}_{\mathscr{G}}V\in\mathfrak{X}(Q)$ is
characterized by $\mathscr{G}(\text{grad}_{\mathscr{G}}V,X)=dV(X),\;\mbox{ for
every }X\in\mathfrak{X}(Q)$. Observe that if the potential function vanishes,
then the trajectories of the mechanical system are just the geodesics with
respect to the connection $\nabla^{\mathscr{G}}$.

##  3 Virtual affine nonholonomic constraints

In this section we present a detailed construction of virtual affine
nonholonomic constraints. As will be clear from the definition of virtual
affine nonholonomic constraints their existence is essentially linked to a
controlled system, rather than to the affine distribution $\mathscr{A}$
defined by the constraints. As a consequence we give the necessary tools for
controlled mechanical systems.

Given a Riemannian metric $\mathscr{G}$ on $Q$, we can use its non-degeneracy
property to define the musical isomoprhism
$\flat:\mathfrak{X}(Q)\rightarrow\Omega^{1}(Q)$ defined by
$\flat(X)(Y)=\mathscr{G}(X,Y)$ for any $X,Y\in\mathfrak{X}(Q)$. Also, denote
by $\sharp:\Omega^{1}(Q)\rightarrow\mathfrak{X}(Q)$ the inverse musical
isomorphism, i.e., $\sharp=\flat^{-1}$.

Given an external force $F^{0}:TQ\rightarrow T^{*}Q$ and a control force
$F:TQ\times U\rightarrow T^{*}Q$ of the form

$F(q,\dot{q},u)=\sum_{a=1}^{m}u_{a}f^{a}(q)$ (3)

where $f^{a}\in\Omega^{1}(Q)$ with $m<n$, $U\subset\mathbb{R}^{m}$ the set of
controls and $u_{a}\in\mathbb{R}$ with $1\leq a\leq m$ the control inputs,
consider the associated mechanical control system of the form

$\nabla^{\mathscr{G}}_{\dot{q}(t)}\dot{q}(t)=Y^{0}(q(t),\dot{q}(t))+u_{a}(t)Y^{a}(q(t)),$
(4)

with $Y^{0}=\sharp(F^{0})$ and $Y^{a}=\sharp(f^{a})$ the corresponding force
vector fields.

Hence, $q$ is the trajectory of a vector field of the form

$\Gamma(v_{q})=G(v_{q})+u_{a}(Y^{a})_{v_{q}}^{V},$ (5)

where $G$ is the vector field determined by the unactuated forced mechanical
system

$\nabla^{\mathscr{G}}_{\dot{q}(t)}\dot{q}(t)=Y^{0}(q(t),\dot{q}(t))$

and where the vertical lift of a vector field $X\in\mathfrak{X}(Q)$ to $TQ$ is
defined by

$X_{v_{q}}^{V}=\left.\frac{d}{dt}\right|_{t=0}(v_{q}+tX(q)).$

###### Definition 2

The distribution $\mathcal{F}\subseteq TQ$ generated by the vector fields
$\sharp(f_{i})$ is called the input distribution associated with the
mechanical control system (4).

Now we will introduce the concept of virtual affine nonholonomic constraint.

###### Definition 3

A virtual affine nonholonomic constraint associated with the mechanical
control system (4) is a controlled invariant affine distribution
$\mathscr{A}\subseteq TQ$ for that system, that is, there exists a control
function $\hat{u}:\mathscr{A}\rightarrow\mathbb{R}^{m}$ such that the solution
of the closed-loop system satisfies
$\psi_{t}(\mathscr{A})\subseteq\mathscr{A}$, where $\psi_{t}:TQ\rightarrow TQ$
denotes its flow.

Before we proceed to the theorem which gives the necessary conditions for the
existence and uniqueness of a control law that turns an affine distribution
into a controlled invariant affine distribution (virtual affine nonholonomic
constraint), we present some necessary preliminaries.

###### Definition 4

If $W$ is an affine subspace of the vector space $V$ modelled on the vector
subspace $W_{0}$, then the dimension of the affine subspace $W$ is defined to
be the dimension of the model vector subspace $W_{0}$.

Two affine subspaces $W_{1}$ and $W_{2}$ of a vector space $V$ are transversal
and we write $W_{1}\pitchfork W_{2}$ if

  1. 1.

$V=W_{1}+W_{2}$.

  2. 2.

$\dim V=\dim W_{1}+\dim W_{2}$, i.e., the dimensions of $W_{1}$ and $W_{2}$
are complementary with respect to the ambient space dimension.

###### Remark 1

If $W_{1}$ and $W_{2}$ are subspaces of $V$ then the previous definition
implies that $V=W_{1}\oplus W_{2}$.

###### Remark 2

If $W_{1}$ and $W_{2}$ are affine spaces modelled on vector subsspaces
$W_{10}$ and $W_{20}$, respectively, then $W_{1}\pitchfork W_{2}$ if and only
if $V=W_{10}\oplus W_{20}$.

###### Proposition 1

For two distributions $\mathscr{A}$ and $\mathscr{F}$ on a manifold $Q$ where
$\mathscr{A}$ is an affine distribution with $\mathscr{D}$ the associated
model distribution and $X$ a vector field on $Q$ satisfying
$\mathscr{A}=X+\mathscr{D}$, the tranversality condition for $\mathscr{A}$ and
$\mathscr{F}$ is an inherited property from the transversality of the model
distribution $\mathscr{D}$ and vice versa, namely,
$\mathscr{A}\pitchfork\mathscr{F}\Leftrightarrow\mathscr{D}\pitchfork\mathscr{F}$.

###### Proof

First consider that $\mathscr{A}\pitchfork\mathscr{F}$ which means that for
every $q\in Q$ we have

$T_{q}Q=\mathscr{A}_{q}+\mathscr{F}_{q}=X(q)+\mathscr{D}_{q}+\mathscr{F}_{q}$

hence, for every $v_{q}\in T_{q}Q$ there exist vectors
$d_{q}\in\mathscr{D}_{q}$ and $f_{q}\in\mathscr{F}_{q}$ such that

$v_{q}=X(q)+d_{q}+f_{q}\;\Leftrightarrow\;v_{q}-X(q)=d_{q}+f_{q}.$

Since $v_{q}-X(q)\in T_{q}Q$ and $v_{q}$ is arbitrary we have
$T_{q}Q=\mathscr{D}_{q}+\mathscr{F}_{q}$ for every $q\in Q$. Together with the
fact that $\dim\mathscr{D}_{q}+\dim\mathscr{F}_{q}=\dim T_{q}Q$, we have that
$\mathscr{D}\pitchfork\mathscr{F}.$
Now suppose that $\mathscr{D}\pitchfork\mathscr{F}$. Note that this is the
same as $TQ=\mathscr{D}\oplus\mathscr{F}$. Hence, as before, for $v_{q}\in
T_{q}Q$ there are $d_{q}\in\mathscr{D}_{q}$ and $f_{q}\in\mathscr{F}_{q}$ such
that

$v_{q}=d_{q}+f_{q}\;\Leftrightarrow\;v_{q}+X(q)=X(q)+d_{q}+f_{q}.$

By the same argument as above, $v_{q}+X(q)\in T_{q}Q$ and $v_{q}$ is
arbitrary, thus, together with the dimension condition, we conclude that
$\mathscr{A}\pitchfork\mathscr{F}.$

###### Proposition 2

Consider two distributions $\mathscr{A}$ and $\mathscr{F}$ where the first is
an affine distribution as defined previously i.e. $\mathscr{A}=X+\mathscr{D}$,
with $X\in\mathfrak{X}(Q)$ and $\mathscr{D}$ its associated model
distribution. For $v_{q}\in\mathscr{A}$ we have

$\mathscr{A}\pitchfork\mathscr{F}\Longrightarrow
T_{v_{q}}(TQ)=T_{v_{q}}\mathscr{A}\oplus\mathscr{F}^{V}_{v_{q}},$

where $\mathscr{F}^{V}_{v_{q}}$ is the vertical lift of $\mathscr{F}_{v_{q}}.$

###### Proof

From the structure of $\mathscr{A}$, i.e., from the fact that each
$v_{q}\in\mathscr{A}_{q}$ can be written as $v_{q}=Z(q)+d_{q}$ where
$d_{q}\in\mathscr{D}_{q}$, we may conclude that $\mathscr{A}_{q}$ is a
$r$-dimensional manifold, where $r$ is the rank of the distribution
$\mathscr{D}$. Thus $\mathscr{A}$ is a fiber bundle whose base space is the
$n$ dimensional manifold $Q$ and whose fibers are $r$ dimensional affine
subspaces. Hence,

$\dim(T_{v_{q}}\mathscr{A})=\dim(T_{d_{q}}\mathscr{D})=n+r$

and since $\dim\mathscr{F}^{V}_{v_{q}}=n-r=m$, we have that

$\dim T_{v_{q}}(TQ)=\dim
T_{v_{q}}\mathscr{A}+\dim\mathscr{F}^{V}_{v_{q}}=n+r+m=2n.$

So, in order to prove that both subspaces are transversal it suffices to prove
that their intersection contains only the zero tangent vector. Indeed, suppose
that $v_{q}\in\mathscr{A}$ and $X_{v_{q}}\in T_{v_{q}}\mathscr{A}$. Since
$\mathscr{A}$ is defined to be the set of vectors satisfying the equation
$\phi=0$, the tangent vector satisfies $T_{v_{q}}\phi(X_{v_{q}})=0$. If, in
addition, $X_{v_{q}}\in\mathscr{F}^{V}_{v_{q}}$, then it can be written as

$X_{v_{q}}=c^{i}\sharp(f_{i})^{V}_{v_{q}}.$

However,

$T_{v_{q}}\phi(\sharp(f_{i})^{V}_{v_{q}})=(S(q)\sharp(f_{i}))^{V}_{v_{q}}$

from where it follows that if $c^{i}\sharp(f_{i})^{V}_{v_{q}}$ was in the null
space of the linear map $T_{v_{q}}\phi$, then $c^{i}\sharp(f_{i})$ would be in
the null space of $S(q)$ which is false, since these are vectors in
$\mathscr{D}_{q}$ and $\mathscr{F}$ and $\mathscr{D}$ are transversal using
the previous proposition. Thus $X_{v_{q}}=0$.

###### Theorem 3.1

If the affine distribution $\mathscr{A}$ and the control input distribution
$\mathscr{F}$ are transversal, then there exists a unique control function
making the distribution a virtual affine nonholonomic constraint associated
with the mechanical control system (4).

###### Proof

Suppose that $\mathscr{A}\pitchfork\mathscr{F}$ and that trajectories of the
contol system (4) may be written as the integral curves of the vector field
$\Gamma$ defined by (5). From Proposition 2 we have

$T_{v_{q}}(TQ)=T_{v_{q}}\mathscr{A}\oplus\mathscr{F}^{V}_{v_{q}},$

where $v_{q}\in\mathscr{A}$ and
$\mathscr{F}^{V}_{v_{q}}=\operatorname{span}\\{(Y^{a})^{V}_{v_{q}}\\}$. Using
the uniqueness decomposition property arising from transversality, we conclude
there exists a unique vector
$\tau^{*}(v_{q})=(\tau_{1}^{*}(v_{q}),\cdots,\tau_{m}^{*}(v_{q}))\in\mathbb{R}^{m}$
such that $\Gamma(v_{q})=G(v_{q})+\tau_{a}^{*}(v_{q})(Y^{a})_{v_{q}}^{V}\in
T_{v_{q}}\mathscr{A}.$ Next, we show that $\Gamma$ depends smoothly on
$v_{q}.$ If $\mathscr{A}$ is defined by $m$ constraints of the form
$\phi^{b}(v_{q})=0$, $1\leq b\leq m$, then the condition above may be
rewritten as $d\phi^{b}(G(v_{q})+\tau_{a}^{*}(v_{q})(Y^{a})_{v_{q}}^{V})=0,$
which is equivalent to

$\tau_{a}^{*}(v_{q})d\phi^{b}((Y^{a})_{v_{q}}^{V})=-d\phi^{b}(G(v_{q})).$

Note that, the equation above is a linear equation of the form
$P(v_{q})\tau=b(v_{q})$, where $b(v_{q})$ is the vector
$(-d\phi^{1}(\Gamma(v_{q})),\dots,-d\phi^{m}(\Gamma(v_{q})))\in\mathbb{R}^{m}$
and $P(v_{q})$ is the $m\times m$ matrix with entries
$P^{b}_{a}(v_{q})=d\phi^{b}((Y^{a})_{v_{q}}^{V})=\mu^{b}(q)(Y^{a})$, where the
last equality may be deduced by computing the expressions in local
coordinates. That is, if $(q^{i}\dot{q}^{i})$ are natural bundle coordinates
for the tangent bundle, then

$\begin{split}d\phi^{b}((Y^{a})_{v_{q}}^{V})&=\left(\frac{\partial\mu^{b}_{i}}{\partial
q^{j}}\dot{q}^{i}dq^{j}+\frac{\partial Z_{i}}{\partial
q^{j}}dq^{j}+\mu^{b}_{i}d\dot{q}^{i}\right)\left(Y^{a,k}\frac{\partial}{\partial\dot{q}^{k}}\right)\\\
&=\mu^{b}_{i}Y^{a,i}=\mu^{b}(q)(Y^{a}).\end{split}$

In addition, $P(v_{q})$ has full rank, since its columns are linearly
independent. In fact suppose that

$c_{1}\begin{bmatrix}\mu^{1}(Y^{1})\\\ \vdots\\\
\mu^{m}(Y^{1})\end{bmatrix}+\cdots+c_{m}\begin{bmatrix}\mu^{1}(Y^{m})\\\
\vdots\\\ \mu^{m}(Y^{m})\end{bmatrix}=0,$

which is equivalent to

$\begin{bmatrix}\mu^{1}(c_{1}Y^{1}+\cdots+c_{m}Y^{m})\\\ \vdots\\\
\mu^{m}(c_{1}Y^{1}+\cdots+c_{m}Y^{m})\end{bmatrix}.$

However, from Proposition 1 we have $\mathscr{D}\cap\mathscr{F}=\\{0\\}$ which
implies that $c_{1}Y^{1}+\cdots+c_{m}Y^{m}=0$. Since $\\{Y_{i}\\}$ are
linearly independent we conclude that $c_{1}=\cdots=c_{m}=0$ and $P$ has full
rank. But, since $P$ is an $m\times m$ matrix, and $\mathscr{D}$ is a regular
distribution, it must be invertible. Therefore, there is a unique vector
$\tau^{*}(v_{q})$ satisfying the matrix equation and
$\tau^{*}:\mathscr{D}\rightarrow\mathbb{R}^{m}$ is smooth since it is the
solution of a matrix equation depending smoothly on $v_{q}$. Hence, $\Gamma$
is a smooth vector field tangent to $\mathscr{A}$ and its flow remains in
$\mathscr{A}.$ $\square$

##  4 An example

Consider a boat with a payload on the see with a position-dependent stream.
The position of the boat’s center of mass is modeled by the configuration
manifold $\mathbb{R}^{2}$ to which we add an orientation to obtain a complete
description of its location in space, so that the system total configuration
manifold is $\mathbb{R}^{2}\times\mathbb{S}$ with local coordinates
$q=(x,y,\theta)$. The sea’s current is modeled by the vector field
$C:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$, $C=(C^{1}(x,y),C^{2}(x,y))$.

The boat is well modeled by a forced mechanical system with Lagrangian
function $L=\frac{m}{2}(\dot{x}^{2}+\dot{y}^{2})+\frac{I}{2}\dot{\theta}^{2},$
where $m$ is the boat’s mass, $I$ is the moment of inertia, and the external
force is denoted by $F^{ext}=W^{1}dx+W^{2}dy$ accounting for the action of the
current on the center of mass of the boat and to which we add a control force
$F=u(\sin\theta dx-\cos\theta dy+d\theta).$

The functions $W^{1}$ and $W^{2}$ are defined according to

$\begin{cases}W^{1}&=m\ d\left(\sin^{2}\theta C^{1}-\sin\theta\cos\theta
C^{2}\right)(\dot{q})\\\ W^{2}&=m\ d\left(-\sin\theta\cos\theta
C^{1}+\cos^{2}\theta C^{2}\right)(\dot{q}).\end{cases},$

where $d$ represents the differential of the functions inside the parenthesis.
The external force assures that in the absence of controls, the dynamics of
the boat satisfies the following kinematic equations

$\begin{cases}\dot{x}=&\sin^{2}\theta C^{1}-\sin\theta\cos\theta C^{2}\\\
\dot{y}=&-\sin\theta\cos\theta C^{1}+\cos^{2}\theta C^{2},\end{cases}$

whenever the initial velocities in the $x$ and $y$ direction vanish. The
corresponding controlled forced Lagrangian system is

$m\ddot{x}=u\sin\theta+W^{1},\quad m\ddot{y}=-u\cos\theta+W^{2},\quad
I\ddot{\theta}=u,$

and, as we will show, it has the following virtual affine nonholonomic
constraint

$\sin\theta\dot{x}-\cos\theta\dot{y}=C^{2}\cos\theta-C^{1}\sin\theta.$

The input distribution $\mathcal{F}$ is generated just by one vector field

$Y=\frac{\sin\theta}{m}\frac{\partial}{\partial
x}-\frac{\cos\theta}{m}\frac{\partial}{\partial
y}+\frac{1}{I}\frac{\partial}{\partial\theta},$

while the virtual nonholonomic constraint is the affine space $\mathcal{A}$
modelled on the distribution $\mathcal{D}$ defined as the set of tangent
vectors $v_{q}\in T_{q}Q$ where $\mu(q)(v)=0,$ with $\mu=\sin\theta dx-
cos\theta dy$. Thus, we may write it as

$\mathcal{D}=\hbox{span}\Big{\\{}X_{1}=\cos\theta\frac{\partial}{\partial
x}+\sin\theta\frac{\partial}{\partial
y},\,X_{2}=\frac{\partial}{\partial\theta}\Big{\\}}.$

The affine space is given as the zero set of the function
$\phi(q,v)=\mu(q)(v)+Z(q)$ with $Z(q)=\cos\theta C^{2}(x,y)-\sin\theta
C^{1}(x,y)$ or, equivalently, as the set of vectors $v_{q}$ satisfying
$v_{q}-C(q)\in\mathcal{D}_{q}$.

We may check that $\mathcal{A}$ is controlled invariant for the controlled
Lagrangian system above. In fact, the control law
$\hat{u}(x,y,\theta,\dot{x},\dot{y},\dot{\theta})=-m\dot{\theta}(\cos\theta\dot{x}+\sin\theta\dot{y})$
makes the affine space invariant under the closed-loop system, since in this
case, the dynamical vector field arising from the controlled Euler-Lagrange
equations given by

$\Gamma=\dot{x}\frac{\partial}{\partial x}+\dot{y}\frac{\partial}{\partial
y}+\dot{\theta}\frac{\partial}{\partial\theta}+\left(\frac{\hat{u}\sin\theta+W^{1}}{m}\right)\frac{\partial}{\partial\dot{x}}+\left(-\frac{\hat{u}\cos\theta-W^{2}}{m}\right)\frac{\partial}{\partial\dot{y}}+\frac{\hat{u}}{I}\frac{\partial}{\partial\dot{\theta}}$

is tangent to $\mathcal{A}$. This is deduced from the fact that

$\Gamma(\sin\theta\dot{x}-\cos\theta\dot{y}+\cos\theta C^{2}(x,y)-\sin\theta
C^{1}(x,y))=0.$

## References

  * [1] A. Anahory Simoes, E. Stratoglou, A. Bloch, L. Colombo. Virtual Nonholonomic Constraints: A Geometric Approach. arXiv e-prints, arXiv: 2207.01299 (2022)
  * [2] F. Fasso, N. Sansonetto. Conservation of energy and momenta in nonholonomic systems with affine constraints. Regular and Chaotic Dynamics, 20(4), 449-462 (2015).
  * [3] F. Fasso, L. García-Naranjo, N. Sansonetto. Moving energies as first integrals of nonholonomic systems with affine constraints. Nonlinearity, 31(3), 755 (2018).
  * [4] B. Griffin, J. Grizzle Nonholonomic virtual constraints for dynamic walking. 54th IEEE Conference on Decision and Control 4053–4060 (2015).
  * [5] K. Hamed, A. Ames. Nonholonomic hybrid zero dynamics for the stabilization of periodic orbits: Application to underactuated robotic walking. IEEE Transactions on Control Systems Technology, 28(6), 2689–2696 (2019).
  * [6] J. Horn, A. Mohammadi, K. Hamed, R. Gregg. Nonholonomic virtual constraint design for variable-incline bipedal robotic walking. IEEE Robotics and Automation Letters, 5(2), 3691–3698 (2020).
  * [7] J. Horn, A. Mohammadi, K. Hamed, R. Gregg. Hybrid zero dynamics of bipedal robots under nonholonomic virtual constraints. IEEE Control Systems Letters, 3(2), 386–391 (2018).
  * [8] J. Horn, R.Gregg. Nonholonomic Virtual Constraints for Control of Powered Prostheses Across Walking Speeds. IEEE Transactions on Control Systems Technology. (2021).
  * [9] A. Moran-MacDonald. Energy injection for mechanical systems through the method of Virtual Nonholonomic Constraints. University of Toronto (2021).