# Finslerian convolution metrics and their special classes

Gilbert Nibaruta  Ecole Normale Supérieure
P. O. Box 6983
Bujumbura
Burundi<EMAIL_ADDRESS>
###### Abstract.

Here, it is introduced a concept of convolution metric in Finslerian Geometry.
This convolution metric is a kind of function obtained by a given mathematical
operation between two Finslerian metrics. Some basic properties of the
Finslerian convolution metrics are studied. Then it is characterized
Finslerian convolution metrics which are of type Riemannian, Minkowskian as
well as Randers. Furthermore, some examples of the Finslerian convolutions are
given.

###### Key words and phrases:

Finslerian metric; Riemannian metric; Randers metric; convolution metric.

###### 2010 Mathematics Subject Classification:

primary 53B40, secondary 58B20

##  1\. Introduction

Convolution metrics play an important role in Differential Geometry and have
applications that include General Relativity since they are natural extensions
of warped product metrics. As a concrete example of application, many basic
solutions of the Einstein field equations on Riemannian manifolds have
solutions when the Riemannian metrics are warped products [4] at most.

Nowadays, a Finslerian metric is of considerable interest due to its rich
structure including mainly Riemann, Randers, Minkowski and Berwald type
metrics. For the details, the study of each of these special Finslerian
metrics and some important examples can be found in [1, 5] and in [6].

In Riemannian Geometry, the notion of convolution metrics was introduced by
Bang-Yen Chen to establish and characterize the Euclidean version of Segre
embedding [2]. The following Proposition shows that the notion of convolution
metrics arises very naturally from tensor product immersions.

######  Proposition 1.1.

[3] Let $M_{1}$ and $M_{2}$ be two $C^{\infty}$ manifolds of dimensions
$n_{1}$ and $n_{2}$ respectively. Consider $g_{1}$ and $g_{2}$ two Riemannian
metrics on $M_{1}$ and $M_{2}$ respectively. If
$\varphi_{1}:(M_{1},g_{1})\longrightarrow\mathbb{R}^{n_{1}}\backslash\\{0\\}$
and
$\varphi_{2}:(M_{2},g_{2})\longrightarrow\mathbb{R}^{n_{2}}\backslash\\{0\\}$
are isometric immersions then the mapping

$\displaystyle\varphi=\varphi_{1}\otimes\varphi_{2}:M_{1}\times
M_{2}\longrightarrow\mathbb{R}^{n_{1}}\otimes\mathbb{R}^{n_{2}}=\mathbb{R}^{n_{1}n_{2}}$

defined by $\varphi(x_{1},x_{2})=\varphi_{1}(x_{1})\varphi_{2}(x_{2})$ gives
rise to a convolution manifold $M_{1}\bigstar M_{2}$ endowed with the
convolution metric

${}_{\rho_{2}}g_{1}\star_{\rho_{1}}g_{2}=\rho_{2}^{2}g_{1}+\rho_{1}^{2}g_{2}+2\rho_{1}\rho_{2}d\rho_{1}\otimes
d\rho_{2}$

whenever ${}_{\rho_{2}}g_{1}\star_{\rho_{1}}g_{2}$ is non-degenerate where
$\rho_{1}=\|\varphi_{1}\|$ and $\rho_{2}=\|\varphi_{2}\|$.

In this work, we construct an extension of the Riemannian convolution metric
on a product manifold (See the Definition 2.5). More precisely, we define a
new notion of convolution metric by taking account of, not only the Riemannian
metrics but also the Minkowski, the Randers and the Berwald metrics. We have
the following:

######  Proposition 1.2.

Let $(M_{1},F_{1})$ and $(M_{2},F_{2})$ be two Finslerian manifolds. Consider
$f_{k}:M_{k}\longrightarrow(0,\infty)$, for $k\in\\{1,2\\}$, some smooth
functions. If ${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ is a convolution Finslerian
metric on a manifold $M_{1}\times M_{2}$ then one of the following assertions
holds

  * (1)

the gradient of one of the functions $f_{k}$ is orthogonal to any tangent
vector on $M_{1}\times M_{2}$.

  * (2)

one of the $f_{k}$ is a nonzero constant. In this last case, when one of the
$f_{k}\equiv 1$, the convolution Finslerian metric
${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ reduces to a Finslerian warped product
metric.

Main results regard the characterization of convolution metrics of type
Riemannian, Minkowskian as well as Randers among Finslerian ones. For example,
we prove the following.

######  Theorem 1.3.

Let $(M_{1},F_{1})$ and $(M_{2},F_{2})$ be two Finslerian manifolds of Randers
type that is, for $k\in\\{1,2\\}$, $F_{k}:=\alpha_{k}+\beta_{k}$ where
$\alpha_{k}$ is a Riemannian metric and $\beta_{k}$ is a $1$-form on $M_{k}$
with $||\beta_{k}||_{\alpha_{k}}<1$. If one of the
$f_{k}:M_{k}\longrightarrow(0,\infty)$ is a nonzero constant function, then
the Finslerian convolution metric ${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ on
$M_{1}\times M_{2}$ is a Randers metric if and only if

$\displaystyle\frac{\alpha_{1}}{\alpha_{2}}=\frac{\beta_{1}}{\beta_{2}}.$
(1.1)

One of the examples of special Finslerian convolution is given by

######  Example 1.4.

Consider
$M=\\{x=(x^{1},x^{2},x^{3},x^{4})\in\mathbb{R}^{4}:x_{1}>0,x_{3}>0\\}$ and
$T_{x}M=\\{(y^{1},y^{2},y^{3},y^{4})\in\mathbb{R}^{4}:y_{1}>0,y_{3}>0\\}$. For
$\lambda\in[2,4]$ and for a positive integer $k$, define
$F:TM\longrightarrow(0,\infty)$ by

$\displaystyle F(x^{1},x^{2},x^{3},x^{4},y^{1},y^{2},y^{3},y^{4})$
$\displaystyle=$ $\displaystyle\left\\{x_{3}^{2}\left(y_{1}^{4}+\lambda
y_{1}^{2}y_{2}^{2}+y_{2}^{4}\right)^{\frac{1}{2}}+8x_{1}^{3}x_{3}^{3}y_{1}y_{3}\right.$
(1.2)
$\displaystyle+\left.x_{3}^{2}\left[y_{3}^{2}+y_{4}^{2}+\lambda\left(y_{3}^{2k}+y_{4}^{2k}\right)^{\frac{1}{k}}\right]\right\\}^{\frac{1}{2}}.$

Then the function $F$ defined in (1.2) is a Finslerian convolution metric.

The outline of this document is organized as follows. In Section 2, we
establish some basic notions on Finslerian convolution manifolds. The Section
3 is devoted to the characterizations of the convolution of Riemannian,
Minkowskian and Randers metrics among the Finslerian convolution metrics. In
the Section 4 we provide examples of Finslerian convolution metrics.

Notation. In the following, $M$ is an $n-$dimensional product of the
$C^{2}$-manifolds $M_{1}$ and $M_{2}$ of dimensions $n_{1}$ and $n_{2}$
respectively. We denote by $T_{x}M$ the tangent space at $x=(x_{1},x_{2})\in
M$ and by $TM:=\bigcup_{x\in M}T_{x}M$ the tangent bundle of $M$. Set
$\mathring{T}M=\\{(x_{1},x_{2},y_{1},y_{2})\in TM:y_{1}\neq 0,y_{2}\neq
0\text{ where }0\text{ is a zero vector }\in T_{x_{k}}M_{k}\text{ for
}k\in\\{1,2\\}\\}$ and $\pi:TM\longrightarrow M:\pi(x,y)\longmapsto x$ the
natural projection. Let
$(x^{1},...,x^{n})=(x_{1}^{1},...,x_{1}^{n_{1}},x_{2}^{n_{1}+1},...,x_{2}^{n})$
be a local coordinate on an open subset $U$ of $M$ and
$(x^{1},...,x^{n},y^{1},...,y^{n})$ be the local coordinate on
$\pi^{-1}(U)\subset TM$. The local coordinate system $(x^{i})_{i=1,...,n}$
produces the coordinate bases $\\{\frac{\partial}{\partial
x^{i}}\\}_{i=1,...,n}$ and $\\{dx^{i}\\}_{i=1,...,n}$ respectively, for $TM$
and cotangent bundle $T^{*}M$.

Convention. We use Einstein summation convention: repeated upper and lower
indices will automatically be summed, from $1$ to $n$, unless otherwise will
be noted. For example $g_{ij}v^{i}v^{j}=\sum_{i,j=1}^{n}g_{ij}v^{i}v^{j}.$

##  2\. Finslerian convolution manifold

######  Definition 2.1.

Let $M$ be an $n$-dimensional manifold. A function
$F:TM\longrightarrow[0,\infty)$ is a Finslerian metric on $M$ if :

  * (i)

$F$ is $C^{\infty}$ on the entire slit tangent bundle $\mathring{T}M$,

  * (ii)

$F$ is positively $1$-homogeneous on the fibers of $TM$, that is
$\forall c>0,~{}F(x,cy)=cF(x,y),$

  * (iii)

the Hessian matrix $(g_{ij}(x,y))_{1\leq i,j\leq n}$ with elements

$\displaystyle g_{ij}(x,y):=\frac{1}{2}\frac{\partial^{2}F^{2}(x,y)}{\partial
y^{i}\partial y^{j}}$ (2.1)

is positive-definite at every point $(x,y)$ of $\mathring{T}M$.

######  Remark 2.2.

The Hessian matrix $g$ whose elements are defined in (2.1) is a natural
Riemannian metric on the pulled-back bundle $\pi^{*}TM$ over the manifold
$\mathring{T}M$. $\pi^{*}TM$ is a vector bundle and its fiber at a point
$(x,y)\in\mathring{T}M$ is

$\displaystyle\pi^{*}TM|_{(x,y)}:=\\{(x,y;v):v\in T_{x}M\\}|_{(x,y)}\cong
T_{x}M.$ (2.2)

######  Definition 2.3.

Let $(M,F)$ be an $n$-dimensional Finslerian manifold and $g$ the Hessian
matrix associated with $F$. The gradient of a smooth function $u$ on $M$ is,
the section of the vector bundle $\pi^{*}TM$ denoted by $\nabla u$, given by

$\displaystyle g_{(x,y)}(\nabla
u_{(x,y)},\xi_{(x,y)})=du_{(x,y)}(\xi_{(x,y)})$ (2.3)

for any $\xi\in\Gamma(\pi^{*}TM)$ and for every $(x,y)\in\mathring{T}M$.
Locally,

$\displaystyle\nabla u_{(x,y)}=g^{ij}(x,y)\frac{\partial u}{\partial
x^{i}}\frac{\partial}{\partial x^{j}}.$ (2.4)

######  Remark 2.4.

Let $M_{1}$ and $M_{2}$ be two $C^{\infty}$ manifolds. For every
$(x_{1},x_{2})\in M_{1}\times M_{2}$, we have the following properties derived
from $M_{1}$ and $M_{2}$.

  * (1)

The projections

$\displaystyle\sigma_{1}$ $\displaystyle:$ $\displaystyle M_{1}\times
M_{2}\longrightarrow M_{1}\text{ such that }\sigma_{1}(x_{1},x_{2})=x_{1}$
$\displaystyle\sigma_{2}$ $\displaystyle:$ $\displaystyle M_{1}\times
M_{2}\longrightarrow M_{1}\text{ such that }\sigma_{2}(x_{1},x_{2})=x_{2}$

are $C^{\infty}$ submersions.

  * (2)

$dim(M_{1}\times M_{2})=dimM_{1}+dimM_{2}$.

######  Definition 2.5.

Let $(M_{1},F_{1})$ and $(M_{2},F_{2})$ be two Finslerian manifolds. Consider
$f_{k}:M_{k}\longrightarrow(0,\infty)$, for $k\in\\{1;2\\}$, some smooth
functions. A convolution of $(M_{1},F_{1})$ and $(M_{2},F_{2})$, via $f_{1}$
and $f_{2}$, is the product manifold $M_{1}\times M_{2}$ endowed with a
function $F_{f_{1},f_{2}}:T(M_{1}\times M_{2})\longrightarrow(0,\infty)$
defined by

$\displaystyle
F_{f_{1},f_{2}}=\sqrt{h_{2}^{2}(\sigma_{1}^{*}F_{1})^{2}+h_{1}^{2}(\sigma_{2}^{*}F_{2})^{2}+\frac{1}{2}h_{1}h_{2}\frac{\partial(\sigma_{1}^{*}F_{1})^{2}}{\partial
y^{i}}\frac{\partial(\sigma_{2}^{*}F_{2})^{2}}{\partial y^{j}}\big{(}\nabla
h_{1}\big{)}^{i}\big{(}\nabla h_{2}\big{)}^{j}}$ (2.5)

where $h_{k}:=f_{k}\circ\sigma_{k}$. The function $F_{f_{1},f_{2}}$ is called
a Finslerian convolution of $F_{1}$ and $F_{2}$, via $f_{1}$ and $f_{2}$, and
is denoted by ${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$. Given a Finslerian
convolution function ${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ on $M_{1}\times
M_{2}$, the pair $(M_{1}\times M_{2},_{f_{2}}F_{1}\star_{f_{1}}F_{2})$ is
called a Finslerian convolution manifold and is denoted by
${}_{f_{2}}M_{1}\bigstar_{f_{1}}M_{2}$.

######  Remark 2.6.

For any $x\equiv(x_{1},x_{2})\in M\equiv M_{1}\times M_{2}$ and
$y\equiv(y_{1},y_{2})\in T_{x}M$,

  * (1)

$F_{f_{1},f_{2}}(x,y)$ exists for some positive-real functions $f_{1}$ and
$f_{2}$ on $M_{1}$ and $M_{2}$ respectively.

  * (2)

$F_{f_{1},f_{2}}$ is not $C^{\infty}$ on the tangent vectors of the form
$(y_{1},0)$ nor $(0,y_{2})$ at a point $(x_{1},x_{2})\in M_{1}\times M_{2}$
but on $(y_{1},y_{2})$ for $y_{1}\neq 0$ and $y_{2}\neq 0$.

  * (3)

${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ is positively homogeneous of degree one
with respect to $(y_{1},y_{2})$.

  * (4)

from the relation (2.5), we have

$\displaystyle F^{2}_{f_{1},f_{2}}(x,y)$ $\displaystyle=$ $\displaystyle
2f_{1}(x_{1})f_{2}(x_{2})F_{1}(x_{1},y_{1})F_{2}(x_{2},y_{2})\frac{\partial
F_{1}(x_{1},y_{1})}{\partial y^{i}}\frac{\partial F_{2}(x_{2},y_{2})}{\partial
y^{j}}\big{(}\nabla f_{1}(x_{1})\big{)}^{i}\big{(}\nabla
f_{2}(x_{2})\big{)}^{j}$ (2.6)
$\displaystyle+f^{2}_{2}(x_{2})F^{2}_{1}(x_{1},y_{1})+f_{1}(x_{1})^{2}F^{2}_{2}(x_{2},y_{2}).$

Since $(x_{1},x_{2})$ is arbitrary a point of $M$, the expression (2.6) can be
written as

$\displaystyle
F^{2}_{f_{1},f_{2}}=f_{2}^{2}F_{1}^{2}+f_{1}^{2}F_{2}^{2}+2f_{1}f_{2}F_{1}F_{2}\frac{\partial
F_{1}}{\partial y^{i}}\frac{\partial F_{2}}{\partial y^{j}}\big{(}\nabla
f_{1}\big{)}^{i}\big{(}\nabla f_{2}\big{)}^{j}.$ (2.7)

######  Lemma 2.7.

Let $M_{1}$ and $M_{2}$ be two $C^{\infty}$ manifolds of dimensions $n_{1}$
and $n_{2}$, respectively, with $n_{1}+n_{2}=n$. If
$F=_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ is a Finslerian convolution function on
the manifold $M_{1}\times M_{2}$ then the Hessian matrix associated with $F$
is given by

$\displaystyle\Big{(}g_{ij}(x_{1},x_{2},y_{1},y_{2})\Big{)}=\left(\begin{array}[]{cc}f_{2}^{2}(x_{2})\Big{(}g_{1\alpha\beta}(x_{1},y_{1})\Big{)}&2f_{1}(x_{1})f_{2}(x_{2})\Big{(}\frac{\partial
f_{1}(x_{1})}{\partial x_{1}^{\alpha}}\frac{\partial f_{2}(x_{2})}{\partial
x_{2}^{\nu}}\Big{)}\\\
\Big{(}0\Big{)}&f_{1}^{2}(x_{1})\Big{(}g_{2\mu\nu}(x_{2},y_{2})\Big{)}\end{array}\right)$
(2.10)

at every point $(x_{1},x_{2},y_{1},y_{2})$ of
$\mathring{T}M_{1}\times\mathring{T}M_{2}$ where
$i,j\in[1,...,n],\alpha,\beta\in[1,...,n_{1}],\mu,\nu\in[n_{1}+1,...,n]$ and,
$g_{1}$ and $g_{2}$ are the Hessian matrices of $F_{1}$ and $F_{2}$
respectively.

Proof. The proof follows from a straightforward calculation using the
Definition 2.1 and the relation (2.7).

######  Remark 2.8.

In general, the components $g_{ij}$ of the convolution matrix $g$ are not
symmetric in all indices.

######  Proposition 2.9.

Let $M_{1}$ and $M_{2}$ be two $C^{\infty}$ manifolds of dimensions $n_{1}$
and $n_{2}$, respectively, with $n_{1}+n_{2}=n$. If
${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ is a Finslerian convolution metric on a
manifold $M_{1}\times M_{2}$ and if $v\neq 0$ is a tangent vector on
$M_{1}\times M_{2}$ then the quadratic form
$(_{f_{2}}g_{1}\star_{f_{1}}g_{2})(v,v)$ associated with $F$ is positive-
definite if

$\displaystyle\frac{1}{f_{1}^{2}}g_{1\alpha\beta}v^{\alpha}v^{\beta}+\frac{1}{f_{2}^{2}}g_{2}v^{\mu}v^{\nu}>-\frac{2}{f_{1}f_{2}}\frac{\partial
f_{1}}{\partial x_{1}^{\alpha}}\frac{\partial f_{2}}{\partial
x_{2}^{\nu}}v^{\alpha}v^{\nu}$ (2.11)

where $\alpha,\beta\in[1,...,n_{1}],\mu,\nu\in[n_{1}+1,...,n]$ and, $g_{1}$
and $g_{2}$ are the Hessian matrices of $F_{1}$ and $F_{2}$ respectively.

Proof. From the Lemma 2.7, for any tangent vector $v\neq 0$ on $M_{1}\times
M_{2}$, we have

$\displaystyle g(v,v)$ $\displaystyle=$
$\displaystyle\sum_{i,j=1}^{n}g_{ij}v^{i}v^{j}$ $\displaystyle=$
$\displaystyle\sum_{\alpha,\beta=1}^{n_{1}}g_{\alpha\beta}v^{\alpha}v^{\beta}+\sum_{\alpha=1}^{n_{1}}\sum_{\mu=n_{1}+1}^{n}g_{\alpha\mu}v^{\alpha}v^{\mu}+\sum_{\mu=n_{1}+1}^{n}\sum_{\alpha=1}^{n_{1}}g_{\mu\alpha}v^{\mu}v^{\alpha}+\sum_{\mu,\nu=n_{1}+1}^{n}g_{\mu\nu}v^{\mu}v^{\nu}.$

Since $g\equiv_{f_{2}}g_{1}\star_{f_{1}}g_{2}$, it follows that

$\displaystyle g_{ij}(v,v)$ $\displaystyle=$
$\displaystyle\sum_{\alpha,\beta=1}^{n_{1}}f_{2}^{2}g_{1\alpha\beta}v^{\alpha}v^{\beta}+\sum_{\alpha=1}^{n_{1}}\sum_{\mu=n_{1}+1}^{n}2f_{1}f_{2}\frac{\partial
f_{1}}{\partial x_{1}^{\alpha}}\frac{\partial f_{2}}{\partial
x_{2}^{\mu}}v^{\alpha}v^{\mu}+0+\sum_{\mu,\nu=n_{1}+1}^{n}f_{1}^{2}g_{\mu\nu}v^{\mu}v^{\nu}$

with $g_{\alpha\mu}:=2f_{1}f_{2}\frac{\partial f_{1}}{\partial
x_{1}^{\alpha}}\frac{\partial f_{2}}{\partial x_{2}^{\mu}}$.

By the fact that $F_{1}$ and $F_{2}$ are Finslerian metrics, the first and the
last terms in the right hand-side of (2) are positives. Hence $g(v,v)>0$ if
the equation (2.11) holds.

When the matrix given in the Lemma 2.7 is positive-definite symmetric on
$\mathring{T}(M_{1}\times M_{2})$, it defines a Finslerian metric on
$M_{1}\times M_{2}$. In this case, ${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ is
called a convolution Finslerian metric and the Finslerian convolution manifold
${}_{f_{2}}M_{1}\bigstar_{f_{1}}M_{2}$ is called a convolution Finslerian
manifold.

Proof of the Proposition 1.2. The proof is obtained by the Lemma 2.7 together
with the Proposition 2.9.

##  3\. Important results

######  Definition 3.1.

A Riemannian metric $g$ on an $n$-dimensional manifold $M$ is a family
$\\{g_{x}\\}_{x\in M}$ where $g_{x}$ is an inner product on the tangent space
$T_{x}M$, such that the elements
$g_{ij}(x):=g_{x}\left(\frac{\partial}{\partial
x^{i}},\frac{\partial}{\partial x^{j}}\right)$ are $C^{\infty}$ in local
coordinates. Since $g_{x}$ is an inner product, the matrix $(g_{ij}(x))$ is
positive-definite at every point $x\in M$. Thus every Riemannian metric $g$ is
a Finslerian metric which arises in the following manner

$\displaystyle F(x,y):=\sqrt{g_{ij}(x)y^{i}y^{j}}.$ (3.1)

In that case, $\frac{1}{2}\frac{\partial^{2}F^{2}(x,y)}{\partial y^{i}\partial
y^{j}}$ is simply $g_{ij}(x)$ which is independant of $y$.

######  Theorem 3.2.

A convolution Finslerian metric ${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ is
Riemannian if and only if $F_{1}$ as well as $F_{2}$ are Riemannian metrics.

Proof. Assume that $F=_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ is a Riemannian metric.
By the Definition 3.1 the components $g_{ij}$ associated with $F$ are
independent of $(y_{1},y_{2})$. Hence, from the Lemma 2.7, we have

$g_{1\alpha\beta}(x_{1},y_{1})=g_{1\alpha\beta}(x_{1})\text{ and
}g_{2\mu\nu}(x_{2},y_{2})=g_{2\mu\nu}(x_{2}).$

That is $F_{1}(x_{1},y_{1})$ and $F_{2}(x_{2},y_{2})$ are Riemannian metrics.

The converse is obvious.

######  Definition 3.3.

A Finslerian metric $F$ on a manifold $M$ is locally Minkowski metric if, for
every $x\in M$, there exists an neighborhood $V$ of $x$ such that the
components $g_{ij}$ of the Hessian matrix, given in (2.1), satisfy

$\displaystyle g_{ij}(x,y)=g_{ij}(y).$ (3.2)

######  Theorem 3.4.

A convolution Finslerian metric ${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ is locally
Minkowskian if and only if $F_{1}$ as well as $F_{2}$ are locally Minkowskian
metrics and the functions $f_{k},k\in\\{1,2\\}$ are both nonzero constants.

Proof. The Theorem 3.4 is obtained by a short calculation of the locally
Minkowskian metric from the Lemma 2.7 by using the expression (3.2).

######  Definition 3.5.

A Randers metric $F$ on an $n$-dimensional manifold $M$ is a Finsler metric
which has the following form

$\displaystyle F:=\alpha+\beta$ (3.3)

where, for every $(x,y)\in\mathring{T}M$,
$\alpha(x,y):=\sqrt{a_{ij}(x)y^{i}y^{j}}$ is a Riemannian metric and
$\beta(x,y):=b_{i}(x)y^{i}$ is a $1$-form on $M$ with
$||\beta||_{\alpha}:=\sqrt{a^{ij}b_{i}b_{j}}<1$.

Proof of the Theorem 1.3. Suppose that $f_{1}$ is a nonzero constant function.
Then, since $F_{1}$ and $F_{2}$ are Randers metrics, the third term of the
right-hand side of the (2.7) vanishes. Hence, the convolution metric becomes

${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}=f_{2}^{2}(\alpha_{1}^{2}+2\alpha_{1}\beta_{1}+\beta_{1}^{2})+f_{1}^{2}(\alpha_{2}^{2}+2\alpha_{2}\beta_{2}+\beta_{2}^{2}).$
(3.4)

Since $f_{1}$ and $f_{2}$ are positive functions, denote
$f_{2}\alpha_{1},f_{1}\alpha_{2},f_{2}\beta_{1}$ and $f_{1}\beta_{2}$ by
$\alpha_{1}^{*},\alpha_{2}^{*},\beta_{1}^{*}$ and $\beta_{2}^{*}$
respectively. We obtain

${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}=\alpha_{1}^{*2}+\alpha_{2}^{*2}+2(\alpha_{1}^{*}\beta_{1}^{*}+\alpha_{2}^{*}\beta_{2}^{*})+\beta_{1}^{*2}+\beta_{2}^{*2}.$
(3.5)

By setting $\alpha=\sqrt{\alpha_{1}^{*2}+\alpha_{2}^{*2}}$ and
$\beta=\sqrt{\beta_{1}^{*2}+\beta_{2}^{*2}}$, by using the Definition 3.5 and
the relation (3.5), we get
$\alpha\beta=\alpha_{1}^{*}\beta_{1}^{*}+\alpha_{2}^{*}\beta_{2}^{*}$. Hence
the relation (1.1) follows from the fact that
$\alpha_{1}^{*}=f_{2}\alpha_{1},\beta_{1}^{*}=f_{2}\beta_{1},\alpha_{2}^{*}=f_{1}\alpha_{2}$
and $\beta_{2}^{*}=f_{1}\beta_{2}$.

######  Corollary 3.6.

A convolution Finslerian metric ${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$ is
Euclidean if and only if $F_{1}$ as well as $F_{2}$ are Euclidean metrics and
the functions $f_{k},k\in\\{1,2\\}$ are both nonzero constants.

##  4\. Some examples of Finslerian convolution metrics

######  Example 4.1.

(Finslerian convolution metric of type Riemannian) Denote by $\|.\|$ and
$\langle.,.\rangle$ the standard Euclidean norm and the inner product in
$\mathbb{R}^{n}$ respectively. Let
$\mathbb{B}_{1}=\\{x_{1}\in\mathbb{R}^{3}:\|x_{1}\|<1\\}$ and
$\mathbb{B}_{2}=\\{x_{2}\in\mathbb{R}^{2}:\|x_{2}\|<1\\}$ be two Euclidean
balls. Define the functions

$\displaystyle F_{k}:T\mathbb{B}_{k}$ $\displaystyle\longrightarrow$
$\displaystyle(0,\infty)$ $\displaystyle(x_{k},y_{k})$ $\displaystyle\to$
$\displaystyle\frac{\sqrt{\|y_{k}\|^{2}-(\|x_{k}\|^{2}\|y_{k}\|^{2}-\langle
x_{k},y_{k}\rangle^{2})}}{1-\|x_{k}\|^{2}}$ (4.1)

for $k\in\\{1,2\\}$. Then the function $F_{1}$ and $F_{2}$ are Riemannian
metrics on $\mathbb{B}_{1}$ and $\mathbb{B}_{2}$ respectively known as Klein
metrics.

Consider two smooth functions
$\rho_{k}:\mathbb{B}_{k}\longrightarrow\mathbb{R}$. Then the function

$\displaystyle F_{f_{1},f_{2}}(x_{1},x_{2},y_{1},y_{2})$ $\displaystyle=$
$\displaystyle\left\\{\frac{e^{2\rho_{2}}}{\left(1-\sum_{r=1}^{3}\big{(}x_{1}^{r}\big{)}^{2}\right)^{2}}\left[\sum_{r=1}^{3}\big{(}y_{1}^{r}\big{)}^{2}\left(1-\sum_{r=1}^{3}\big{(}x_{1}^{r}\big{)}^{2}\right)+\left(\sum_{r=1}^{3}x_{1}^{r}y_{1}^{r}\right)^{2}\right]\right.$
$\displaystyle+\frac{e^{2\rho_{1}}}{\left(1-\sum_{s=4}^{5}\big{(}x_{2}^{s}\big{)}^{2}\right)^{2}}\left[\sum_{s=4}^{5}\big{(}y_{2}^{s}\big{)}^{2}\left(1-\sum_{s=4}^{5}\big{(}x_{2}^{s}\big{)}^{2}\right)+\left(\sum_{s=4}^{5}x_{2}^{s}y_{2}^{s}\right)^{2}\right]$
$\displaystyle+\left.2e^{2\rho_{1}}e^{2\rho_{2}}\sum_{r=1}^{3}\sum_{s=4}^{5}\frac{\partial\rho_{1}}{\partial
x_{1}^{r}}\frac{\partial\rho_{2}}{\partial
x_{2}^{s}}y_{1}^{r}y_{2}^{s}\right\\}^{\frac{1}{2}}.$

is the Riemannian convolution of $F_{1}$ and $F_{2}$, via
$f_{1}:=e^{\rho_{1}}$ and $f_{2}:=e^{\rho_{2}}$.

When one the $\rho_{k}$ is a constant function or when $\nabla\rho_{k}\perp v$
for each tangent vector $v$ on $\mathbb{B}_{k}$, the Riemannian convolution
function (4.1) becomes a convolution Riemannian metric on
$\mathbb{B}_{1}\times\mathbb{B}_{2}$.

######  Example 4.2.

(Finslerian convolution metric of type Minkowskian) Let consider the
Finslerian function $F$ defined by the relation (1.2) in the Example 1.4.

It is easy to show that $F$ has the form ${}_{f_{2}}F_{1}\star_{f_{1}}F_{2}$,
with $f_{k}:U\subset\mathbb{R}^{2}\longrightarrow(0,\infty)$ such that
$f_{k}(x_{1},x_{2})=c_{k}$. That is $F$ is a Minkowskian convolution function
on $M_{1}\times M_{2}$ where $M_{1}=\\{(y^{1},y^{2})\in\mathbb{R}^{2}\\}$,
$M_{2}=\\{(y^{3},y^{4})\in\mathbb{R}^{2}\\}$, $F_{1}=\left(y_{1}^{4}+\lambda
y_{1}^{2}y_{2}^{2}+y_{2}^{4}\right)^{\frac{1}{4}}$ and
$F_{2}=\left[y_{3}^{2}+y_{4}^{2}+\lambda\left(y_{3}^{2k}+y_{4}^{2k}\right)^{\frac{1}{k}}\right]^{\frac{1}{2}}$
via $c_{1}$ and $c_{2}$.

######  Example 4.3.

(Finslerian convolution metric of type Randers) Denote by $\|.\|_{V}$ the norm
on a $3$-vector space $V$, by $\|.\|$ the standard Euclidean norm and by
$\langle.,.\rangle$ the inner product in $\mathbb{R}^{n-3}$. Let
$\mathbb{B}=\\{x_{2}\in\mathbb{R}^{n-3}:\|x_{2}\|<1\\}$ be an Euclidean ball.
Then, for any $\epsilon\in[0,1)$, the following function
$F:T\big{(}V\times\mathbb{B}\big{)}\longrightarrow(0,\infty)$ defined by

$\displaystyle F(x_{1},x_{2},y_{1},y_{2})$ $\displaystyle=$ $\displaystyle
F(x_{1}^{1},x_{1}^{2},x_{1}^{3},x_{2}^{4},...,x_{2}^{n},y_{1}^{1},y_{1}^{2},y_{1}^{3},y_{2}^{4},...,y_{2}^{n})$
(4.3) $\displaystyle=$
$\displaystyle\left\\{\|x_{2}\|^{2}\left(\|y_{1}\|_{V}^{2}+\epsilon
y_{1}^{3}\right)^{2}+2\sum_{\alpha=1}^{3}\sum_{\nu=4}^{n}x_{1}^{\alpha}x_{2}^{\nu}y_{1}^{\alpha}y_{2}^{\nu}\right.$
$\displaystyle+\left.\|x_{1}\|_{V}^{2}\left[\frac{\sqrt{\|y_{2}\|^{2}-(\|x_{2}\|^{2}\|y_{2}\|^{2}-\langle
x_{2},y_{2}\rangle^{2})}+\langle
x_{2},y_{2}\rangle}{1-\|x_{2}\|^{2}}\right]^{2}\right\\}^{\frac{1}{2}}$

for any $x_{1}=(x_{1}^{1},x_{1}^{2},x_{1}^{3})\in V$,
$x_{2}=(x_{2}^{4},...,x_{2}^{n})\in\mathbb{B},y_{1}=(y_{1}^{1},y_{1}^{2},y_{1}^{3})\in
T_{x_{1}}V\cong V$ and .$y_{2}=(y_{2}^{4},...,y_{2}^{n})\in
T_{x_{2}}\mathbb{B}\cong\mathbb{R}^{n-3}$ such that

$\|x_{2}\|^{2}\left(\|y_{1}\|_{V}^{2}+\epsilon
y_{1}^{3}\right)^{2}+\|x_{1}\|_{V}^{2}\left[\frac{\sqrt{\|y_{2}\|^{2}-(\|x_{2}\|^{2}\|y_{2}\|^{2}-\langle
x_{2},y_{2}\rangle^{2})}+\langle
x_{2},y_{2}\rangle}{1-\|x_{2}\|^{2}}\right]^{2}>-2\sum_{\alpha=1}^{3}\sum_{\nu=4}^{n}x_{1}^{\alpha}x_{2}^{\nu}y_{1}^{\alpha}y_{2}^{\nu}.$

Then the function $F$ defined in (4.3) is Randers convolution metric of
$\|y_{1}\|_{V}^{2}+\epsilon y_{1}^{3}$ and
$\frac{\sqrt{\|y_{2}\|^{2}-(\|x_{2}\|^{2}\|y_{2}\|^{2}-\langle
x_{2},y_{2}\rangle^{2})}+\langle x_{2},y_{2}\rangle}{1-\|x_{2}\|^{2}}$, via
$\|x_{1}\|$ and $\|x_{2}\|$.

## Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this
work.

## References

  * [1] D. Bao, S.-S. Chern and Z. Shen, An Introduction to Riemann-Finsler Geometry, Springer-Verlang New York, 200, (2000), 1-417.
  * [2] B.-Y. Chen, Convolution of Riemannian manifolds and its applications, Bull. Austral. Math. Soc. 66 (2002) 177-191.
  * [3] B.-Y. Chen, Differential Geometry of Warped Product Manifolds and Submanifolds, World Scientific ISBN 978-981-3208-92-6, (2016) 1-78.
  * [4] W. Kühnel and H.-B. Rademacher, Conformal vector fields on pseudo-Riemannian spaces, Diff. Geom. Appl. 7 (1997) 237-250.
  * [5] G. Nibaruta, M. Karimumuryango, A. Nibirantiza and D. Ndayirukiye, Twisted products Berwaldian metrics of polar type, Differ. Geom.-Dyn. Syst. 22 (2020), to apper.
  * [6] Y.-B Shen and Z. Shen, Introduction to Modern Finsler Geometry, Higher Education Press Limited Company and World Scientific Publishing Co. Pte. Ltd (2016) 3-38.


# Terahertz magnon frequency comb

Xianglong Yao  Zhejunyu Jin  Zhenyu Wang  Zhaozhuo Zeng  Peng Yan
<EMAIL_ADDRESS>School of Physics and State Key Laboratory of Electronic Thin
Films and Integrated Devices, University of Electronic Science and Technology
of China, Chengdu 610054, China

###### Abstract

Magnon frequency comb (MFC), the spin-wave spectra composing of equidistant
coherent peaks, is attracting much attention in magnonics. A terahertz (THz)
MFC, combining the advantages of the THz and MFC technologies, is highly
desired because it would significantly advance the MFC applications in
ultrafast magnonic metrology, sensing, and communications. Here, we show that
the THz MFC can be generated by nonlinear interactions between spin waves and
skyrmions in antiferromagnets [Z. Jin _et al_., arXiv:2301.03211]. It is found
that the strength of the three-wave mixing between propagating magnons and
breathing skyrmions follows a linear dependence on the driving frequency and
the MFC signal can be observed over a broad driving frequency range. Our
results extend the working frequency of MFC to the THz regime, which would
have potential applications in ultrafast spintronic devices and promote the
development of nonlinear magnonics in antiferromagnets.

##  I Introduction

The far-infrared electromagnetic spectrum in the frequency region of $0.1-30$
terahertz (THz), known as the THz band, has long been considered the last
remaining scientific gap in the electromagnetic spectrum Masayoshi2007 ;
Sirtori2002 . THz technology has been widely applied to many fields, such as
the wireless communication Rappaport2019 , medical imaging Pickwell2006 ,
security inspection Davies2008 , etc. Photonic and/or electronic THz devices,
like quantum cascade laser Li2014 , uni-travelling-carrier photodiode Ito2005
, Schottky-diode-based multipliers Maestrini2010 , and transistor-based
integrated circuits John2020 , have already been realized. Very recently, the
THz technology starts making its way into the field of spintronics. For
example, a novel THz emitter utilizing the spin degree of freedom in magnetic
materials has emerged Feng2021 ; Seifert2022 , which presents several
unprecedented advantages, such as the ultra-broad bandwidth (up to 30 THz)
Seifert2016  and flexible tunability by external magnetic fields Yang2016  or
internal magnetic textures Wu2022 .

Merging THz technology with other powerful techniques can yield unique
multidimensional insight into fundamental processes at ultrafast time scales.
An optical frequency comb is a spectrum consisting of a sequence of discrete
and equally-spaced spectral lines that represent precise marks in frequency,
like an optical ruler Udem2002 . The optical frequency comb technique
revolutionized optical frequency metrology and spectroscopy Kippenberg2011 ;
Fortier2019 ; Picque2019  and enabled optical atomic clocks Ludlow2015 ;
Hall2006 ; Hansch2006 . A THz optical frequency comb capable of high-
resolution measurement was recently demonstrated Burghoff2014 ; Shin2023 ,
which would significantly advance THz technology applications in spectroscopy,
metrology, sensing, and high-speed wireless communications. Most recently, it
has been reported that a magnon frequency comb (MFC) can be generated by the
nonlinear scattering between magnons and topological solitons, like skyrmion
Wang2021 , vortex Wang2022  and domain wall Beining2018 ; Zhou2021 . However,
the working frequency of MFC in ferromagnets often lies in the GHz regime,
which would not be able to keep up the demand of higher-frequency operations.
Extending the frequency range of the MFC to THz is therefore of fundamental
interest and necessary as well for ultrafast magnonics.

Antiferromagnets with two opposite magnetic sublattices have unique
advantages, such as the full freedom in magnon polarization, vanishingly small
stray field, and ultrafast magnetization dynamics typically in the THz region
Jungwirth2016 ; Baltz2018 ; Gomonay2018 . These properties make
antiferromagnets a promising platform for magnonics Rezende2019 . Inspired by
the three-wave mixing mechanism in ferromagnets, we propose to generate MFC
with a THz central frequency by nonlinear interaction between breathing
skyrmions and propagating magnons in antiferromagnets. By performing
systematical micromagnetic simulations, we find the following differences
compared with its ferromagnetic counterpart: (i) It is difficult to induce the
skyrmion breathing in antiferromagnets merely by propagating magnons and an
additional driving field is demanded. (ii) There is no frequency window to
observe the MFC. Modeling calculations reveal an unexpectedly large critical
microwave field and show a linear dependence of the three-wave coupling
strength on the microwave frequency, which explain the observed two features.

The paper is organized as follows. In Sec. II, we present the analytical model
describing the nonlinear interaction between magnons and breathing skyrmions.
The linear dispersion relation of magnons propagating in antiferromagnets is
given. The dynamical equations of magnon modes involved in three-magnon
processes are derived and numerically solved. Section III gives the full
micromagnetic simulations to verify our theoretical analysis. Conclusions are
drawn in Sec. IV.

##  II Theoretical model

We consiser an antiferromagnetic film hosting a N$\acute{\text{e}}$el-type
skyrmion stabilized by the interfacial Dzyaloshinskii-Moriya interaction
(DMI), as shown in Fig. 1(a). The Hamiltonian of the antiferromagnetic system
can be written as

$\mathcal{H}=J\sum\limits_{<i,j>}\mathbf{S}_{i}\cdot\mathbf{S}_{j}-K\sum\limits_{i}(\mathbf{S}_{i}\cdot\hat{z})^{2}+\sum\limits_{<i,j>}\mathbf{D}_{ij}\cdot(\mathbf{S}_{i}\times\mathbf{S}_{j}),$
(1)

where $\mathbf{S}_{i}$ is the unit vector of the spin at site $i$. On the
right-hand side of Eq. (1), the first term is the antiferromagnetic exchange
interaction with coefficient $J>0$. The second term is the uniaxial
magnetocrystalline anisotropy with the easy-axis along $\hat{z}$ direction and
$K$ being the anisotropic constant. The third term describes the DMI where the
DM vector $\mathbf{D}_{ij}=D(\hat{z}\times\mathbf{r}_{ij})$ with $D$ the DMI
strength and $\mathbf{r}_{ij}$ representing the position vector connecting two
neighboring spins. Under the continuum description, the antiferromagnetic
system is characterized by the net magnetization
$\mathbf{m}=(\mathbf{S}_{A}+\mathbf{S}_{B})/2$ and the staggered magnetization
$\mathbf{n}=(\mathbf{S}_{A}-\mathbf{S}_{B})/2$, where $\mathbf{S}_{A}$ and
$\mathbf{S}_{B}$ are the unit magnetization of two sublattices. Then, the
system Hamiltonian in Eq. (1) can be recast as

$\displaystyle\mathcal{H}=$ $\displaystyle\int\bigg{\\{}\frac{\lambda}{2}{\bf
m}^{2}+\frac{A}{2}\Big{[}(\nabla{\bf n})^{2}+\partial_{x}{\bf
n}\cdot\partial_{y}{\bf n}\Big{]}+L{\bf m}\cdot(\partial_{x}{\bf n}$ (2)
$\displaystyle+\partial_{y}{\bf
n})-\frac{K_{c}}{2}n_{z}^{2}+\frac{D_{c}}{2}\Big{[}n_{z}\nabla\cdot{\bf
n}-({\bf n}\cdot\nabla)n_{z}\Big{]}\bigg{\\}}d{\bf r},$

where $\lambda=8J$, $A=2Jd^{2}$, $L=2\sqrt{2}Jd$, $K_{c}=2K$, and $D_{c}=2Dd$
are the homogeneous exchange, inhomogeneous exchange, parity-breaking,
magnetic anisotropy, and DMI constants, respectively, and $d$ is the lattice
constant.

It has been demonstrated that the MFC can be generated by nonlinear coupling
between magnons and breathing skyrmion in ferromagnets Wang2021 . It is
naturally expected that this mechanism could also work in antiferromagnets. To
investigate the nonlinear interaction between antiferromagnetic magnons and
skyrmion, we express the dynamical staggered magnetization in terms of the
magnon creation and annihilation operators $a^{{\dagger}}(\textbf{r},t)$ and
$a(\textbf{r},t)$ by using the Holstein-Primakoff transformation. Then the
Hamiltonian can be rewritten as
$\mathcal{H}=\mathcal{H}^{(0)}+\mathcal{H}^{(2)}+\mathcal{H}^{(3)}+\mathcal{H}^{(4)}+\cdots$,
where $\mathcal{H}^{(0)},\mathcal{H}^{(2)},\mathcal{H}^{(3)}$, and
$\mathcal{H}^{(4)}$ are the ground state energy, two-, three-, and four-magnon
processes, respectively. The dispersion relation of antiferromagnetic magnons
is described by the second-order Hamiltonian $\mathcal{H}^{(2)}$, which can be
expressed as

$\omega=2\gamma\sqrt{(2J+K)^{2}-J^{2}[\cos(k_{x}d)+\cos(k_{y}d)]^{2}},$ (3)

with the gyromagnetic ratio $\gamma=1.76\times 10^{11}\text{rad
s}^{-1}\text{T}^{-1}$, and $\mathbf{k}=(k_{x},k_{y})$ being the wave vector of
spin waves, in the lattice model. In the continue limit, the dispersion
relation reads

$\omega=\gamma\sqrt{\lambda\bigg{(}\frac{A}{2}k^{2}+K_{c}\bigg{)}},$ (4)

with $k=|\mathbf{k}|$. Micromagnetic simulations agree well with the
analytical formula Eq. (3) [see the dashed black line in Fig. 1(b)], which
justify the validity of the antiferromagnetic spin-wave dispersion relation.
Simulation details can be found in Sec. III. For Eq. (4), the numerical
dispersion shows a good agreement in a low-$k$ region, but deviates in the
high-$k$ region [see the dashed red line in Fig. 1(b)]. It is because the
spin-wave dispersion relation Eq. (4) is valid only when the wavelength of
spin waves is much larger the lattice distance ($2\pi/k\gg d$). According to
Eq. (3) [or Eq. (4)], the antiferromagnetic resonance frequency is given by
$\omega_{\text{AFMR}}/2\pi=(\gamma/\pi)\sqrt{4JK+K^{2}}=0.369$ THz, below
which spin waves cannot propagate [see the gay region in Fig. 1(b)].

Figure 1: (a) Schematic diagram of an antiferromagnetic film hosting a
skyrmion. (b) The dispersion relation of antiferromagnetic spin waves obtained
by the FFT of the dynamical magnetization from micromagnetic simulations. The
dashed black and red lines are the analytical formulas Eq. (3) and Eq. (4),
respectively. (c) The FFT spectra of the antiferromagnetic film with (violet
curve) and without (orange curve) a skyrmion (SK). The gray regions in (b) and
(c) correspond to the band gap of spin waves. (d) The time evolution of the
skyrmion center position ($x_{s}$, $y_{s}$) and radius ($R_{s}$) under the
driving field with frequency 0.095 THz.

The nonlinear interaction between magnons and breathing skyrmion is described
by the third-order Hamiltonian $\mathcal{H}^{(3)}$, which involves four modes:
the incident magnon mode $a_{k}$, skyrmion breathing mode $a_{r}$, confluence
mode $a_{p}$, and splitting mode $a_{q}$. The dynamical equations of four
magnon modes can be derived as Slobodianiuk2019

$\displaystyle
i\frac{da_{k}}{dt}=(\omega_{k}-i\Gamma_{k})a_{k}+g_{q}a_{r}a_{q}+g_{p}a_{r}^{{\dagger}}a_{p}+h_{1}e^{-i\omega_{1}t},$
(5) $\displaystyle
i\frac{da_{r}}{dt}=(\omega_{r}-i\Gamma_{r})a_{r}+g_{q}a_{k}a_{q}^{{\dagger}}+g_{p}a_{k}^{{\dagger}}a_{p}+h_{2}e^{-i\omega_{2}t},$
$\displaystyle
i\frac{da_{p}}{dt}=(\omega_{p}-i\Gamma_{p})a_{p}+g_{p}a_{k}a_{r},$
$\displaystyle
i\frac{da_{q}}{dt}=(\omega_{q}-i\Gamma_{q})a_{q}+g_{q}a_{k}a_{r}^{{\dagger}},$

where $\Gamma_{v}=\alpha_{v}\omega_{v}\ (v=k,r,p,q)$ are the damping rate of
the magnon modes with the effective Gilbert damping constants $\alpha_{v}$ and
mode frequencies $\omega_{v}$. $g_{p}$ and $g_{q}$ are the coupling strength
of the three-magnon confluence and splitting, respectively. The incident spin
wave is excited by the microwave field $h_{1}e^{-i\omega_{1}t}$. Since the
skyrmion breathing mode in antiferromagnets is difficult to be excited merely
by the incident spin waves in a non-resonant manner, as shall be discussed
later, we apply another microwave field $h_{2}e^{-i\omega_{2}t}$ to resonantly
excite the skyrmion breathing mode. Here, $h_{i}$ and $\omega_{i}$ ($i=1,2$)
are the microwave field amplitude and frequency, respectively. Then, the
incident spin-wave mode ($\omega_{k}$) and breathing mode ($\omega_{r}$) of
the skyrmion would mix with each other and generate the sum-frequency
($\omega_{p}=\omega_{k}+\omega_{r}$) and difference-frequency
($\omega_{q}=\omega_{k}-\omega_{r}$) modes. These two secondary signals
further hybridize with the skyrmion breathing mode to generate higher-order
frequency modes, eventually leading to the MFC.

##  III Micromagnetic simulations

To verify the above picture, we perform full micromagnetic simulations using
MUMAX3 Vansteenkiste2014 . We consider a G-type antiferromagnetic film with
dimensions $1000\times 1000\times 1$ nm3. The cell size of $1\times 1\times 1$
nm3 is used to discrete the film in simulations. Magnetic parameters of KMnF3
are adopted Barker2016 ; Shen2020 : $M_{s}=3.76\times 10^{5}$ A/m,
$A^{\mathrm{sim}}=AM_{s}/4=6.59$ pJ/m,
$K^{\mathrm{sim}}=K_{c}M_{s}/2=1.16\times 10^{5}$ J/m3,
$D^{\mathrm{sim}}=D_{c}M_{s}/2=1$ mJ/m2, and $\alpha=1\times 10^{-3}$.
Absorbing boundary conditions are used to avoid the spin-wave reflection by
film edges.

To characterize the intrinsic modes of antiferromagnetic skyrmion, we apply a
sinc-function field $\mathbf{h}(t)=h_{0}\mathrm{sinc}(\omega_{H}t)\hat{z}$
with amplitude $h_{0}=10$ mT and cutoff frequency $\omega_{H}/2\pi=2.5$ THz
over the antiferromagnetic film. By carrying a standard fast Fourier
transformation (FFT) for each cell and then averaging over the whole film, we
find one main peak at 0.095 THz in the band gap [see the violet curve in Fig.
1(c)]. For the film without a skyrmion, this peak is absent [see the orange
curve in Fig. 1(c)]. This indicates that the peak at 0.095 THz comes from the
intrinsic mode of the antiferromagnetic skyrmion. By analyzing the skyrmion
motion under a sinusoidal microwave field with frequency 0.095 THz, we
identify this mode as the breathing mode ($\omega_{r}$) of the
antiferromagnetic skyrmion, as shown in Fig. 1(d). Here, the position and
radius of skyrmion are obtained by a circular curve fitting of the $m_{z}=0$
contour.

Figure 2: (a) Snapshot of the interaction between propagating magnons and
skyrmion in AFM. The incident spin waves and breathing skyrmion are excited by
the microwave fields $h_{1}\sin(\omega_{1}t)\hat{x}$ (black bar) and
$h_{2}\sin(\omega_{2}t)\hat{z}$ (blue box), respectively. The red rectangle is
the detection region. Black arrows label the Hall trajectory of magnons. (b)
The spin-wave spectrum of the detection region in (a) when only the exciting
field $h_{1}=500$ mT is applied. (c) The spin-wave spectra as a function of
the driving field amplitude ($h_{1}$). The driving frequency is fixed at 1.3
THz. (d) The FFT spectrum at three representative fields $h_{1}=0.02$ mT
(upper panel), 0.08 mT (middle panel), and 10 mT (lower panel), respectively.
In (c) and (d), the skyrmion breathing mode is excited by the second microwave
field with amplitude $h_{2}=5$ mT and frequency $\omega_{2}/2\pi=0.095$ THz.

In ferromagnets, when the amplitude of the local exciting field increases
above a threshold, the skyrmion breathing would be excited via a three-wave
splitting process. Subsequent wave-mixing would result in the MFC Wang2021 .
Inspired by this idea, we first apply a linearly-polarized microwave field
$\mathbf{h}(t)=h_{1}\sin(\omega_{1}t)\hat{x}$ with $\omega_{1}/2\pi=1.3$ THz
in a narrow rectangular area [black bar in Fig. 2(a)] to excite the incident
spin waves, which then interact with the skyrmion for generating the MFC.
However, there is no sign of the MFC except the incident spin-wave mode and
its frequency-doubling signal, even when the field amplitude increases up to
500 mT, as shown in Fig. 2(b), which is already beyond the scope of
conventional experiments. It is noted that the skyrmion breathing in chiral
antiferromagnets can be conveniently excited by a modification of the DMI or
magnetocrystalline anisotropy Qiu2021 ; Komineas2022 .

To overcome this obstacle, we apply another microwave field
$\mathbf{h}_{2}=h_{2}\sin(\omega_{2}t)\hat{z}$ with $\omega_{2}/2\pi=0.095$
THz ($\omega_{2}=\omega_{r}$) in a square region [blue box in Fig. 2(a)] to
excite the skyrmion breathing mode. Then, the incident magnons would interact
with the breathing skyrmion. Black arrows represent the trajectories of
scattered magnons, which may include both the linear and nonlinear topological
magnon spin Hall effects Jin2023 . The experimental scheme for detection shall
be discussed below. To determine whether the MFC is generated or not, we
detect the magnon spectrum in a rectangular region behind the skyrmion [red
rectangle in Fig. 2(a)] by performing FFT of the dynamical magnetization.
Figure 2(c) shows the FFT spectra by continuously varying the microwave field
amplitude $h_{1}$. Around 0.5 and 1.3 THz, two sets of magnon signals can be
observed. The modes around 0.5 THz correspond to the frequency multiplication
of the skyrmion breathing mode of frequency $n\omega_{r}$ with $n$ being an
integer. Because of the band gap limited by the antiferromagnetic resonance
$\omega_{\text{AFMR}}$, only the high-order modes above 0.369 THz ($n\geq 5$)
can escape from the skyrmion and are subsequently detected in the red
rectangle region. The modes around 1.3 THz are the MFC modes generated by
nonlinear interaction between the incident magnons and skyrmion breathing
mode. Since the skyrmion breathing mode is directly excited by the microwave
field $\mathbf{h}_{2}$ rather than the incident magnon, the driving field
amplitude is not required to be larger than the aforementioned threshold value
($\gg 500$ mT). Nevertheless, when the incident spin waves are excited by a
weak driving field ($h_{1}=0.02$ mT), the MFC signals can hardly be observed
[see the upper panel in Fig. 2(d)]. This is mainly because the amplitudes of
newly generated MFC modes are very small and decay rapidly during their
propagation in the presence of Gilbert damping. By increasing the driving
field amplitude, the MFC modes with frequency spacing equal to $\omega_{r}$
clearly emerge, as shown in the middle and lower panels of Fig. 2(d).

Figure 3: (a) The amplitudes of three main peaks at $\omega_{k}=\omega_{1}$
and $\omega_{k}\pm\omega_{r}$ as a function of the driving field amplitude
($h_{1}$). The solid lines are the fitting curves with the parameter $g=55.7$
MHz based on Eq. (5). (b) Analytical curves of the mode amplitudes, which are
obtained by numerically solving Eq. (5) with parameters $\omega_{k}/2\pi=1.3$
THz, $\omega_{r}/2\pi=\omega_{2}/2\pi=0.095$ THz, $h_{2}=0$ mT,
$\alpha=0.001$, and $g=55.7$ MHz.

The coupling strength ($g_{p,q}$) between spin waves and skyrmion is crucial
for the formation of MFC, which depends on the mode overlap and is difficult
to calculate in our case due to the lacking of the analytical skyrmion
profile. Hence, we treat them as fitting parameters. For simplicity, we assume
$g_{p,q}=g$ and that all modes have the same damping rates
$\alpha_{k,r,p,q}=\alpha$. By numerically solving Eq. (5) and fitting the
amplitudes of the incident, confluence, and splitting magnon modes at steady
state (see Fig. 6 in Appendix), we obtain the coupling strength $g=55.7$ MHz
[see Fig. 3(a)]. Substituting this coupling strength back into Eq. (5), we can
calculate the threshold field for the excitation of the skyrmion breathing
mode by propagating magnons. We identify the critical value to be
$h_{1}^{\text{cr}}=813$ mT, as shown in Fig. 3(b). It is too high compared to
the case in ferromagnets. This result provides an explanation why propagating
magnons alone is not able to split the skyrmion breathing, as noticed in our
micromagnetic simulations.

Figure 4: (a) The spin-wave spectrum as a function of the driving field
frequency $\omega_{1}$. The driving field amplitude is fixed at $h_{1}=30$ mT.
(b) The FFT amplitudes of three modes at frequencies $\omega_{1}$ (upper
panel), and $\omega_{1}\pm\omega_{r}$ (middle panel) as a function of the
driving field frequency. The lower panel in (b) shows the dependence of the
coupling strength on the driving field frequency. Magenta dots are extracted
from simulation results and the green line is the linear fitting.

Next, we investigate the dependence of the MFC on the driving frequency, as
shown Fig. 4(a). We find that, over a rather broad driving frequency range,
the MFC is still visible and the mode spacing of the MFC is always equal to
the skyrmion breathing frequency. This result significantly differs from the
MFC generated in ferromagnets, where the MFC can only be observed in a
frequency window and we attributed it to the Gaussian profile of the
frequency-dependent three-wave coupling Wang2021 . By the same fitting method
used in Fig. 3(a), we extract the coupling strengths for different driving
frequencies, as shown in Fig. 4(b). One can see that the coupling strength $g$
monotonically decreases with the increase of the driving frequency
$\omega_{1}$, approximately following a linear function $g=c\omega_{1}+g_{0}$
with the dimensionless parameter $c=-1.04\times 10^{-4}$ and intercepting
coupling $g_{0}=0.191$ THz. This result suggests that the nonlinear three-wave
coupling between magnons and breathing skyrmion in ferromagnets and
antiferromagnets has very different frequency dependencies, which might
originate from distinct internal modes of the skyrmion in ferromagnets and
antiferromagnets Kravchuk2019 . It might also be because the dispersion of
magnons is parabolic in ferromagnets, but linear in antiferromagnets. Further
investigations are necessary to elucidate the underlying physical origin
behind such interesting findings, which goes beyond the scope of this paper.

Figure 5: Schematic of experimental detection of nonlinear topological magnon
spin Hall effect. THz magnons are excited by femtosecond laser (purple arrow)
and scattered by breathing skyrmion (green arrows). Meanwhile, incident (black
arrows) and nonlinear magnons (blue/red arrows) are scattered to different
channels (gray background) due to the very different Hall angles. At the end
of each channel, the accumulated magnons with different frequencies are
detected through magneto-optical spectroscopy (blue arrows). Here, the Kerr
rotation and Faraday rotation describe the rotation of the plane of
polarization for reflected and transmitted laser, respectively. It is noted
that left-handed magnons will be scattered to top right channels (not shown).

In Ref. Jin2023 , we have identified a so-called nonlinear topological magnon
spin Hall effect by analyzing the real-space scattering patterns of the MFC.
Below, we propose a multi-channel device to measure it, by minimizing the
spatial overlap between different magnon modes, as shown in Fig. 5.
Experimentally, we can excite THz magnons in antiferromagnets by femtosecond
laser Hortensius2021 . When the generated magnons interact with the breathing
skyrmion, the MFC emerges and the incident and nonlinear magnons are scattered
to different channels due to their very different Hall angles. One can
distinguish the accumulated magnons with different frequencies by magneto-
optical effects (e.g., Faraday or Kerr rotations), as plotted in Fig. 5.

##  IV Conclusion

In summary, we theoretically demonstrated that the THz MFC can be generated in
antiferromagnetic film by nonlinear interactions between magnons and
skyrmions. Although the mechanism of the MFC generation in both ferromagnets
and antiferromagnets comes from the three-wave mixing, there exist significant
differences between them. First, the skyrmion breathing in antiferromagnets
can hardly be excited merely by propagating magnons because of an unexpectedly
large threshold field amplitude of microwaves. An additional low-frequency
driving source is therefore needed to assist the MFC generation. Second, the
dependence of the coupling strength on the driving frequency in two magnetic
systems is different. The frequency dependence follows a Gaussian profile in
ferromagnets Wang2021 , while it exhibits a linear dependence in
antiferromagnets. As a consequence, the MFC in antiferromagnets is visible
over a broad driving frequency range, in contrast to a narrow window in
ferromagnets. From the application point of view, the THz MFC we predicted can
be utilized to detect magnetic textures or defects in antiferromagnets, which
is difficult to realize by conventional means because of the vanishingly small
net magnetization. Our findings bring the MFC to the THz regime, which would
advance MFC technology applications in ultrafast magnonic metrology, sensing,
and communications.

##  V Acknowledgments

###### Acknowledgements.

We thank H. Yang, L. Song, and X. Liu for helpful discussions. This work was
funded by the National Key R$\&$D Program under Contract No. 2022YFA1402802
and the National Natural Science Foundation of China (NSFC) (Grants No.
12374103 and No. 12074057). Z.W. acknowledges financial support from the NSFC
under Grant No. 12204089. X.Y., Z.J., and Z.W. contributed equally.

## APPENDIX

In general, the dynamical equations of magnon modes Eq. (5) is hard to solve
analytically. Thus, we numerically solve it by using the ode45 solver in
MATLAB. The initial amplitude of four modes are set as $a_{k}=0$,
$a_{r}=0.001$, $a_{p}=0$, and $a_{q}=0$. After a long-time evolution (20 ns),
a stationary solution can be obtained, as shown in Fig. 6.

Figure 6: The time-evolution of four magnon modes obtained from numerically
solving Eq. (5). The exciting field of the incident magnon and skyrmion
breathing modes are set as $h_{1}=40$ mT,
$\omega_{1}/2\pi$=$\omega_{k}/2\pi=1.3$ THz, $h_{2}=5$ mT, and
$\omega_{2}/2\pi$=$\omega_{r}/2\pi=0.095$ THz. The coupling strength is
$g=55.7$ MHz and the damping constant is $\alpha=0.001$.

## References

  * (1) M. Tonouchi, Cutting-edge terahertz technology, Nat. Photon. 1, 97 (2007).
  * (2) C. Sirtori, Bridge for the terahertz gap, Nature (London) 417, 132 (2002).
  * (3) T. S. Rappaport, Y. Xing, O. Kanhere, S. Ju, A. Madanayake, S. Mandal, A. Alkhateeb, and G. C. Trichopoulos, Wireless Communications and Applications Above 100 GHz: Opportunities and Challenges for 6G and Beyond, IEEE Access 7, 78729 (2019).
  * (4) E. Pickwell and V. P. Wallace, Biomedical applications of terahertz technology,  J. Phys. D: Appl. Phys. 39, R301 (2006).
  * (5) A. G. Davies, A. D. Burnett, W. Fan, E. H. Linfield, and J. E. Cunningham, Terahertz spectroscopy of explosives and drugs, Mater. Today 11, 18 (2008).
  * (6) L. Li, L. Chen, J. Zhu, J. Freeman, P. Dean, A. Valavanis, A. G. Davies, and E. H. Linfield, Terahertz Quantum Cascade Lasers with $>1$ W Output Powers, Electron. Lett. 50, 309 (2014).
  * (7) H. Ito, F. Nakajima, T. Furuta, and T. Ishibashi, Continuous THz-wave generation using antenna-integrated uni-travelling-carrier photodiodes, Semicond. Sci. Technol. 20, S191 (2005).
  * (8) A. Maestrini, B. Thomas, H. Wang, C. Jung, J. Treuttel, Y. Jin, G. Chattopadhyay, I. Mehdi, and G. Beaudin, Schottky diode-based terahertz frequency multipliers and mixers, C. R. Physique 11, 480 (2010).
  * (9) L. John, A. Tessmann, A. Leuther, P. Neininger, T. Merkle, and T. Zwick, Broadband 300-GHz Power Amplifier MMICs in InGaAs MHEMT Technology, IEEE Trans. THz Sci. Technol. 10, 309 (2020).
  * (10) Z. Feng, H. Qiu, D. Wang, C. Zhang, S. Sun, B. Jin, and W. Tan, Spintronic terahertz emitter, J. Appl. Phys. 129, 010901 (2021).
  * (11) T. S. Seifert, L. Cheng, Z. Wei, T. Kampfrath, and J. Qi, Spintronic sources of ultrashort terahertz electromagnetic pulses, Appl. Phys. Lett. 120, 180401 (2022).
  * (12) T. Seifert, S. Jaiswal, U. Martens, J. Hannegan, L. Braun, P. Maldonado, F. Freimuth, A. Kronenberg, J. Henrizi, I. Radu, E. Beaurepaire, Y. Mokrousov, P. M. Oppeneer, M. Jourdan, G. Jakob, D. Turchinovich, L. M. Hayden, M. Wolf, M. Münzenberg, M. Kläui, and T. Kampfrath, Efficient metallic spintronic emitters of ultrabroadband terahertz radiation, Nat. Photon. 10, 483 (2016).
  * (13) D. Yang, J. Liang, C. Zhou, L. Sun, R. Zheng, S. Luo, Y. Wu, and J. Qi, Powerful and Tunable THz Emitters Based on the Fe/Pt Magnetic Heterostructure, Adv. Opt. Mater. 4, 1944 (2016).
  * (14) W. Wu, S. Lendinez, M. T. Kaffash, R. D. Schaller, H. Wen, and M. B. Jungfleisch, Controlling polarization of spintronic THz emitter by remanent magnetization texture, Appl. Phys. Lett. 121, 052401 (2022).
  * (15) Th. Udem, R. Holzwarth, and T. W. Hansch, Optical frequency metrology, Nature (London) 416, 233 (2002).
  * (16) T. J. Kippenberg, R. Holzwarth, and S. A. Diddams, Microresonator-based optical frequency combs, Science 332, 555 (2011).
  * (17) T. Fortier and E. Baumann, 20 years of developments in optical frequency comb technology and applications, Commun. Phys. 2, 153 (2019).
  * (18) N. Picqué and T. W. Hänsch, Frequency comb spectroscopy, Nat. Photon. 13, 146 (2019).
  * (19) A. D. Ludlow, M. M. Boyd, J. Ye, E. Peik, and P. O. Schmidt, Optical atomic clocks, Rev. Mod. Phys. 87, 637 (2015).
  * (20) J. L. Hall, Nobel Lecture: Defining and measuring optical frequencies, Rev. Mod. Phys. 78, 1279 (2006).
  * (21) T. W. Hänsch, Nobel Lecture: Passion for precision, Rev. Mod. Phys. 78, 1297 (2006).
  * (22) D. Burghoff, T.-Y. Kao, N. Han, C. W. I. Chan, X. Cai, Y. Yang, D. J. Hayton, J.-R. Gao, J. L. Reno, and Q. Hu, Terahertz laser frequency combs, Nat. Photon. 8, 462 (2014).
  * (23) D.-C. Shin, B. S. Kim, H. Jang, Y.-J. Kim, and S.-W. Kim, Photonic comb-rooted synthesis of ultra-stable terahertz frequencies, Nat. Commun. 14, 790 (2023).
  * (24) Z. Wang, H. Y. Yuan, Y. Cao, Z.-X. Li, R. A. Duine, and P. Yan, Magnonic Frequency Comb through Nonlinear Magnon-Skyrmion Scattering, Phys. Rev. Lett. 127, 037202 (2021).
  * (25) Z. Wang, H. Y. Yuan, Y. Cao, and P. Yan, Twisted Magnon Frequency Comb and Penrose Superradiance, Phys. Rev. Lett. 129, 107203 (2022).
  * (26) B. Zhang, Z. Wang, Y. Cao, P. Yan, and X. R. Wang, Eavesdropping on spin waves inside the domain-wall nanochannel via three-magnon processes, Phys. Rev. B 97, 094421 (2018).
  * (27) Z.-W. Zhou, X.-G. Wang, Y.-Z. Nie, Q.-L. Xia, and G.-H. Guo, Spin wave frequency comb generated through interaction between propagating spin wave and oscillating domain wall, J. Magn. Magn. Mater. 534, 168046 (2021).
  * (28) T. Jungwirth, X. Marti, P. Wadley, and J. Wunderlich, Antiferromagnetic spintronics,  Nat. Nanotech. 11, 231 (2016).
  * (29) V. Baltz, A. Manchon, M. Tsoi, T. Moriyama, T. Ono, and Y. Tserkovnyak, Antiferromagnetic spintronics, Rev. Mod. Phys. 90, 015005 (2018).
  * (30) O. Gomonay, V. Baltz, A. Brataas, and Y. Tserkovnyak, Antiferromagnetic spin textures and dynamics, Nat. Phys. 14, 213 (2018).
  * (31) S. M. Rezende, A. Azevedo, and R. L. Rodríguez-Suárez, Introduction to antiferromagnetic magnons, J. Appl. Phys. 126, 151101 (2019).
  * (32) D. V. Slobodianiuk, G. A. Melkov, K. Schultheiss, H. Schultheiss, and R. V. Verba, Nonlinear ferromagnetic resonance in the presence of three-magnon scattering in magnetic nanostructures, IEEE Magn. Lett. 10, 6103405 (2019).
  * (33) A. Vansteenkiste, J. Leliaert, M. Dvornik, M. Helsen, F. Garcia-Sanchez, and B. V. Waeyenberge, The design and verification of MUMAX3,  AIP Adv. 4, 107133 (2014).
  * (34) J. Barker and O. A. Tretiakov, Static and Dynamical Properties of Antiferromagnetic Skyrmions in the Presence of Applied Current and Temperature, Phys. Rev. Lett. 116, 147203 (2016).
  * (35) L. Shen, J. Xia, X. Zhang, M. Ezawa, O. A. Tretiakov, X. Liu, G. Zhao, and Y. Zhou, Current-Induced Dynamics and Chaos of Antiferromagnetic Bimerons, Phys. Rev. Lett. 124, 037202 (2020).
  * (36) L. Qiu, L. Shen, X. Zhang, Y. Zhou, G. Zhao, W. Xia, H.-B. Luo, and J. P. Liu, Interlayer coupling effect on skyrmion dynamics in synthetic antiferromagnets, Appl. Phys. Lett. 118, 082403 (2021).
  * (37) S. Komineas and P. E. Roy, Breathing skyrmions in chiral antiferromagnets, Phys. Rev. Research 4, 033132 (2022).
  * (38) Z. Jin, X. Yao, Z. Wang, H. Y. Yuan, Z. Zeng, W. Wang, Y. Cao, and P. Yan, Nonlinear topological magnon spin Hall effect, arXiv:2301.03211 (Phys. Rev. Lett. in press).
  * (39) V. P. Kravchuk, O. Gomonay, D. D. Sheka, D. R. Rodrigues, K. Everschor-Sitte, J. Sinova, J. van den Brink, and Y. Gaididei, Spin eigenexcitations of an antiferromagnetic skyrmion, Phys. Rev. B 99, 184429 (2019).
  * (40) J. R. Hortensius, D. Afanasiev, M. Matthiesen, R. Leenders, R. Citro, A. V. Kimel, R. V. Mikhaylovskiy, B. A. Ivanov, and A. D. Caviglia, Coherent spin-wave transport in an antiferromagnet, Nat. Phys. 17, 1001 (2021).


# RWN: A Novel Neighborhood-Based Method for Statistical Disclosure Control

Noah Perry, Department of Statistics, University of California, Davis
Norman Matloff, Department of Computer Science, University of California,
Davis
Patrick Tendick, FINRA

###### Abstract

A novel variation of the data swapping approach to statistical disclosure
control is presented, aimed particularly at preservation of multivariate
relations in the original dataset. A theorem is proved in support of the
method, and extensive empirical investigation is reported.

##  1 Introduction

The field of statistical disclosure control (SDC)—maintaining privacy of
individual records in a dataset while retaining the statistical utility of the
data—has long been the subject of arcane technical analysis, conducted mainly
by statisticians. The advent of differential privacy (DP) methods in 2006 [7]
brought computer scientists into the field, and the arcane nature of the field
changed with the well-publicized adoption of DP by the United States Census
Bureau in 2017. That move by the Bureau has become somewhat controversial [3]
[19] and though we do not address that controversy, we note the salutary
effect of greatly increasing public awareness of SDC.

Different SDC tools may be appropriate for different databases in different
settings, not just in terms of numeric degree of protection afforded by a
tool, but also in terms of usability, interpretatibility and transparency, for
end users [18]. Here we develop new methodology that we believe database
administrators (DBAs) will find useful in a variety of settings.

Our proposed method, Randomization within Neighborhoods (RWN), to be presented
below, is inspired by data swapping, a classic approach to SDC. However, RWN
differs from previous data swapping approaches, in that it exploits a certain
statistical independence property, to be described in detail in the next
section.

This paper is organized as follows. We give a brief overview of SDC methods in
Section 2, followed by a discussion of considerations in SDC specific to
databases used for statistical analysis in Sections 3 and 4. The RWN method is
presented in Section 5. The underlying theory is given in Section 6. Tuning
parameter selection is covered in Section 7. Our empirical investigation is
discussed in Section 8, and computational issues are discussed in Section 9.

Note by the way that while we refer for convenience at some points to census
data as a concerete example, our methodology is meant to be general, not
specific to the census. It could be applied to medical data, employee data and
so on.

##  2 SDC Methods

Good surveys of SDC methods are in [6] [11] [4]. In order to contrast with
RWN, we give a brief overview here.

Previously the Census Bureau had used the popular data swapping method for SDC
[10] [6] [9]. A set of key variables is defined that may render certain
individual records in the data vulnerable to disclosure. Some records,
especially those deemed most at risk to disclosure, will have the values of
their key variables swapped with those in other records, say drawn from the
same geographic region.

There are variants, notably data shuffling [17]. These methods are typically
applied one variable at a time, thus creating the concern of attenuation of
multivariate relations.

The Bureau has also used cell suppression, in which any query concerning a
very small number of database records is denied. However, in 2017, the Bureau,
after performing various simulations, decided that the swapping approach was
in danger of reconstruction attacks and turned to DP.

Another standard SDC approach is data perturbation, in which random noise is
added to achieve privacy. The Census Bureau has used this method in the past
as well,

DP is also a perturbation method. Its novelty, though, is in its ability to be
able to quantify the degree of privacy, in a manner having certain
mathematical traits, such as composability. Another difference from classic
noise addition methods is that in DP, the noise is typically added to a
statistic of interest, say a mean (global DP), rather than to the microdata
itself (local DP).

Two of the present authors, NM and PT, have a data perturbation background
background [15] [13] [22] [21], and were interested in modernizing that
approach, making a proposal in 2016 [14] for a novel SDC method inspired by
data swapping. The present work, joined by author NP, develops those ideas.

##  3 Statistical Views and Goals

Leo Breiman’s famous essay [5] on predictive modeling described a “cultural”
difference between researchers in statistics and computer science, the former
viewing the world more in terms of probabilistic behavior, the latter in terms
of algorithms. Curiously, as pointed out by statistician Larry Wasserman [23],
a similar difference later arose in the SDC field. Referring to the need in
(global) DP to develop a separate DP-compliant method for each statistical
procedure used (mean, regression analysis etc.), he noted the contrasting
views:

> CS view: Receive a query for a [specific statistical procedure], return a
> private answer.
>
> Statistics view: Give me data. Then I can: draw plots, fit models, test fit,
> estimate parameters, make predictions …

An amusement park metaphor will be useful. Under the CS approach, one must buy
a separate ticket for each ride. Some rides won’t be available at all, pending
development of tickets tailored to those rides. With statistics, one purchases
a day pass, good for all rides. Our focus here is on settings in which one
wishes to have a “statistical day pass.” We protect the data in some way, say
perturbation, then let users conduct whatever types of statistical analyses
they wish in an open-ended way.

Good arguments can be made for either view. As noted, the one-query-at-a-time
nature of the DP/CS approach enables the setting of precise guarantees of
privacy, which would be difficult or impossible in the statistics approach. On
the other hand, this means the user is restricted to only the types of queries
for which a DP version has already been developed and implemented. Linear
regression may be available, say, but not quantile regression.

(We note in passing that unlike most applications of DP, the US Census
Bureau’s DP methodology does amount to a “day pass.” This is because they add
noise to the raw data, which are cell counts in a huge contingency table,
rather than to the output of a statistical analysis, such as a mean.)

Due to our goal of developing a statistical analysis “day pass,” in the sense
Wasserman described—estimating parameters, fitting models and so on, open-
endedly—we look at the data in the usual statistical manner, i.e. as a sample
from a population. This is in contrast to many SDC applications in which the
data themselves are of primary interest, say the total count of people in a
given income range for a given census block.

A typical example might be that of a medical database, in which the privacy of
individual patients is required, but with which medical researchers can still
conduct statistical analyses, making population inferences.

We suppose here that there is some population value $\theta$ for which we wish
to obtain a sample estimate $\widehat{\theta}$, performing statistical
operations such as inference (confidence intervals, hypothesis tests). These
operations will be conducted on the perturbed data obtained by applying RWN to
our original microdata. Typically $\theta$ will be vector valued, such as a
vector of regression coefficients.

##  4 Desirable Statistical “Day Pass” Characteristics

In developing a new SDC procedure, such as our proposed RWN, these goals are
key:

  * •

Ability to handle mixed continuous and discrete/categorical data.

  * •

Preservation, to the degree possible, of not only univariate but also
multivariate distributions/relations.

  * •

Limiting the increase in size of the standard errors of $\widehat{\theta}$.

  * •

Preservation, to the degree possible of statistical inference levels related
to $\widehat{\theta}$.

On the other hand, as noted, we do not take as a goal the preservation of
marginal totals as in Census data.

Let’s elaborate a bit on these goals.

###  4.1 Handling Mixed Continuous and Categorical Data

A major obstacle to data perturbation methods is their inability to handle
discrete/categorical data. Consider a variable such as Number of Children in
Family. After noise addition, a value may become negative, an unacceptable
situation. A similar difficulty arises with categorical variables, after they
are converted to dummy (one-hot) form.

Indeed, the vast majority of the Census Bureau’s TopDown algorithm [1] is
devoted to making adjustments to negative values, and satisfying certain
constraints involving marginal totals.

RWN will be seen to handle mixed continuous and discrete/categorical data in a
simple, natural manner.

###  4.2 Preservation of Multivariate Relations

Absent some compensating feature, any change to the data arising from applying
a “day pass” SDC procedure, say perturbation or swapping, will result in
distortions of the relations between variables in the data. This will also
occur with cell suppression methods. Since analysis of multivariate relations
comprise the very core of statistics, we take as a major goal at least
approximately preserving such relations.

We are of course willing to let those relations be one aspect of the
utility/privacy tradeoff that is necessary to any disclosure avoidance
technique. Let’s call this property Multivariate Relations Attenuation
Resistance (MRAR).

The goal then is to develop an SDC method that includes MRAR, with the method
providing the user a “lever” that she can use to choose her desired
utility/privacy tradeoff level.

Comparatively little work in the SDC field has focused on MRAR. It is
mentioned only briefly in [11] and [6] — and no wonder, as MRAR is a
challenging condition to meet.

Consider noise addition methods. One actually can preserve second-order moment
structure by setting the covariance matrix of the noise to that of the data
[12] [15] [21]. But higher-order moments are lost and other distortions can
occur. And there are no obvious techniques for extending this property with
noise addition in mixed continuous/categorical variable settings.

##  5 RWN: Randomization within Neighborhoods

The method works roughly as follows. For each record in the data, we define a
neighborhood using either a Euclidean distance-based radius or k-nearest
neighbors. Then, for each record $r$ we randomly choose a subset of the
variables to perturb. For each such variable, we replace its original value by
its counterpart in a randomly chosen record in the neighborhood of $r$. A key
point is that a different random neighbor record is used for each of the
variables to be perturbed in $r$.

More formally:

Let $W=(w_{ij}),i=1,...,n,j=1,...,p$ denote our original data on $n$
individuals and $p$ variables and
$W^{\prime}=(w^{\prime}_{ij}),i=1,...,n,j=1,...,p$ be the released (i.e.
perturbed) data.

Choose neighborhood radius $\epsilon>0$, or number of nearest neighbors $k$,
and modification probability $q$. Then we form our released data $W^{\prime}$
as follows:

> For $i=1,...n$:
>
>   1. 1.
>
> Consider record $i$ in the database:
>
> $r_{i}=(w_{i1},...,w_{ip})$ (1)
>
>   2. 2.
>
> Find the set $S_{i}$ of records within the neighborhood of $r_{i}$ other
> than $r_{i}$ itself. Each neighborhood is defined to be either the
> $k$-nearest neighbors of $r_{i}$ or the set of neighbors within $\epsilon$
> distance of $r_{i}$, whichever set is larger.
>
>   3. 3.
>
> For $j=1,...,p$:
>
> With probability $1-q$, leave $w_{ij}$ unmodified, but with probability $q$,
> modify it. For variables that are chosen to be modified, we replace $w_{ij}$
> with the value in variable $j$ of a random record in the neighborhood
> $S_{i}$. As noted, there will be a different such random record for each
> $j$. This results in a perturbed data point $w^{\prime}_{ij}$. For
> unmodified variables, $w^{\prime}_{ij}=w_{ij}$.
>
>   4. 4.
>
> Store the released, modified version of $r_{i}$ as
>
> $r_{i}^{\prime}=(w_{i1}^{\prime},...,w_{ip}^{\prime})$ (2)
>
>

A key point is that in Step 3, the $p$ actions here are taken independently of
each other. In other words, the process acts as if the $p$ variables in the
data are statistically independent of each other. This would at first seem to
violate our goal of MRAR, but it is all resolved in the theorem in Section 6
below.

We call this technique Randomization Within Neighborhoods (RWN).

###  5.1 Comparison to Rank Swapping

Rank swapping [16] does data swapping using ranks rather than data values, in
order to facilitate dealing with discrete variables. It’s implemented as the
function rankSwap() in the popular R package for SDC, sdcMicro [20]. It too is
neighborhood-based. RWN differs from rank swapping in several respects:

  * •

Rank swapping’s MRAR feature is focused on bivarate relations, while RWN is
fully multivariate, exploiting the theorem in Section 6.

  * •

RWN’s neighborhoods are more complex than rank swap’s swapping ranges due the
following features:

    * –

The size of the RWN neighborhoods can differ for each record.

    * –

The distances are computed based on all variables in a record instead of one
variable at a time, so the records within a given neighborhood can be
”similar” in many different ways.

    * –

RWN duplicates data values, rather than actually swapping them. This actually
helps reduce k-anonymity issues [20].

    * –

Setting $q<1$ allows for the possibility that some values may be the original
ones.

##  6 Rationale and Theoretical Basis

Since we are using neighborhoods, one might ask, “Why not just replace entire
data rows—a given row is replaced by a neighboring row—rather than do
replacement component by component, taking different components from different
neighboring rows?” That would achieve our MRAR goal, but at the expense to too
much increase in $Var(\widehat{\theta})$.

By independently perturbing each element of a row, we can limit the increase
in standard errors of $\widehat{\theta}$. However, we then must ask whether
our method has the all-important MRAR property. The following theorem shows
that it does.

> Let $f_{X}$ be a density function for the $p$-variate vector
> $X=(X_{1},...,X_{p})$. Consider the conditional distribution of $X$, given
> that $X$ is in a small neighborhood of $t$. Then the $X_{i}$ are
> approximately independent in this distribution.

For expositional convenience, the theorem and proof will be stated for the
case $p=2$.

Theorem: Consider a bivariate random vector $(X,Y)$ having a joint density,
and set $\epsilon>0$. For any $t$ in $\mathcal{R}^{2}$, let $A_{t,\epsilon}$
denote the $\epsilon$ neighborhood of $t$. Let $F$ denote the joint cdf of
$(X,Y)$. Given $(X,Y)=t$, define $G_{t,\epsilon}$ to be the cdf of $(X,Y)$,
given that that vector is in $A_{t,\epsilon}$. Finally, given $(X,Y)=t$,
define independent random variables $U$ and $V$ to be drawn randomly from the
first- and second-coordinate marginal distributions of $G_{t,\epsilon}$,
respectively. Then

$\lim_{\epsilon\rightarrow 0}P\left(U\leq a\textrm{ and }V\leq
b\right)=F(a,b)$ (3)

for all $-\infty<a,b<\infty$.

In other words, as $\epsilon$ goes to 0, the bivariate distribution of $(U,V)$
goes to that of $(X,Y)$, even though $U$ and $V$ are independent while $X$ and
$Y$ are not independent.

$\square$

Note that (3) concerns the unconditional distribution of $(U,V)$. The latter
is a random vector in $A_{(X,Y),\epsilon}$.

Proof:

First,

$\lim_{\epsilon\rightarrow 0}U=X$ (4)

and

$\lim_{\epsilon\rightarrow 0}V=Y$ (5)

Using the Bounded Convergence Theorem, we have

$\displaystyle\lim_{\epsilon\rightarrow 0}P\left(U\leq a\textrm{ and }V\leq
b\right)$ $\displaystyle=$ $\displaystyle\lim_{\epsilon\rightarrow
0}E\left[P\left(U\leq a\textrm{ and }V\leq b~{}|~{}X,Y\right)\right]$ (6)
$\displaystyle=$ $\displaystyle\lim_{\epsilon\rightarrow 0}E\left[P(U\leq
a~{}|~{}X,Y)\cdot P(V\leq b~{}|~{}X,Y)\right]$ (7) $\displaystyle=$
$\displaystyle E\left[1_{X\leq a}\cdot 1_{Y\leq b}\right]$ (8)
$\displaystyle=$ $\displaystyle E\left[1_{X\leq a\textrm{ and }Y\leq
b}\right]$ (9) $\displaystyle=$ $\displaystyle P(X\leq a\textrm{ and }Y\leq
b)$ (10) $\displaystyle=$ $\displaystyle F(a,b)$ (11)

##  7 Neighborhoods and Tuning Parameters

The neighborhoods are formed using both the Euclidean distance-based radius
$\epsilon$ and the number of nearest neighbors $k$, which must be specified by
the user. In short, $\epsilon$ provides control over the similarity of the
data points within a neighborhood to each other while the nearest neighbor
parameter $k$ controls the minimum size for the neighborhood.

For many datasets, there are typical records as well as records with more
unusual or extreme values. Typical records will have many neighbors even for
small values of $\epsilon$ while unusual records may have zero neighbors
unless $\epsilon$ is large. If $\epsilon$ alone were used to form
neighborhoods, RWN would set these unusual records with empty neighborhoods to
contain missing values in all variables. This would protect their privacy,
which is important considering that these unusual records may correspond to
more identifiable individuals, but would result in a complete loss in their
utility.

To avoid this, the user may choose to increase $\epsilon$ until these more
extreme records have non-empty neighborhoods. However, as seen in Figure 1, an
increase in $\epsilon$ can substantially increase the size of the neighborhood
for typical records as well, causing the values in a typical record to be
mixed with very dissimilar records in the perturbation process, potentially
leading to a decrease in utility of the perturbed data.

On the other hand, using $k$ alone would impose both an upper and lower bound
on the neighborhood size. For instance, for small $k$, the neighborhood size
may be suitable for unusual records but unnecessarily small for typical
records.

Thus, having two neighborhood size parameters $\epsilon$ and $k$ gives the
Data Stewardship Organization (DSO) finer control over the perturbation of the
data to balance utility with privacy for a specific dataset.

In Figure 1 we illustrate how RWN interacts with a specific dataset and choice
of tuning parameters to form neighborhoods for perturbation. Using the bodyfat
data, we calculate the minimum distance for each record and plot it alongside
the neighborhood size (i.e. the number of records within the neighborhood) for
multiple choices of $\epsilon$ while holding constant $q=1$ and $k=5$. In the
charts, the red horizontal line depicts the value of $\epsilon$ and the gray
vertical line depicts the value of $k$.

Figure 1: Body Fat Data, Plots of Minimum Distance vs. Neighborhood Size

##  8 Empirical Investigations

We will present the results of our empirical investigation shortly. But first,
we describe our criteria.

Any SDC method is a balance of statistical utility and degree of privacy.
Empirical investigation of the method must then define measures for these two
criteria.

###  8.1 Statistical Utility

As noted, statistical analysis is at its core a matter of identifying
relations between variables. The question to consider in the SDC context is
whether the relationships that exist in the original data tend to remain
intact in the released data. Here we follow [8], who note that a reasonable
measure is to assess whether the released dataset “can obtain approximately
the same substantive [relational] results while simultaneously protecting the
privacy.” This is one of the approaches we take below, using color correlation
plots such as in Figure 2. We find that we can perturb data while retaining
broad correlation structure, including to a large extent the strength of the
correlations.

Another aspect of utility is validity of standard errors for statistical
inference purposes. For SDC methods affecting only a small portion of the
data, this is less of an issue. For RWN, one can prove, say, that the standard
errors are asymptotically valid, as $q\rightarrow 0$, and find this to hold in
our empirical work below. (See [8] for a DP solution, under certain
assumptions. As usual, the problem of discrete/categorical variables remains a
challenge.)

Since we emphasize MRAR, a utility measure is needed toward that end. We
address this by investigating how correlations, regression coefficients, and
principal components are affected by perturbation.

Most SDC methodology is aimed at estimating relationships, rather than using
those relationships for prediction of new cases. Here, we investigate that
latter aspect.

###  8.2 Privacy

A number of measures of privacy have been used in the SDC literature. The
reader is referred to the references on this, but here our choices were guided
primarily by two examples. (As already mentioned, RWN by its nature is
resistant to k-anonymity problems.)

  * (a)

The bodyfat dataset (see later section for details) has several outliers. If
an intruder to the database knows, for instance, that a certain individual has
the highest Body Mass Index of all the subjects in the data, then this
individual is at risk of the intruder identifying sensitive variables (say if
this dataset contained private health records).

  * (b)

Problems may arise involving “inliers.” A hypothetical example that has been
used in the SDC literature is that of an employee database, in which an
intruder knows that there is just one female electrical engineer. The intruder
then queries the total salaries of all female electrical engineers, and
thereby illicitly learns her salary.

A number of measures could be used to gauge potential problems of these sorts.
We address outlier issues such as in Example (a) by using Mahalanobis
distance, and via Cook’s distance for linear regression analysis. We measure
inlier problems via distance to closest neighbor.

###  8.3 Experiments

The first dataset used for the empirical experiments is the bodyfat dataset
from the mfp R package [2]. The raw dataset contains body measurements of 252
adult males and two estimates of their body fat percentage calculated using
the Brozek and Siri equations. For our analysis, we use only the body fat
percentage based on the Siri equation, eliminate 11 observations that contain
values that appear to be erroneous or biologically implausible such as body
fat percentages less than four percent, and calculate Body Mass Index (BMI)
for each individual.

For the experiments in this section, we use the cleaned bodyfat dataset as
well as several perturbed versions of this dataset created using RWN with
tuning parameters $k=5$, $q=1$, and varying values of $\epsilon$. In each
figure where a horizontal blue line is present, it denotes the result for the
unperturbed data.

Correlation:

We calculated Pearson correlation coefficients pairwise for all variables in
the bodyfat dataset. The color correlation plots in Figure 2 show that for
small $\epsilon$, the correlation coefficients in the perturbed and
unperturbed data for a given variable are of the same sign and of very similar
magnitude. This provides evidence of the MRAR property of RWN.

Figure 2: Body Fat Data, Correlation Plot

Figure 3: Body Fat Data, Correlation Plot (Large eps)

Regression:

We regressed the body fat percentage on BMI and the neck, chest, abdomen, and
hip measurements. After estimating the model on each datasets, we compared the
estimated coefficients, standard errors, and Cook’s distances.

As seen in Figure 4, there is non-trivial variation in the estimated
coefficients on the intercept, BMI, and neck variables over varying
$\epsilon$. We also note that there are numerous cases where the sign of the
estimated coefficient on the perturbed data differs from the sign of the
corresponding estimated coefficient in the model estimated on the unperturbed
data. To some extent, these results are unsurprising. It is clear from the
correlation plot that many of the body measurements have a high positive
correlation with each other. Consequently, we expect some instability in the
estimated coefficients due to multicollinearity.

Figure 4: Body Fat Data, Regression Coefficients

As shown in Figure 5, with the exception of the intercept for larger
$\epsilon$ values, the standard errors in the perturbed data appear to be
relatively stable and similar to the standard errors calculated using the
unperturbed data.

Figure 5: Body Fat Data, Regression Standard Errors

The cleaned bodyfat dataset contains one individual who is substantially
larger than all the others. Consequently, this individual would be one of the
most easily identifiable in the unperturbed data. In Figure 6, this individual
corresponds to the largest Cook’s distance in the unperturbed data. However,
even with minimal perturbation, the maximum Cook’s distance becomes much
lower, suggesting that this individual is no longer as easily identifiable in
the perturbed data.

Figure 6: Body Fat Data, Cook’s Distances

Mahalanobis distance and minimum Euclidean distance:

Figures 7 and 8 show a similar result as the Cook’s distance plots; privacy is
provided to identifiable data points for even small values of $\epsilon$.

Figure 7: Body Fat Data, Mahalanobis Distances

Figure 8: Body Fat Data, Boxplots of Minimum Euclidean Distances Before and
After Perturbation

Principal component analysis:

After scaling the data, we performed principal component analysis on both the
unperturbed and perturbed datasets. For small values of $\epsilon$, we find
that the standard deviation of the first principal component and proportion of
variance corresponding to the first principal component are slightly higher
than in the unperturbed data. As $\epsilon$ increases the variation appears to
be spread over more principal components.

Figure 9: Body Fat Data, PCA Variance Plots

###  8.4 Prediction-Oriented Assessment

Statistical applications tend to fall into one of two general categories,
which we will refer to as Description and Prediction. The former has the goal
of understanding some entity, process, effect and so on, while the latter
concerns predicting new data.

Much of the SDC literature has been aimed at the Description side of things,
estimating means, totals, regression coeficients and the like. The Prediction
side has rarely been the focus, and we turn to that aspect in this section,
particularly classification settings, the form of many modern applications.

Intuitively, SDC methods should do fairly well in classification settings.
Consider 2-class problems, for instance, where we are predicting $Y=0,1$ from
a vector of covariates $X$, say with a continuous distribution. Assuming equal
misclassification costs, the set of $t$ for which $P(Y=1|X=t)=0.5$ forms the
decision boundary. What is the effect of data perturbation?

Consider a nonparametric regression method, say random forests. If a data
point in the training set is far from the boundary, perturbation will have
little or no effect on future predictions; only points near the boundary have
much impact. One may conjecture, then, that SDC methods will not compromise
prediction ability much, at least for nonparametric regression methods. This
was confirmed in our experiments.

We first looked at pef, a dataset included in the R regtools package. This is
data on programmers and engineers, from the 2000 US census. We predict a
variable occ, which codes one of six occupations, from variables such as age,
income and education. Here we took $q=0.5$, for $k=5,10,25,50$.
Misclassification rates are as follows:

method | no SDC | k = 5 | k = 10 | k = 25 | k = 50
---|---|---|---|---|---
RF | 0.628 | 0.627 | 0.645 | 0.682 | 0.659

This dataset does not lend itself to strong predictability, with an error rate
about about 63%. However, that rate increases only slightly under SDC. The
above results were based on 25 replications, with a holdout set size of 1000.

Here is the same analysis on a second dataset, the well-known Pima diabetes
study. It’s quite different from the census data, in that it is much smaller
(768 rows, vs. 20090 for pef), thus requiring more privacy protection. On the
other hand, greater predictive accuracy is possible for this data.

method | no SDC | k = 5 | k = 10 | k = 25 | k = 50
---|---|---|---|---|---
RF | 0.241 | 0.242 | 0.233 | 0.245 | 0.241

Again, performance appears not to decline due to the privacy action, and may
even help, due to salutory smooting effects.

##  9 Large Data Sets and Computational Complexity

Data sets are growing rapidly in both size and dimension. This is driven by
many factors, including

  * •

A proliferation of data sources, including apps, back-office systems like ERP,
smart devices, and sensors.

  * •

The growth of the systems themselves as businesses and other organizations
achieve global scale.

  * •

Increases in data collection and storage capacity through networks and mass
storage.

  * •

Larger computing devices and cloud scale computing.

  * •

The transition from recording data about entities (like people) to storing
transactions like purchases to storing events like clicks. This has happened
as analysts go from trying to understand people to understanding what they do
transactionally to how they behave and where they go.

  * •

The transition from structured to semi-structured to unstructured data.
Structured data, like that found in traditional database tables, is strictly
constrained in terms of dimension. Semi-structured data like JSON objects,
which is often captured from running applications, is unconstrained in
dimension. Unstructured data, like documents, photos, or videos, is of almost
unlimited dimension. Unstructured data might not seem like a candidate for the
method described in this paper, but documents, photos, and videos are often
reduced to a set of features described by categorical, binary, or numeric
variables.

As a result, datasets can now easily be in the billions of rows with hundreds
or thousands of variables or features. But with the growth of datasets comes a
growth in risk. Larger datasets put more people at risk, and increased
dimensionality increases damage per person. Also, the increase in
dimensionality makes it easier to identify someone’s record in a dataset. We
need to be able to apply the method to these huge datasets, so the algorithm
must be reasonably efficient from a computational complexity perspective. As
we will see in the next section, the basic algorithm is computationally
expensive, but with slight modification can handle large datasets.

##  10 Computational Complexity of the Basic Method

To assess computational complexity of the method, going forward we will assume
$p$ is large but fixed, $q=1$, and $n$ is increasing. We will also assume that
the distance metric used to define closeness is arbitrary. This is reasonable,
since we need wide latitude to define closeness in different ways for very
different datasets. Under these conditions, whether we are selecting nearby
points based on epsilon neighborhood or $k$-nearest neighbor, the
computational complexity of the basic method is $O(n^{2})$, since we have to
calculate the distance between all possible points. This is probably not
tenable for large datasets, e.g., $=10^{9}$, for which the number of
calculations would be of the order $10^{18}$.

##  11 Alternative Methods

###  11.1 Method 1: Draw Neighbors from a Sample of Points

For the first method, we will just take a smaller sample of data points to use
as neighbors. That is, we will take a sample S of size $m<<n$. For each data
point in the original data set, select the neighbors from the sample, then
apply the method as before. For this method, the complexity is $O(mn)$,
considerably better than the original method. For $n=10^{9}$, $m=10^{4}$,
complexity is of order $10^{13}$ instead of $10^{18}$, an improvement of five
orders of magnitude. Assuming that all the data is in memory, the cost of
generating a sample of size $m$ is $O(m)$, so if we were to use a distinct
random sample for each point, the complexity would still be $O(mn)$. For
$n=10^{9}$, $m=10^{4}$, complexity is still order $10^{13}$. The advantage is
that we potentially get slightly richer data and better protection.

###  11.2 Method 2: Sample Randomly from the Distance Matrix

In Method 2, we are actually sampling twice as many points as we need to,
since when we sample point $x_{j}$ to be a possible neighbor of the point
$x_{i}$, we can also use $x_{i}$ as a potential neighbor of $x_{j}$. This is
due to the symmetry of distance metrics and the distance matrix $D=[d_{ij}]$
where $d_{ij}$ is the distance between points $i$ and $j$. So instead of
drawing a new sample for each data point, we could simply draw a sample of
elements from the entire distance matrix. Then for each data point, we could
look at the distances in the sample and find those that are sufficiently small
or find the smallest $k$ elements.

We will now describe the method more formally and derive the complexity. Let
$W=(x_{ij})$ be the original data set. We want to sample pairs
$(x_{i},x_{j})$, $i<j$ at random with a sample size of $n_{s}=n{m}/2$ where
$m<<n$ is the desired sample size per data point. There are $n(n-1)/2$ such
pairs total, and we are sampling $n_{s}<<n(n-1)/2$. To generate the sample, we
can use a mixed congruential random number generator to generate numbers in
$[1,n(n-1)/2]$.

Once we have $n_{s}$ random integers in the range $[1,n(n-1)/2],$ we will map
them to $(i,j)$ pairs, $1<=i<j<=n$ to obtain the set $S_{d}=\\{(i,j)\\}$. We
then calculate the distances $d_{ij}$ = $<xi,xj>$ and store them in an
undirected graph $G_{d}$ with nodes representing the $n$ data points and edges
$\\{e_{ij}:(i,j)\in S_{d}\\}$ To map a random integer $r$ to an $(i,j)$ pair,
we will use the following algorithm:

  1. 1.

Calculate $v=\sqrt{(8r+1)}$ as an integer operation. That is, the calculation
of $v$ produces an integer plus a remainder that indicates whether $v$ is the
exact square root.

  2. 2.

If $v$ is odd and the exact square root,
$j\leftarrow(1+v)/2$
$i\leftarrow j-1$
and we’re done.

  3. 3.

Else

    1. (a)

If $v$ is odd but not the exact root,
$j\leftarrow(3+v)/2$

    2. (b)

Else ($v$ even)
$j\leftarrow(2+v)$
$i\leftarrow r–(j-1)(j-2)/2$

Since the complexities of generating the random numbers, mapping the random
numbers to $(i,j)$ pairs, calculating distances, and finding neighbors are all
$O(nm)$, the complexity of the entire algorithm is $O(nm)$.

###  11.3 Method 3: Partitioning

With the advent of cloud computing, it has become much more feasible to
harness the power of many computers or virtual machines, each with gigabytes
of memory. Under this scenario, the problem (and data) are partitioned and
spread across instances. There are several ways this approach could be applied
to the problem at hand:

  1. 1.

Simply partition the dataset into $u$ equal partitions and then apply the
original method. In this case, the size of each partition is $n/u$, so the
complexity is $O(n^{2}/u^{2})=O(n^{2})/u^{2}$, a reduction by a factor of
$u^{2}$ over the original method.

  2. 2.

Partition the dataset into u equal partitions and then apply the Method 2
above. In this case, the size of each partition is $n/u$, so the complexity is
$O(mn/u^{2})=O(n^{2})m/nu^{2}$, a reduction by a factor of $nu^{2}/m$ over the
original method.

Partitioning has two advantages: It both spreads the work across multiple
instances, thereby increasing the computing power that can be brought to bear
on the problem, and it also reduces the complexity of the problem being solved
on each node. However, the resulting combined dataset will be somewhat
different, since it is a combination of perturbed subsets.

## References

  * [1] Abowd, J., Ashmean, R., Garfinkel, S., Kifer, D., Leclerc, P., Machanavajihala, A., Moran, B., and Sexton, W. Census topdown: Differentially private data, incremental schemas, and consistency with public knowledge, 2019.
  * [2] Ambler, G. mfp: Multivariable Fractional Polynomials, 2015.  R package version 1.5.2.
  * [3] Bahrampour, T., and Lang, M. J. New system to protect census data may compromise accuracy, some experts say.  Washington Post (6 2021).
  * [4] Brand, R. Microdata Protection through Noise Addition.  Springer Berlin Heidelberg, Berlin, Heidelberg, 2002, pp. 97–116.
  * [5] Breiman, L. Statistical Modeling: The Two Cultures (with comments and a rejoinder by the author).  Statistical Science 16, 3 (2001), 199 – 231.
  * [6] Duncan, G., Elliot, M., and Salazar, G. Statistical Confidentiality: Principles and Practice.  Statistics for Social and Behavioral Sciences. Springer New York, 2011\.
  * [7] Dwork, C. Differential privacy.  In Proceedings of the 33rd International Conference on Automata, Languages and Programming - Volume Part II (Berlin, Heidelberg, 2006), ICALP’06, Springer-Verlag, p. 1–12.
  * [8] Evans, G., King, G., Schwenzfeier, M., and Thakurta, A. Statistically valid inferences from privacy protected data.  https://gking.harvard.edu/dp, 2021.
  * [9] Fienberg, S. E., and McIntyre, J. Data swapping: Variations on a theme by dalenius and reiss.  In Privacy in Statistical Databases (Berlin, Heidelberg, 2004), J. Domingo-Ferrer and V. Torra, Eds., Springer Berlin Heidelberg, pp. 14–29.
  * [10] Fischetti, M., and Salazar-González, J.-J. Partial cell suppression: A new methodology for statistical disclosure control.  Statistics and Computing 13, 1 (Feb. 2003), 13–21.
  * [11] Hundepool, A., Domingo-Ferrer, J., Franconi, L., Giessing, S., Nordholt, E., Spicer, K., and de Wolf, P. Statistical Disclosure Control.  Wiley Series in Survey Methodology. Wiley, 2012.
  * [12] Kim, J. J. Method for limiting disclosure in microdata based on random noise and transformation.  In Proceedings of the Section on Survey Research Methods (JSM) (1986), pp. 303–308.
  * [13] Matloff, N., and Tendick, P. The ’curse of dimensionality’ in database security.  In Database Security, II: Status and Prospects. Result of the IFIP WG 11.3 Workshop on Database Security, Kingston, Ontario, Canada, October 5-7, 1988 (1988), C. E. Landwehr, Ed., North-Holland, pp. 225–232.
  * [14] Matloff, N., and Tendick, P. A new method for avoiding data disclosure while automatically preserving multivariate relations, 2015.
  * [15] Matloff, N. S. Another look at the use of noise addition for database security.  In 1986 IEEE Symposium on Security and Privacy (1986), pp. 173–173.
  * [16] Moore, R. Controlled data-swappng techniques for masking public use microdata sets.  https://www.census.gov/srd/CDAR/rr96-04_Controlled_DataSwapping.pdf, 1996.
  * [17] Muralidhar, K., and Sarathy, R. Data shuffling: A new masking approach for numerical data.  Management Science 52, 5 (2006), 658–670.
  * [18] Oberski, D. L., and Kreuter, F. Differential privacy and social science: An urgent puzzle.  Harvard Data Science Review 2, 1 (1 2020).
  * [19] Rathindra Sarathy, K. M. Evaluating laplace noise addition to satisfy differential privacy for numeric data.  Transactions on Data Privacy 4 (4 2011).
  * [20] Templ, M., Kowarik, A., and Meindl, B. Statistical disclosure control for micro-data using the r package sdcmicro.  Journal of Statistical Software 67, 4 (2015), 1–36.
  * [21] Tendick, P. Optimal noise addition for preserving confidentiality in multivariate data.  Journal of Statistical Planning and Inference 27 (1991), 341–353.
  * [22] Tendick, P., and Matloff, N. A modified random perturbation method for database security.  ACM Trans. Database Syst. 19, 1 (Mar. 1994), 47–63.
  * [23] Wasserman, L. A statistical view of differential privacy.  www.cs.cmu.edu/afs/cs/usr/wing/www/class/15-895/LarryWasserman.pdf, 2011\.