∎

e1e-mail<EMAIL_ADDRESS>
11institutetext: School of Basic Sciences
Indian Institute of Technology Bhubaneswar,
Jatni, Khurda, Odisha 752050 India

# Warp Drives and Martel–Poisson charts

Abhishek Chowdhurye1,addr1

(Received: date / Accepted: date)

###### Abstract

We extend the construction of Alcubierre–Natário class of warp drives to an
infinite class of spacetimes with similar properties. This is achieved by
utilising the Martel–Poisson charts which closely resembles the Weak
Painlevé–Gullstrand form for various background metrics (Mink, AdS, dS). The
highlight of this construction is the non–flat intrinsic metric which in three
dimensional spacetimes introduces conical singularities at the origin and in
higher dimensions generates a non–zero Ricci scalar for the spatial
hypersurfaces away from the origin. We analyse the expansion/contraction of
space and $(-ve)$ energy densities associated with this class of warp drives
and find interesting deviations due to the global imprints of the conical
defects. Several generalizations are also discussed.

###### Keywords:

Warp Drive Conical defects Martel–Poisson ChartsPainlevé–Gullstrand

††journal: Eur. Phys. J. C

##  1 Introduction

The idea of “effective” superluminal travel is of perpetual interest focusing
mainly on “gedanken–experiments” pushing the foundations of general relativity
and communication. Like wormhole geometries 10.1119/1.15620 , warp drive
spacetime first proposed by Alcubierre in 1994 Alcubierre_1994  and then
further elaborated upon by Natário in Nat_rio_2002  allows for arbitrary large
velocities within the framework of general relativity. It appears as a
“bubble” propagating on some spacetime background such that the observers
inside the bubble are in an inertial reference frame requiring no external
energy sources to accelerate. These “reaction–less drives” have a severe
limitation, sustaining them violates all known energy conditions even for
non–relativistic velocities and would probably require exotic matter Lobo_2004
. Another interesting limitation is the lack of control experienced by an
observer inside a superluminal warp bubble which has been addressed by the
Krasnikov tube Everett_1997 .
In this paper we have generalized the construction of Alcubierre like warp
drive in three dim Minkowski spacetime by leveraging the Martel–Poisson (MP)
family of charts Martel_2001  of which the Painlevé–Gullstrand (PG)
coordinates are a special case. This construction has been extended to
$\text{AdS}_{3}$, dS backgrounds and their conical defect counterparts. A
generic key feature is the appearance of a non–flat intrinsic metric with
conical singularity which affect the global properties of the spacetime.
Another important realization comes from the flow velocities of the background
metrics where the expansion and $(-ve)$ energy density of the warp drive is
considerably reduced if the drive travels along the background flow
velocities. This will be useful in analog gravity setups which imitate the
warp metric Fischer_2003 ; finazzi2012analogue  and can be studied
experimentally PhysRevLett.105.240401 . Throughout the paper we have used
natural units $G=c=1$ and metric signature $(-,+,+,\cdots)$.
The rest of the paper is organized as follows. In section 2 we review and
extend the construction of MP–charts for spacetimes of the form (1) in $D$
dimensions. Interpretation of the flow coordinates are discussed including the
special case of PG–coordinates. Section 3 highlights the conical singularities
of spatial hypersurfaces in three dimensional Minkowski and AdS spacetimes and
their conical defect cousins. Regularization of the singularities and
embedding of the warp drive in these spacetimes are discussed. We conclude in
section 4 with comments on possible explorations to higher dimensions,
analogue gravity, quantum effects and optimizations. The paper ends with the
relevant references.

##  2 Martel–Poisson charts

For a large class of static spherically symmetric spacetimes in $D\geq 3$
dimensions the metric takes the form

$ds^{2}=-fdt^{2}+\frac{dr^{2}}{f}+r^{2}d\Omega_{(D-2)}^{2}$ (1)

where $d\Omega_{(D-2)}^{2}$ is the line element on an unit $(D-2)$ sphere and
$f(r)=(1+\tilde{f}(r))$ 111See later sections 3 and 4 for generalizations..
The main idea is to describe the spacetime in terms of radially infalling or
outgoing timelike geodesic observers with D–velocities
$u^{\mu}=\frac{dx^{\mu}}{d\tau}=(\dot{t},\dot{r},0,\cdots,0)$ where $\tau$ is
the proper time along the geodesic. For static spacetimes, the Killing vector
$\xi=(\partial/\partial t)$ leads to a conserved energy per unit mass
$\epsilon=-\xi_{\mu}u^{\mu}=f\dot{t}$ and utilizing the normalization of the
D–velocities $u_{\mu}u^{\mu}=-f\dot{t}^{2}+(\dot{r}^{2}/f)=-1$, we can finally
write

$u_{\mu}=\left(-\epsilon,\mp\frac{\sqrt{\epsilon^{2}-f}}{f},0,\cdots,0\right)\,,$
(2)

where the ($-ve$) and ($+ve$) signs are for the infalling and outgoing
geodesics respectively. So far, we have not established that these
D–velocities are indeed geodesics. We demand that this congruence of
D–velocities are hypersurface orthogonal to timelike surfaces $T(t,r)$ =
constant i.e. $u_{\mu}=-(\partial T/\partial x^{\mu})$ where $T$ is now the
proper time for these free falling (geodesic) observers. Redefining the time
coordinate as 222Note that $T$ is an exact differential i.e.
$\displaystyle\frac{\partial\epsilon}{\partial
r}=\frac{\partial\left(\frac{\sqrt{\epsilon^{2}-f}}{f}\right)}{\partial t}\,$.

$dT=\epsilon\mskip 1.5mudt\pm\frac{\sqrt{\epsilon^{2}-f}}{f}dr\,,$ (3)

the metric in (1) takes the form

$ds^{2}=-dT^{2}+\frac{1}{\epsilon^{2}}\left(dr+v(r)\mskip
1.5mudT\right)^{2}+r^{2}d\Omega_{(D-2)}^{2}\,,$ (4)

where $v(r)=\pm\sqrt{\epsilon^{2}-f}$ is the flow velocity. The coordinate
charts are valid and cover regions for which $v\in\mathbb{R}$. The geodesic
D–velocities of these Eulerian observers now take the form

$u^{\mu}=\left(\dot{T},\dot{r},0,\cdots,0\right)=(1,-v,0,\cdots,0)\,.$ (5)

Note that except for the case $\epsilon=1$, the induced metric on the
hypersurface is not Ricci flat. In fact, as we will see later the scalar
curvature and other curvature invariants of the induced metric blows up at
$r=0$.
At this point we can have two different interpretations for the parameter
$\epsilon$ which we will now illustrate using the familiar example of
Schwarzschild black holes in $D$ dimensions where $\displaystyle
f(r)=\left(1-\frac{r_{h}}{r^{(D-3)}}\right)$, $r_{h}$ being the location of
the horizon and the Eulerian observers are infalling.

#### Free fall from finite radius

From the D–velocity we have a conservation equation

$\frac{1}{2}\left(\frac{dr}{dT}\right)^{2}-\frac{1}{2}\left(\frac{r_{h}}{r^{(D-3)}}\right)=\frac{\epsilon^{2}-1}{2}=E\,,$
(6)

where $E$ is the total K.E. + P.E. of the observer. We can choose the
observers to start at a finite radius $r_{0}$ with zero velocity, then
$E=-\frac{1}{2}\frac{r_{h}}{r_{0}^{(D-3)}}$ but the coordinate chart is not
valid beyond $r>r_{0}$ as $v(r)$ becomes imaginary. However, if the observers
start at infinity then $E=0$ or $\epsilon=1$ which gives the PG–coordinates
for the Schwarzschild black holes where the induced metric on $T$ = constant
hypersurfaces are flat.

#### Free fall with initial velocity

The observers can start at infinity but now with some initial velocity
$v_{\infty}$ w.r.t. time coordinate $t$ such that
$\epsilon=\frac{1}{\sqrt{1-v_{\infty}^{2}}}\geq 1$. Again, zero initial
velocity gives the PG–coordinates for the Schwarzschild black holes.
There is another interesting limit where we can set $v_{\infty}=1$ i.e. the
speed of light. Scaling $\tilde{T}=(T/\epsilon)$ in (3) and then taking the
limit $\epsilon\rightarrow\infty$ gives the infalling Eddington–Finkelstein
(EF) coordinates Faraoni2020WhenPC .
Of course, it is also possible to have finite radius non–zero initial velocity
observers but for the remainder of this paper we will focus on the initial
velocity interpretation of the parameter $\epsilon$.

##  3 Warp Drive in various spacetimes

While the original Alcubierre warp drive is a deformation of the Minkowski
spacetime, we will now discuss embedding the drive in various spherically
symmetric static spacetimes, see also garattini2023black ;
garattini2024wormholewarp . Our focus will be on the minimal non–trivial
dimension $D=3$ from which it is relatively clear how to generalize for higher
dimensions. Three dimensions is also special as it admits Conical defect
solutions which change the global properties of the spacetime retaining the
local ones.

###  3.1 Minkowski spacetime

The MP–charts for the Minkowski spacetime in three dimensions with a tunable
parameter $\epsilon$ takes the form,

$ds^{2}=-dT^{2}+\frac{1}{\epsilon^{2}}\left(dr+\sqrt{\epsilon^{2}-1}\,dT\right)^{2}+r^{2}d\phi^{2}\,,$
(7)

where the radially infalling geodesic observers have the velocity
$v_{\infty}=\sqrt{1-\frac{1}{\epsilon^{2}}}$ at spatial infinity. Though this
metric has non–vanishing Christoffel symbols, it is just an artifact of the
coordinate choice and the curvature tensors vanish.

#### Conical singularity

It is instructive to analyze the induced metric on the $T$ = constant
hypersurfaces,

$ds^{2}_{C_{\epsilon}}=\frac{1}{\epsilon^{2}}dr^{2}+r^{2}d\phi^{2}\quad
0\leq\phi<2\pi\,.$ (8)

By carrying out the substitutions $r=\epsilon\tilde{r}$ and
$\phi=\epsilon\tilde{\phi}$ such that $0\leq\tilde{\phi}<2\pi\epsilon\,$ it is
straightforward to realize that this is a metric on a conical space
$C_{\epsilon}$ with a deficit angle $2\pi(1-\epsilon)$. Like the plane, the
cone is everywhere flat except at the tip where its curvature $R$ is singular.
Calculations by means of the standard formulas of the Riemannian geometry
based on existence of tangent spaces would fail to reveal this delta–like
singularity. A natural recipe to handle such singularities and their
generalizations can be found in Fursaev_1995 . To summarize, we consider an
embedding of $\mathcal{C}_{\epsilon}$ in a three dimensional Euclidean space
with the parametrization $x=\epsilon\mskip
1.5mu\tilde{r}\cos(\tilde{\phi}/\epsilon),\,y=\epsilon\mskip
1.5mu\tilde{r}\sin(\tilde{\phi}/\epsilon),\,z=\sqrt{\left|1-\epsilon^{2}\right|}\mskip
1.5mu\tilde{r}$ defining the conical surface

$z^{2}-\frac{\left|1-\epsilon^{2}\right|}{\epsilon^{2}}\left(x^{2}+y^{2}\right)=0\quad,\quad
z\geq 0\quad.$ (9)

For $\epsilon\neq 1$, there is a singularity at $z=0$ where we can’t introduce
a tangent space and calculate the curvature in the usual way.
We can however, introduce a regularization by rolling off the cone tip. The
simplest one corresponds to changing $\mathcal{C}_{\epsilon}$ to a hyperbolic
space with a tunable regularization parameter $a$,

$\displaystyle
z^{2}-\frac{\left|1-\epsilon^{2}\right|}{\epsilon^{2}}\left(x^{2}+y^{2}\right)=a^{2}\quad,\quad
z\geq 0,$ (10) $\displaystyle
ds_{H}^{2}=\frac{\tilde{r}^{2}+a^{2}\epsilon^{2}}{\tilde{r}^{2}+a^{2}}d\tilde{r}^{2}+\tilde{r}^{2}d\tilde{\phi}^{2}\,.$
(11)

To evaluate the curvature scalar, we first evaluate the integral curvature on
the regularized hyperbolic space in the $a\rightarrow 0$ limit,

$\lim_{H\rightarrow
C_{\epsilon}}\int_{H}R=4\pi(1-\epsilon)+\int_{C_{\epsilon}-\\{\tilde{r}=0\\}}R\,,$
(12)

where the second integral is zero. In a more general setting one can embed the
conical defect in a non Ricci flat two dimensional manifold and the second
integral would collect the curvature around the defect. We note that as only
the singular point $\tilde{r}=0$ can give rise to the first term in (12), it
is best to introduce a local representation for the curvature on
$C_{\epsilon}\,$,

${}^{(\epsilon)}R=\frac{2(1-\epsilon)}{\epsilon}\delta(\tilde{r})+R\,.$ (13)

#### The warp drive

We can embed the Alcubierre warp drive in the infalling Minkowski spacetime
background (7). It is crucial to note that this embedding can’t be achieved by
a coordinate transformation of the warp drive in the usual $\epsilon=1$
charts. Therefore, it represents an infinite class of different embeddings
parameterized by $\epsilon$ or equivalently $v_{\infty}$.
We begin, by writing the metric in the standard ADM form in the Cartesian
coordinates (warp drive breaks the spherical symmetry),

$ds^{2}=-Ndt^{2}+\gamma_{ij}\left(dx^{i}+N^{i}dt\right)\left(dx^{j}+N^{j}dt\right)\,,$
(14)

where $N$ is the lapse, $N^{i}$ are the shift vectors and $\gamma_{ij}$ are
the components of the induced metric on the $T$ = constant timelike
hypersurfaces. For the metric as in (7), the embedding is as follows,

$\displaystyle N$ $\displaystyle=1$ (15) $\displaystyle N^{i}$
$\displaystyle=(1-f(r_{s}))\sqrt{(\epsilon^{2}-1)}\,\frac{x^{i}}{r}-\delta^{i}_{x}\,f(r_{s})\,v$
(16) $\displaystyle\gamma_{ij}$
$\displaystyle=\delta_{ij}+\left(\frac{1}{\epsilon^{2}}-1\right)\frac{x^{i}x^{j}}{r^{2}}$
(17) $\displaystyle\gamma^{ij}$
$\displaystyle=\delta^{ij}+\left(\epsilon^{2}-1\right)\frac{x^{i}x^{j}}{r^{2}}\,,$
(18)

where for three dimensions, $r=\sqrt{x^{2}+y^{2}}$ is the radial coordinate,
$r_{s}$ is the distance from the center of the drive which is moving along the
($+ve$) $x$–direction with constant velocity $v$. Here, the form function
$f(r_{s})$ decides the size and shape of the wrap drive bubble, its exact form
is of little consequence (unless one is optimizing) as long as it has the
value $f(0)=1$ at the center of the bubble, $f\rightarrow 0$ as
$r\rightarrow\infty$ outside the bubble and drops sharply from $1$ to $0$ at
the wall of the bubble at some distance $R$ from the center of the bubble
333One such specific function popular in the literature is
$f(r_{s})=\frac{\tanh[\sigma(r_{s}+R)]-\tanh[\sigma(r_{s}-R)]}{2\tanh(\sigma
R)}\,,$ (19) where $\sigma$ has the interpretation of being inversely
proportional to the bubble wall thickness and can be taken to be large.. The
choice of shift vectors $N^{i}$ are inspired by ellis2004warp , where the flow
velocity is that of the warp drive inside the bubble and that of the
background spacetime outside the bubble. The Eulerian observers with geodesics
$u_{\mu}=(-1,0,0)$ are still in free fall and hence the above metric can be
considered as a warp drive different from the Natário type Nat_rio_2002 . We
assume the warp drive to be a test particle moving along the timelike curve
$x=x_{0}(T)=vT$, regardless of the value of $v(T)$. One can verify that the
proper time along this curve equals the coordinate time. Note that here, the
induced metric is not flat and as discussed earlier will exhibit conical
singularity at the origin. Warp drives with non flat $\gamma_{ij}$ have been
discussed earlier in Broeck_1999 ; garattini2024wormholewarp . As expected,
for $\epsilon=1$, we get back the original Alcubierre warp drive
Alcubierre_1994 .
Few comments on the distance $r_{s}$ and the induced metric $\gamma_{ij}$ are
in order. We can always choose $r_{s}$ to be the coordinate distance
$r_{s}=\sqrt{(x-x_{0}(T))^{2}+y^{2}}$ from the center of the warp drive. In
general, the induced metric $\gamma_{ij}$ is not $\delta_{ij}$, therefore the
size & shape of the bubble inside the boundary $r_{s}=R$ will evolve over
time. Even without the warp drive, the metric in (7) has non–zero
expansion/contraction and shear for a spacelike patch orthogonal to the
congruence of Eulerian geodesic observers 444Non zero rotation is not possible
as Eulerian geodesics are hypersurface orthogonal.. Our focus will be on the
large $x_{0}$ and $(\epsilon-1)$ small limit where these background expansion
and shear are at check. However, there are other possible choices, either
$r_{s}=\sqrt{\gamma_{ij}(x-x_{0}(T))^{i}(x-x_{0}(T))^{j}}$ or
$r_{s}=\text{Geodesic}(x,x_{0}(T))$ where $\text{Geodesic}(x,x_{0}(T))$ is the
geodesic distance between the center of the warp drive $(x_{0}(T),0)$ and some
point $(x,y)$. Another possible choice for the induced metric is
$\gamma_{ij}=\delta_{ij}+(1-f(r_{s}))\left(\frac{1}{\epsilon^{2}}-1\right)\frac{x^{i}x^{j}}{r^{2}}\,$,
an interpolating metric which is $\delta_{ij}$ inside the bubble and
$\gamma_{ij}$ as in (17) outside 555It is difficult to implement the other
choices of $r_{s}$ with this choice of induced metric as now the definition of
$r_{s}$ depends on $r_{s}$ itself through the function $f(r_{s})$. This puts
severe constraints on possible $f(r_{s})$, which we checked to be incompatible
with the desired properties of $f(r_{s})$.. We will discuss some of these
cases in A.
With our first choice of $r_{s}$ and the induced metric (17) we shall now
discuss the movement of the warp bubble due to the expansion and contraction
of the space around it. The space expands behind the warp bubble and contracts
in front of it, thereby pushing the bubble forward (and maybe a spaceship
inside it) at arbitrary velocity $v$ 666The same mechanism also applies to the
the expanding universe, wherein galaxies move faster than light w.r.t. each
other due to the expansion of the space itself.. We shall also discuss the
energy density as measured by the Eulerian observers i.e.
$\rho=T_{\mu\nu}u^{\mu}u^{\nu}$ where $T^{\mu\nu}$ is the stress energy tensor
of the metric and $u^{\mu}=(1,-N^{x},-N^{y})$ are the 3–velocities.
To evaluate the expansion $\theta$ and energy density $\rho$ associated with
the warp drive spacetime, we shall compute at first the extrinsic curvature
$K_{\mu\nu}=\nabla_{\mu}u_{\nu}=\frac{1}{2}\mathcal{L}_{u}P_{\mu\nu}$ for the
$T$ = constant hypersurfaces embedded in the full spacetime. Here, the lapse
is $N=1$, the timelike geodesics $u^{\mu}$ are hypersurface orthogonal and
$P_{\mu\nu}=g_{\mu\nu}+u_{\mu}u_{\nu}$ is the projector to the hypersurfaces.
The pullback of $K_{\mu\nu}$ to the above hypersurfaces are

$\displaystyle K_{ij}=\frac{1}{2}\mathcal{L}_{u}\gamma_{ij}$
$\displaystyle=\frac{1}{2}\left(\partial_{T}-\mathcal{L}_{N}\right)\gamma_{ij}$
(20)
$\displaystyle=\frac{1}{2}\partial_{T}\gamma_{ij}-\frac{1}{2}\left(\gamma_{kj}D_{i}N^{k}+\gamma_{ik}D_{j}N^{k}\right)\,,$

where $D_{i}$ are the covariant derivatives w.r.t. the induced metric
$\gamma_{ij}$. The expansion is the trace $\theta=K^{i}_{i}$. The Einstein
tensor $G_{\mu\nu}=R_{\mu\nu}-\frac{1}{2}R\,g_{\mu\nu}$ along with the
Gauss–Codazzi relations which connects curvature tensors of the metric to
curvatures on the hypersurface(s) give,

$\displaystyle{}^{(2)}R+\left(K^{2}-K^{ij}K_{ij}\right)=2\,G_{\mu\nu}u^{\mu}u^{\nu}=16\pi\rho+2\Lambda$
(21) $\displaystyle D_{j}K^{j}_{i}-D_{i}K=8\pi p_{i}\,,$

where $\rho$ is the matter energy density and $p_{i}$ are the matter momentum
density and $\Lambda$ is the cosmological constant.
For the background metric in (7) without the warp drive, the extrinsic
curvature and expansion away from the conical singularity at $r=0$ are given
by 777The regularization mentioned in (11) needs modification as now we want
${}^{(2)}R+\left(K^{2}-K^{ij}K_{ij}\right)$ to vanish at $r=0$ as $\rho=0$.
The regularized metric around $r=0$ with $0\leq\phi<2\pi$ is given by
$\displaystyle
ds^{2}=-dT^{2}+\frac{1}{\epsilon^{2}}\left(A(r)\,dr+\sqrt{\epsilon^{2}-1}\,B(r)\,dT\right)^{2}+r^{2}d\phi^{2}$
$\displaystyle
A(r)=\frac{a^{2}\epsilon^{4}+r^{2}}{a^{2}\epsilon^{2}+r^{2}}\quad
B(r)=\frac{\sqrt{\epsilon^{2}(1-\log(A(r)))-1}}{\sqrt{\epsilon^{2}-1}}$ (22) .

$K_{ij}=\sqrt{\epsilon^{2}-1}\begin{pmatrix}\frac{-y^{2}}{r^{3}}&\frac{xy}{r^{3}}\\\\[2.0pt]
\frac{xy}{r^{3}}&\frac{-x^{2}}{r^{3}}\\\
\end{pmatrix}\,,\quad\theta=-\frac{\sqrt{\epsilon^{2}-1}}{r}\,.$ (23)

Closed form expressions for the $K_{ij}$ and $\rho$ can be obtained for the
metric with the warp drive centered at $x_{0}$. However, we are interested in
the limit $x_{0}\gg(x-x_{0})\approx y\approx r_{s}\,$, where the warp drive is
sufficiently far way from the center $r=0$. We shall also take $v_{\infty}\ll
1$ for the Eulerian observers. To $\mathcal{O}(\frac{1}{r})$ and
$\mathcal{O}(v_{\infty})$, the expressions are as follows 888Here,
${}^{(2)}R=0$ on the hypersurfaces, except at $r=0$. This is not true in
$D\geq 4$ dimensions, thereby opening up new possibilities to minimize the
energy density.

$\displaystyle\theta=$
$\displaystyle\,(v+v_{\infty})(x-x_{0})\frac{f^{\prime}(r_{s})}{r_{s}}-(1-f(r_{s}))\,\frac{v_{\infty}}{r}$
$\displaystyle+\frac{y^{2}\,v_{\infty}}{r_{s}r}f^{\prime}(r_{s})$ (24)
$\displaystyle\rho=$
$\displaystyle-\frac{1}{16\pi}\left[\frac{1}{2}(v+v_{\infty})^{2}\,y^{2}\frac{f^{\prime}(r_{s})^{2}}{r_{s}^{2}}\right.$
(25)
$\displaystyle\left.+v\,v_{\infty}(x-x_{0})\left(2\,r_{s}(1-f(r_{s}))-y^{2}f^{\prime}(r_{s})\right)\frac{f^{\prime}(r_{s})}{{r_{s}^{2}\,r}}\right].$

At first, we note that for $v_{\infty}=0$, the expressions falls back to the
usual Alcubierre warp drive in Minkowski spacetime. As expected, the zeroth
order expressions for both $\theta$ and $\rho$ have the combination
$(v+v_{\infty})$ for velocity. In these river flow coordinates, the background
metric is flowing like a fluid along the radially inward direction. However,
the wrap drive is flowing against this “fluid” in the $(+ve)$ $x$–direction,
i.e. it is flowing with velocity $(v+v_{\infty})$ w.r.t. the fluid. Therefore,
it is best to move along the flow to minimize the peak expansions/contractions
and $(-ve)$ energy density. In fact, for a free falling warp drive to zeroth
order they are zero. We also note that, the first
$\mathcal{O}\left(\frac{1}{r}\right)$ term in $\theta$ is that of the
background metric and maybe subtracted, the other one is proportional to
$(y/r)$ and hence very small. By construction $f^{\prime}(r_{s})$ is non zero
only along a very small width of the ring $r_{s}=R$. Hence, the effects of the
warp drive in both $\theta$ and $\rho$ are localized around this ring, the
boundary of the warp bubble.

#### Conical defect backgrounds

Ordinary gravity in three dimensions is special, it is dynamically trivial.
Effects of localized sources are on the global properties whereas the outside
metric is flat. For a point source of mass $M$ placed at the origin $r=0$, the
spatial metric is a cone and the full metric is given by

$ds^{2}=-dt^{2}+\frac{1}{\alpha^{2}}\mskip 1.5mudr^{2}+r^{2}d\phi^{2}\quad
0\leq\phi<2\pi\,,$ (26)

where $\alpha=1-4M$ DESER1984220  999Owing to Birkhoff’s theorem, in $D\geq 4$
dimensions, point mass solutions with spherical symmetry are Schwarzschild
like. We can however break the spherical symmetry in solutions like the cosmic
string or p–branes where the mass density is concentrated along a line or a
membrane and the traverse spatial metric is still conical Kibble:1976sj ;
Randall_1999 .. We would like $\alpha$ to be positive which restricts the mass
$M<(1/4)$. Negative mass or $\alpha>1$ can be considered as unphysical. It is
instructive to scale $t\rightarrow\alpha\mskip 1.5mut$ and write the metric as
MP–charts with the parameter $\epsilon$,

$ds^{2}=-dT^{2}+\frac{1}{\epsilon^{2}}\left(dr+\sqrt{\left(\epsilon^{2}-\alpha^{2}\right)}\,dT\right)^{2}+r^{2}d\phi^{2}\,,$
(27)

such that
$v_{\infty}=\alpha^{2}\sqrt{\left(1-\frac{1}{\epsilon^{2}}\right)\,}$. Keeping
in mind that the charts are invalid for
$\left(\epsilon^{2}-\alpha^{2}\right)<0$, there are two natural choices for
the parameter $\epsilon$ to realize an embedding of the warp drive in these
defect backgrounds.

##### Case $\epsilon=1$

: Though $v_{\infty}$ is zero, the flow velocity of the background metric is
non–zero. However, the spatial metric $\gamma_{ij}=\delta_{ij}$ and
non–singular at the origin $r=0$. The expressions for the $\theta$ and $\rho$
are as follows,

$\displaystyle\theta=$
$\displaystyle\left(v+\sqrt{1-\alpha^{2}}\right)(x-x_{0})\frac{f^{\prime}(r_{s})}{r_{s}}$
(28) $\displaystyle-$
$\displaystyle(1-f(r_{s}))\frac{\sqrt{1-\alpha^{2}}}{r}+\sqrt{1-\alpha^{2}}\,y^{2}\frac{f^{\prime}(r_{s})}{r_{s}r}+\mathcal{O}(\frac{1}{r^{2}})$
$\displaystyle\rho=$
$\displaystyle-\frac{1}{16\pi}\left[\frac{1}{2}\left(v+\sqrt{1-\alpha^{2}}\right)^{2}\frac{f^{\prime}(r_{s})^{2}}{r_{s}^{2}}+\mathcal{O}(\frac{1}{r})\right]\,.$
(29)

Here, the first $\mathcal{O}(\frac{1}{r})$ term in $\theta$ is the background
metric contribution and may be subtracted. Again the $(-ve)$ energy density is
localized around the boundary of the warp bubble, even for the
$\mathcal{O}(\frac{1}{r})$ terms which we omit writing down.

##### Case $\epsilon=\alpha$

: Here, $v_{\infty}=\alpha\sqrt{\left(\alpha^{2}-1\right)}$. This is only
possible if we allow the defect to be sourced by a $(-ve)$ mass $M$. This is a
particularly simple case where the flow velocity of the background metric
vanishes. The expressions for the $\theta$ and $\rho$ are as follows,

$\displaystyle\theta=$ $\displaystyle v\mskip
1.5mu(x-x_{0})\frac{f^{\prime}(r_{s})}{r_{s}}$ (30) $\displaystyle\rho=$
$\displaystyle-\frac{1}{16\pi}\left[\frac{1}{2}v^{2}y^{2}\frac{f^{\prime}(r_{s})^{2}}{\alpha^{2}\,r_{s}^{2}}\right.$
(31) $\displaystyle\left.+\mskip
1.5muv^{2}(x-x_{0})y^{2}(\alpha^{2}-1)\frac{f^{\prime}(r_{s})^{2}}{\alpha^{2}\mskip
1.5mur_{s}^{2}\mskip 1.5mur}+\mathcal{O}(r^{-2})\right]\,.$

We note that the expansion $\theta$ is the same as the original Alcubierre
warp drive whereas the $(-ve)$ energy is suppressed by a factor of
$(1/\alpha^{2})$.

###  3.2 AdS spacetime

As discussed in Faraoni2020WhenPC  the MP–charts for anti–de Sitter space
covers only part of the spacetime. The $D$ dimensional metric in global
coordinates take the form

$ds^{2}=-\left(1+H^{2}r^{2}\right)dt^{2}+\frac{dr^{2}}{1+H^{2}r^{2}}+r^{2}d\Omega^{2}_{(D-2)}$
(32)

where $H=\frac{1}{\ell}$ is the Hubble constant, $\ell$ is the AdS length
scale and $\Lambda=-\frac{(D-1)(D-2)}{2\,\ell^{2}}$ is the cosmological
constant. Following the arguments in section 2 one can construct the MP–charts
with parameter $\epsilon$, however now the Eulerian observers are outgoing and
travels with initial velocity $v_{\infty}$ from $r=0$ with proper time $T$
defined as

$dT=\epsilon\mskip
1.5mudt-\frac{\sqrt{\epsilon^{2}-\left(1+H^{2}r^{2}\right)}}{\left(1+H^{2}r^{2}\right)}dr\,.$
(33)

Our focus is mainly on $\text{AdS}_{3}$ with the line element

$\displaystyle
ds^{2}=-dT^{2}+\frac{1}{\epsilon^{2}}\left(dr-\sqrt{\epsilon^{2}-\left(1+H^{2}r^{2}\right)}dT\right)^{2}+r^{2}d\phi^{2}$
$\displaystyle 0\leq r\leq\frac{\sqrt{\epsilon^{2}-1}}{H}\equiv r_{+}\,,$ (34)

where due to the $(-ve)$ Misner–Sharp–Hernandez mass Misner1964RELATIVISTICEF
; Hernandez:1966zia  of AdS spacetimes the usual $\epsilon=1,\,v_{\infty}=0$,
PG–coordinates doesn’t exist 101010We note that for de–sitter dS spacetimes,
the metric is as in (32) with $H^{2}\rightarrow-H^{2}$ and PG–coordinates
exists along with the MP–charts covering the full spacetime..
Here too, the $T$ = constant hypersurfaces exhibit conical singularity at
$r=0$ as discussed in section 3.1. With the assumption that $H\ll 1$ and
$\epsilon$ is some reasonable number greater than 1, we can push the charts to
cover large regions of spacetime except for regions near the AdS boundary
$r\rightarrow\infty$. The extrinsic curvature tensor away from the singularity
is as follows,

$\displaystyle K_{ij}=-$ (35)
$\displaystyle\begin{pmatrix}\frac{H^{2}r^{2}\left(x^{2}+y^{2}\epsilon^{2}\right)-y^{2}\epsilon^{2}\left(\epsilon^{2}-1\right)}{r^{3}\epsilon^{2}\sqrt{\epsilon^{2}-\left(1+H^{2}r^{2}\right)}}&\frac{xy\left(\epsilon^{2}-1\right)\left(\epsilon^{2}-H^{2}r^{2}\right)}{r^{3}\epsilon^{2}\sqrt{\epsilon^{2}-\left(1+H^{2}r^{2}\right)}}\\\\[3.0pt]
\frac{xy\left(\epsilon^{2}-1\right)\left(\epsilon^{2}-H^{2}r^{2}\right)}{r^{3}\epsilon^{2}\sqrt{\epsilon^{2}-\left(1+H^{2}r^{2}\right)}}&\frac{H^{2}r^{2}\left(x^{2}\epsilon^{2}+y^{2}\right)-x^{2}\epsilon^{2}\left(\epsilon^{2}-1\right)}{r^{3}\epsilon^{2}\sqrt{\epsilon^{2}-\left(1+H^{2}r^{2}\right)}}\\\
\end{pmatrix}$
$\displaystyle\theta=\frac{\epsilon^{2}-\left(1+2H^{2}r^{2}\right)}{r\sqrt{\epsilon^{2}-\left(1+H^{2}r^{2}\right)}}\,.$
(36)

We shall now embed the warp drive centered at $x_{0}$ in this background.
Closed form expressions for the $K_{ij}$ and $\rho$ can be obtained but as in
the Minkowski case, we are mostly interested in the limit
$x_{0}\gg(x-x_{0})\approx y\approx r_{s}\,$, where the warp drive is
sufficiently far way from the singularity at the center $r=0$. As we have
finite range charts, we have to treat the order of limits carefully. The
radial distance $r$ in these MP–charts is cutoff at $r_{+}(\epsilon,H)\,$
therefore it is best to scale it as $r=(1-s_{r})\mskip 1.5mur_{+}$ such that
$0\leq s_{r}\leq 1$ and then focus on the region near $s_{r}\approx 0\,$. To
$\mathcal{O}(s_{r})$ and $\mathcal{O}(H)$, the expansion and energy density
are as follows,

$\displaystyle\theta=$ $\displaystyle\left(v-\sqrt{2(\epsilon^{2}-1)\mskip
1.5mus_{r}}\right)(x-x_{0})\frac{f^{\prime}(r_{s})}{r_{s}}$ (37)
$\displaystyle-\frac{(1-f(r_{s})}{\sqrt{2}}\left(\frac{1}{\sqrt{s_{r}}}-\frac{11}{4}\sqrt{s_{r}}\right)H$
$\displaystyle-\sqrt{2}\mskip
1.5muy^{2}\sqrt{s_{r}}\frac{f^{\prime}(r_{s})}{r_{s}}\mskip 1.5muH$
$\displaystyle\rho=$
$\displaystyle-\frac{1}{16\pi}\left[\frac{1}{2}\left(v-\sqrt{2\mskip
1.5mu(\epsilon^{2}-1)s_{r}}\right)^{2}y^{2}\frac{f^{\prime}(r_{s})^{2}}{r_{s}^{2}\epsilon^{2}}\right.$
(38)
$\displaystyle+v^{2}y^{2}(x-x_{0})\sqrt{(\epsilon^{2}-1)}\frac{f^{\prime}(r_{s})^{2}}{r_{s}^{2}\epsilon^{2}}H$
$\displaystyle\left.+\mathcal{O}(\sqrt{s_{r}})\mskip
1.5muH+\mathcal{O}(H^{2})\vphantom{\frac{1}{2}}\right]\,.$

Again, we notice the combination $\left(v-\sqrt{2\mskip
1.5mu(\epsilon^{2}-1)\mskip 1.5mus_{r}}\right)$ appearing at the zeroth order,
however for AdS the background flow velocity and the warp drive velocity are
both along $(+ve)$ $x$–direction. The first $\mathcal{O}(H)$ term in (37) is
the contribution from the background AdS and maybe subtracted and the second
term is very small. Similarly, in (38) the background contribution from the
cosmological constant has been subtracted in the $\mathcal{O}(H^{2})$ terms.
After the appropriate subtractions remaining terms are again proportional to
$f^{\prime}(r_{s})^{\\#}$ such that the change in $\theta$ and $\rho$ are
localized around the boundary of the warp bubble. Interestingly, the zeroth
order term in the $(-ve)$ energy density is suppressed by $(1/\epsilon^{2})$
thereby reducing the energy as $\epsilon>1$ 111111This is also true for the
Minkowski space. It is not visible in (3.1) as we truncated the expressions to
$\mathcal{O}(v_{\infty})$..
The analysis for warp drive embedded in dS spacetime is similar. However, the
MP–charts now cover the full spacetime. Even crossing the apparent horizon at
$r=(1/H)$ is regular. To study the warp drive near this horizon, we can use
the scaling $r=((1-s_{r})/H)$ and follow the steps as discussed above. the
results are similar and as expected, nothing special happens because of the
horizon.

#### Conical defect backgrounds

Similar to the case of conical defect solutions in three dimensional Minkowski
space, $\text{AdS}_{3}$ too have conical defect solutions for a point mass $M$
placed at the origin $r=0$. The spatial metric is a cone and the full metric
is given by

$ds^{2}=-(\alpha^{2}+r^{2}H^{2})dt^{2}+\frac{dr^{2}}{(\alpha^{2}+r^{2}H^{2})}+r^{2}d\phi^{2}\,,$
(39)

where $\alpha=1-4M$ and $0\leq\phi<2\pi$ DESER1984405  121212 $\text{dS}_{3}$
is bit tricky. In all three spacetimes, there is a second source at
$r=\infty$. In Minkowski and $\text{AdS}_{3}$ spacetimes the geodesic distance
between the source at $r=0$ and $r=\infty$ is infinite, so the second source
can be ignored but for $\text{dS}_{3}$ it is finite and contributes.. The
MP–charts with $\epsilon=\alpha$ doesn’t exist. Therefore, we focus on the
$\epsilon=1,\,v_{\infty}=0$ charts given by

$\displaystyle
ds^{2}=-dT^{2}+\left(dr-\sqrt{(1-\alpha^{2})-r^{2}H^{2}}\,dT\right)^{2}+r^{2}d\phi^{2}$
$\displaystyle 0\leq r\leq\sqrt{1-\alpha^{2}}\,H^{-1}\equiv{}^{c}r_{+}\,.$
(40)

We note that for physical particles, we have $0<\alpha\leq 1$. We shall now
embed the warp drive in this conical $\text{AdS}_{3}$ background and do a
similar analysis as done earlier. The important changes are that the spatial
metric is now $\delta_{ij}$ and scaling of the radial coordinate is
$r=(1-s_{r})\,{}^{c}r_{+}$ such that $0\leq s_{r}\leq 1$ and then focus on the
region $s_{r}\approx 0$. The expansion $\theta$ is as given in (37) with the
replacement $(\epsilon^{2}-1)\rightarrow(1-\alpha^{2})$. Parts of the
$\mathcal{O}(H)$ not proportional to $f^{\prime}(r_{s})$ are from the
background metric and should be subtracted. For the energy density $\rho$, the
zeroth order $\mathcal{O}(H^{0})$ terms is the same as in (38) with the
replacement $(\epsilon^{2}-1)\rightarrow(1-\alpha^{2})$ and omitting the
$\epsilon^{2}$ in the denominator. From $\mathcal{O}(H)$, terms are different
starting at $\mathcal{O}\left(\sqrt{s_{r}}\right)$. Again the background
contribution from the cosmological constant needs to be subtracted from the
$\mathcal{O}(H^{2})$ terms. Throughout, we also don’t have the
$(1/\epsilon^{2})$ suppression as $\epsilon=1$.

##  4 Discussion

In this paper we have extended the Alcubierre warp drive spacetime to an
infinite class of warp drives utilizing the Martel–Poisson charts for
Minkowski and AdS spacetimes and their conical defect cousins in three
dimensions. The Minkowski spacetime admits Painlevé–Gullstrand coordinates but
for AdS spacetime the usual PG charts vanish and the warp drive spacetimes
discussed in this paper are the only possibility. Though we have analyzed the
warp drives away from the conical singularities at the origin, the fact that
spatial hypersurfaces are cones have global consequences which have been
captured in the expressions for the expansion and $(-ve)$ energy densities
localized around the warp drive boundary.
There are several possibilities for generalizations. We can look for warped
spherically symmetric background spacetimes

$ds^{2}=-e^{2\Phi}dt^{2}+\frac{dr^{2}}{f}+r^{2}d\Omega_{(D-2)}^{2}$ (41)

which admits MP–charts with partial coverage of the spacetime,

$\displaystyle ds^{2}=$ $\displaystyle-
dT^{2}+\frac{e^{2\Phi}}{\epsilon^{2}f}\left(dr+\frac{\epsilon\sqrt{f}}{e^{\Phi}}\sqrt{1-\frac{e^{2\Phi}}{\epsilon^{2}}}\mskip
1.5mudT\right)^{2}$
$\displaystyle+r^{2}d\Omega_{(D-2)}^{2}\,,\quad\left(1-\frac{e^{2\Phi}}{\epsilon^{2}}\right)\geq
0\,.$ (42)

Another direction would be to look at higher dimensions $D\geq 4$ and find
spacetimes whose MP–charts has spatial hypersurfaces which near the
singularity at the origin takes the form,

$\displaystyle ds_{C}^{2}=$ $\displaystyle
d\tilde{r}^{2}+\tilde{r}^{2}d\tilde{\phi}^{2}$ (43)
$\displaystyle+\sum_{i,j=1}^{D-3}\left(w_{ij}(\theta)+h_{ij}(\theta)\tilde{r}^{2}\right)d\theta^{i}d\theta^{j}+\mathcal{O}(\tilde{r}^{4})$

or generalize the definition of a cone from two to higher dimensions,

$ds_{C}^{2}=d\tilde{r}^{2}+\tilde{r}^{2}d\tilde{\Omega}_{(D-2)}^{2}\,,$ (44)

where unless $d\tilde{\Omega}_{(D-2)}^{2}$ is the metric of the unit sphere
$S^{(D-2)}$, the generalized cone $ds_{C}^{2}$ is singular. For example,
MP–charts of various black hole backgrounds would come under its purview. For
most of these generalizations, the Ricci scalar of the spatial hypersurface
will be non zero thereby opening up new possibilities for the energy density
$\rho$. If we focus on the spherically symmetric cases, the expressions
involving the extrinsic curvatures would generalize easily from the three
dimensional expressions in the previous sections upto some coefficients. We
have to symmetrize the expressions in all spatial coordinates except $x$, the
direction of travel for the warp drive.
While warp drives are beyond the realm of current experiments, various analog
gravity setups imitate the metric Fischer_2003 ; finazzi2012analogue  and can
be studied experimentally PhysRevLett.105.240401 . It might be possible to
generalize these setups to accommodate the Martel–Poisson versions of the warp
drive.
Various other properties of the wrap drive are of interest, for example the
zero–expansion drives of Natário Nat_rio_2002 , the formation of event horizon
and various quantum effects Hiscock_1997 ; Finazzi_2009 ; Barcel__2022  to
optimization of various parameters including the form functions and the idea
of a warp bubble inside another bubble Broeck_1999 ; Helmerich_2023 . We hope
to study some of these aspects in our future works.

###### Acknowledgements.

The work A.C. is supported by IIT Bhubaneswar Seed Grant SP–103.

##  Appendix A Different choices for $r_{s}$ and $\gamma_{ij}$

In this appendix we shall explore, two combinations of choices for the warp
drive radial coordinate $r_{s}$ and the spatial induced metric $\gamma_{ij}$
as mentioned in section 3. Our focus is on the warp drives embedded in the
MP–charts of Minkowski spacetime but similar analysis goes through for the
conical defect background, $\text{AdS}_{3}$ & $\text{dS}_{3}$ spacetimes and
their conical defects cousins.

##### Case I :

We take $r_{s}=\sqrt{(x-x_{0})^{2}+y^{2}}$ and
$\gamma_{ij}=\delta_{ij}+(1-f(r_{s}))\left(\frac{1}{\epsilon^{2}}-1\right)\frac{x^{i}x^{j}}{r^{2}}$.
There are two important deviations from the expressions in (3.1). Firstly, the
effect of $(1-f(r_{s}))$ term in $\gamma_{ij}$ starts to show only at
$\mathcal{O}(v_{\infty}^{2})$, so if we are truncating to
$\mathcal{O}(v_{\infty})$, all results will match. Secondly, the spatial
curvature tensors, particularly the Ricci scalar ${}^{(2)}R$ is not zero at
finite $r$. Its contribution to the energy density $\rho$ starts at order
$\mathcal{O}(v_{\infty}^{2})$,

$\frac{(y^{2}-r_{s}^{2})f^{\prime}(r_{s})-r_{s}y^{2}f^{\prime\prime}(r_{s})}{r_{s}^{3}}v_{\infty}^{2}+\mathcal{O}(\frac{1}{r})$
(45)

Similarly, the $K_{ij}$ in (20) now picks up contributions from the
$\frac{1}{2}\partial_{T}\gamma_{ij}$ term, again modifying expressions from
$\mathcal{O}(v_{\infty}^{2})$ onward. In all, this choice complicates things
but after background subtraction $\theta$ and $\rho$ are proportional to
$f^{\prime}(r_{s})^{\\#}$ and $f^{\prime\prime}(r_{s})^{\\#}$, therefore
localized around the wrap drive bubble. A lot now depends on the choice of
$f(r_{s})$, tweaking it might give a $(+ve)$ ${}^{(2)}R$ thereby reducing the
overall $(-ve)$ energy density required to operate the drive.

##### Case II :

We take $r_{s}=\sqrt{\gamma_{ij}(x-x_{o}(T))^{i}(x-x_{0}(T))^{j}}$ and
$\gamma_{ij}=\delta_{ij}+\left(\frac{1}{\epsilon^{2}}-1\right)\frac{x^{i}x^{j}}{r^{2}}$.
Here, as we are not changing the spatial metric, there is no contribution from
${}^{(2)}R$, neither $\frac{1}{2}\partial_{T}\gamma_{ij}$ contributes to
$K_{ij}$. Extra contributions as compared to results in section 3 are due to
the changes in the derivatives of $f(r_{s})$ or equivalently $r_{s}$ w.r.t.
$(x,\,y)$ coordinates. Again, the changes are from
$\mathcal{O}(v_{\infty}^{2})$ onward and yet again the contributions to
$\theta$ and $\rho$ are localized to the boundary of the warp drive.

## References

  * (1) M.S. Morris, K.S. Thorne, American Journal of Physics 56(5), 395 (1988).  DOI 10.1119/1.15620.  URL https://doi.org/10.1119/1.15620
  * (2) M. Alcubierre, Classical and Quantum Gravity 11(5), L73–L77 (1994).  DOI 10.1088/0264-9381/11/5/001.  URL http://dx.doi.org/10.1088/0264-9381/11/5/001
  * (3) J. Natário, Classical and Quantum Gravity 19(6), 1157–1165 (2002).  DOI 10.1088/0264-9381/19/6/308.  URL http://dx.doi.org/10.1088/0264-9381/19/6/308
  * (4) F.S.N. Lobo, M. Visser, Classical and Quantum Gravity 21(24), 5871–5892 (2004).  DOI 10.1088/0264-9381/21/24/011.  URL http://dx.doi.org/10.1088/0264-9381/21/24/011
  * (5) A.E. Everett, T.A. Roman, Physical Review D 56(4), 2100–2108 (1997).  DOI 10.1103/physrevd.56.2100.  URL http://dx.doi.org/10.1103/PhysRevD.56.2100
  * (6) K. Martel, E. Poisson, American Journal of Physics 69(4), 476–480 (2001).  DOI 10.1119/1.1336836.  URL http://dx.doi.org/10.1119/1.1336836
  * (7) U.R. Fischer, M. Visser, Europhysics Letters (EPL) 62(1), 1–7 (2003).  DOI 10.1209/epl/i2003-00103-6.  URL http://dx.doi.org/10.1209/epl/i2003-00103-6
  * (8) S. Finazzi.  Analogue gravitational phenomena in bose-einstein condensates (2012)
  * (9) O. Lahav, A. Itah, A. Blumkin, C. Gordon, S. Rinott, A. Zayats, J. Steinhauer, Phys. Rev. Lett. 105, 240401 (2010).  DOI 10.1103/PhysRevLett.105.240401.  URL https://link.aps.org/doi/10.1103/PhysRevLett.105.240401
  * (10) V. Faraoni, G. Vachon, The European Physical Journal. C, Particles and Fields 80 (2020).  URL https://api.semanticscholar.org/CorpusID:219956329
  * (11) R. Garattini, K. Zatrimaylov.  Black holes and warp drive (2023)
  * (12) R. Garattini, K. Zatrimaylov.  On the wormhole–warp drive correspondence (2024)
  * (13) D.V. Fursaev, S.N. Solodukhin, Physical Review D 52(4), 2133–2143 (1995).  DOI 10.1103/physrevd.52.2133.  URL http://dx.doi.org/10.1103/PhysRevD.52.2133
  * (14) H.G. Ellis.  The warp drive and antigravity (2004)
  * (15) C.V.D. Broeck, Classical and Quantum Gravity 16(12), 3973–3979 (1999).  DOI 10.1088/0264-9381/16/12/314.  URL http://dx.doi.org/10.1088/0264-9381/16/12/314
  * (16) S. Deser, R. Jackiw, G. ’t Hooft, Annals of Physics 152(1), 220 (1984).  DOI https://doi.org/10.1016/0003-4916(84)90085-X.  URL https://www.sciencedirect.com/science/article/pii/000349168490085X
  * (17) T.W.B. Kibble, J. Phys. A 9, 1387 (1976).  DOI 10.1088/0305-4470/9/8/029
  * (18) L. Randall, R. Sundrum, Physical Review Letters 83(17), 3370–3373 (1999).  DOI 10.1103/physrevlett.83.3370.  URL http://dx.doi.org/10.1103/PhysRevLett.83.3370
  * (19) C.W. Misner, D.H. Sharp, Physical Review 136, 571 (1964).  URL https://api.semanticscholar.org/CorpusID:123043032
  * (20) W.C. Hernandez, C.W. Misner, Astrophys. J. 143, 452 (1966).  DOI 10.1086/148525
  * (21) S. Deser, R. Jackiw, Annals of Physics 153(2), 405 (1984).  DOI https://doi.org/10.1016/0003-4916(84)90025-3.  URL https://www.sciencedirect.com/science/article/pii/0003491684900253
  * (22) W.A. Hiscock, Classical and Quantum Gravity 14(11), L183 (1997).  DOI 10.1088/0264-9381/14/11/002.  URL https://dx.doi.org/10.1088/0264-9381/14/11/002
  * (23) S. Finazzi, S. Liberati, C. Barceló, Physical Review D 79(12) (2009).  DOI 10.1103/physrevd.79.124017.  URL http://dx.doi.org/10.1103/PhysRevD.79.124017
  * (24) C. Barceló, V. Boyanov, L.J. Garay, E. Martín-Martínez, J.M. Sánchez Velázquez, Journal of High Energy Physics 2022(8) (2022).  DOI 10.1007/jhep08(2022)288.  URL http://dx.doi.org/10.1007/JHEP08(2022)288
  * (25) C. Helmerich, J. Fuchs, A. Bobrick, L. Sellers, S. Dangelo, G. Martire, J.F. Agnew, in _AIAA SCITECH 2023 Forum_ (American Institute of Aeronautics and Astronautics, 2023).  DOI 10.2514/6.2023-0553.  URL http://dx.doi.org/10.2514/6.2023-0553


Article Title]Geostatistics in the presence of multivariate complexities: comparison of multi-Gaussian transforms

[1,2]Sultan<EMAIL_ADDRESS>
1,2]Peter A. Dowd

1]Chaoshui Xu

[1]ARC Training Centre for Integrated Operations for Complex Resources, University of Adelaide, Adelaide, Australia

[2]School of Civil, Environmental & Mining Engineering, University of Adelaide, Adelaide, Australia

One of the most challenging aspects of multivariate geostatistics is dealing with complex relationships between variables. Geostatistical co-simulation and spatial decorrelation methods, commonly used for modelling multiple variables, are ineffective in the presence of multivariate complexities. On the other hand, multi-Gaussian transforms are designed to deal with complex multivariate relationships, such as non-linearity, heteroscedasticity and geological constraints. These methods transform the variables into independent multi-Gaussian factors that can be individually simulated. This study compares the performance of the following multi-Gaussian transforms: rotation based iterative Gaussianisation, projection pursuit multivariate transform and flow transformation. Case studies with bivariate complexities are used to evaluate and compare the realisations of the transformed values. For this purpose, commonly used geostatistical validation metrics are applied, including multivariate normality tests, reproduction of bivariate relationships, and histogram and variogram validation. Based on most of the metrics, all three methods produced results of similar quality. The most obvious difference is the execution speed for forward and back transformation, for which flow transformation is much slower.

§ INTRODUCTION

Geostatistical conditional simulation often requires the modelling of multiple cross-correlated variables. For example, multivariate geostatistics is commonly applied to model ore grades and other mineral deposit variables [Journel and Huijbregts, 1978, Wackernagel, 2003]. Among the Gaussian methods, there are three ways to perform this task [Rossi and Deutsch, 2014]. The first is to use conditional co-simulation algorithms such as sequential Gaussian co-simulation [Verly, 1993] and turning bands co-simulation [Emery, 2008]. The second approach is to use hierarchical co-simulation, which enables full, collocated and multi-collocated co-kriging to be applied [Almeida and Journel, 1994]. Finally, the third way is to transform variables into uncorrelated factors and simulate them individually.

Direct and hierarchical co-simulations use a linear model of co-regionalisation (LMC) [Journel and Huijbregts, 1978] or a Markov model [Journel, 1999] to account for cross-correlations between variables. However, inference of cross-variograms can be challenging, whereas transformation to independent factors significantly simplifies multivariate modelling. For example, principal component analysis (PCA) and independent component analysis (ICA) have been applied in geostatistics to transform cross-correlated variables into independent orthogonal factors [Davis and Greenes, 1983, Tercan and Sohrabian, 2013]. However, the validity of extending PCA and ICA decorrelation to non-zero lag distances is an assumption. Alternatively, minimum/maximum autocorrelation factors (MAF) performs a double spectral decomposition at lag zero and a single non-zero lag to achieve spatial decorrelation [Desbarats and Dimitrakopoulos, 2000]. Other spatial decorrelation methods applied in geostatistics are uniformly weighted exhaustive diagonalisation with Gauss iterations and rotational joint diagonalisation [Mueller and Ferreira, 2012, Mueller et al, 2020].

Geostatistical co-simulation is also limited by the multi-Gaussianity assumption, which is often unrealistic as there may be complex multivariate relationships [Barnett et al, 2014]. LMC cannot include these complexities in a covariance matrix, and linear transformations (e.g., PCA, ICA and MAF) are only applicable to linear relationships. This motivated [Leuangthong and Deutsch, 2003] to apply stepwise conditional transformation (SCT) in the geostatistical framework to deal with multivariate complexities, such as non-linearity, heteroscedasticity and inequality constraints. In geostatistics, SCT was the earliest application of a multi-Gaussian transform (MGT). The idea behind MGT approaches is to simply transform variables into standard multivariate Gaussian distributions with zero correlation. SCT removes the complexities by transforming the original variables into multi-Gaussian factors. However, SCT cannot handle high-dimensional datasets and it requires extensive data cleaning and the order of the variables to be predefined [Rossi and Deutsch, 2014]. Recently, [de Figueiredo et al, 2021] proposed a direct multivariate simulation based on SCT and demonstrated its applicability to a six-dimensional dataset with non-linear relationships.

Two of the most popular MGT methods in geostatistics are the projection pursuit multivariate transform (PPMT) [Barnett et al, 2014, Barnett et al, 2016] and flow transformation, also known as flow anamorphosis (FA) [van den Boogaart et al, 2017]. PPMT searches for the direction that has the maximum projection index [Friedman, 1987] and applies a normal score transformation in that direction. It is an iterative algorithm that transforms the data along projections based on the departure from Gaussianity. PPMT can be applied to a higher number of variables than SCT, it can work with smaller datasets and does not need the order of variables to be predefined [Barnett et al, 2014]. On the other hand, FA continuously deforms the original distribution into multi-Gaussian space using Lagrangian mechanics [van den Boogaart et al, 2017]. Furthermore, its affine equivariance makes FA suitable for compositional data analysis [Tolosana-Delgado et al, 2019], in which it can be paired with various log-ratio transforms [Pawlowsky-Glahn et al, 2015].

A case study of the Gol-e-Gohar iron deposit by [Hosseini and Asghari, 2019] suggests that the combination of log-ratio+FA+MAF produces fewer artefacts during back transformations and provides better reproduction of compositional constraints than the combination of PPMT+MAF. In the case study, the additive log-ratio transformation was chained only with FA because of its affine equivariance. However, [Manchuk et al, 2017] successfully chained PPMT with an isometric log-ratio transformation [Egozcue et al, 2003] to reproduce the sum constraint. Similarly, additive and fraction ratios (i.e., not logarithmic) have been used on PPMT factors to reproduce sum and fractional constraints, also known as inequality constraints [Bassani et al, 2018]. Nevertheless, it is also essential to chain these methods with MAF or any other spatial decorrelation to ensure that factors are independent at non-zero lags. For example, the combination of PPMT+MAF shows significantly better variogram reproduction compared to the transformation by PPMT or by MAF [Erten and Deutsch, 2021].

This study compared the performance of PPMT, FA and a relatively new method to geostatistics, rotation based iterative Gaussianisation (RBIG) [Laparra et al, 2011]. RBIG is similar to PPMT, but it rotates the data using either PCA or ICA and applies the normal score transformation after each rotation. The following sections provide more details on selected methods together with a comprehensive comparison based on different metrics. The comparison is based on three bivariate case studies from undisclosed mining deposits with strong multivariate complexities. The results were carefully assessed using different statistical and qualitative metrics.

§ MATERIALS AND METHODS

§.§ Case studies with multivariate complexities

This paper applies selected methods to three bivariate case studies with complex relationships (Figure <ref>). These are confidential mining datasets and are, therefore, undisclosed. Case A consists of 6,074 drill hole samples for which there is an inequality constraint between total and soluble copper grades. The average horizontal spacing between samples is 123 m with a composite length of 2 m. This type of multivariate complexity occurs when one variable is a fraction of another variable, and there have been attempts to model this relationship [Hosseini and Asghari, 2015, Bassani et al, 2018, Abildin et al, 2019]. Inequality constraints can also occur in iron ore deposits between iron and elements such as silica and aluminium oxide when there is a linear inequality between variables [Madani and Abulkhair, 2020, Abulkhair and Madani, 2021].

Cross plots of three bivariate case studies with an inequality constraint (case A), non-linearity (case B) and heteroscedasticity (case C). Number of samples, Pearson $\rho$ and Spearman $\rho_S$ correlation coefficients are shown in text boxes.

For case B, there is a non-linearity between 9,990 iron and magnesium oxide samples. In this case, the sample points are irregularly spaced with an average of 4.2 m for horizontal and 2.4 m for vertical spacing. Non-linearities are not rare in drill hole datasets, and several geostatistical case studies have focussed on modelling such datasets [Leuangthong and Deutsch, 2003, Barnett et al, 2014, Barnett et al, 2016, de Figueiredo et al, 2021, Erten and Deutsch, 2021].

Finally, the case C data come from an underground mine with an average of 14.5 m spacing and a composite length of 2 m. It is a bivariate dataset with 33,021 samples of titanium and zirconium with a heteroscedastic relationship between them. Heteroscedasticity is similar to a linear relationship, expressed by a non-constant variance of one variable across the range of values of another. Heteroscedastic relationships has been reported in many multivariate geostatistical studies [Barnett et al, 2014, Barnett et al, 2016, de Figueiredo et al, 2021, Erten and Deutsch, 2021]. However, case C demonstrates a more obvious heteroscedastic relationship with a moderately high correlation.

As all three cases are mining datasets prone to sampling irregularities, duplicates and outliers, they underwent careful data cleaning. In addition, cell declustering [Deutsch and Journel, 1992] was applied for cases A and B. Cells of 85m × 85m × 85m were chosen for case A and 18m × 18m × 18m for case B after evaluating the effect of different cell sizes. As a result, mean and standard deviation values decreased after declustering for both cases. Cell declustering was not performed for case C because the sum of weights was not equal to the number of samples. Table <ref> shows the descriptive statistics for all three case studies. As multivariate complexities are involved, both Spearman and Pearson correlation coefficients are reported.

Descriptive statistics of three bivariate datasets

2*Parameter     2cCase A     2cCase B     2cCase B

CuT (%)     CuS (%)     Fe (%)     MgO (%)     Ti (ppm)     Zr (ppm)

Mean     0.42     0.12     29.01     8.86     2,928.88     128.75

Standard deviation     0.64     0.39     16.10     10.60     2,357.84     64.01

Pearson correlation     2c0.71     2c-0.86     2c0.68

Spearman correlation     2c0.77     2c-0.90     2c0.86

§.§ Multivariate transformation

§.§.§ Rotation based iterative Gaussianisation

RBIG is an iterative algorithm that applies marginal Gaussianisation followed by an orthonormal rotation. Although [Laparra et al, 2011] demonstrated that any orthonormal rotation, even a simple random rotation matrix, could be used in RBIG, PCA and ICA are more suitable. The choice of a rotation matrix is important and is usually based on multiple factors. PCA provides a suboptimal convergence rate compared to ICA and requires more iterations. However, convergence in ICA takes more time, especially in higher dimensional cases. In this study, we used RBIG with both PCA (RBIGP) and ICA (RBIGI) rotations for comparison. For RBIGP, the original MATLAB code from [Laparra et al, 2011] was implemented with histogram equalisation for marginal Gaussianisation. For RBIGI, a Python implementation was used with a normal score transformation for marginal Gaussianisation together with fast ICA [Hyvarinen, 1999] from the Scikit-learn library [Pedregosa et al, 2011] for rotation.

The steps in RBIG are as follows:

* First marginal Gaussianisation

\begin{equation}
\end{equation}

where $\Psi_{(0)}(Z)$ is a normal score transformation for RBIGI or histogram equalisation for RBIGP applied to each dimension of the original data $Z$ before the first iteration.

* ICA or PCA rotation expressed by an orthonormal rotation matrix $R_{(i)}$ at each iteration $i$

\begin{equation}
\end{equation}

* Marginal Gaussianisation of rotated variables at each iteration $i$

\begin{equation}
\end{equation}

* Repeat steps 2 and 3 for a predefined number of iterations.

In this study, the algorithm runs for a fixed number of iterations without any stopping criteria. In addition, RBIG saves rotation matrices and Gaussian tables at each iteration, including the first marginal Gaussianisation. Finally, the back transformation to the original state is performed in reverse order of the data saved at each iteration.

§.§.§ Projection pursuit multivariate transform

PPMT methodology is based on iteratively searching for interesting projections followed by normal score transformation along those projections [Barnett et al, 2014, Barnett et al, 2016]. An interesting projection is one that has a maximum departure from Gaussianity based on the projection index developed by [Friedman, 1987]. The input parameters for the original algorithm are a maximum number of iterations and stopping criteria. The algorithm terminates after reaching the targeted projection index percentile based on the bootstrapping algorithm or after a set number of iterations [Barnett et al, 2014]. A detailed description and visualization of PPMT and its stopping criteria can be found in [Barnett et al, 2016].

Similar to RBIG, PPMT runs for a fixed number of iterations in this study. The purpose of doing so is to check the multivariate normality of the transformed RBIG and PPMT factors after the same number of iterations. PPMT records rotation matrices and Gaussian tables at each step, which are used during the back transformation following the above steps in reverse order. We used a Python implementation of PPMT based on the original Fortran code from [Barnett et al, 2014].

§.§.§ Flow transformation

The FA methodology differs from that of RBIG and PPMT, which are based on orthonormal rotations and marginal Gaussianisations. FA continuously deforms the original kernel density function into standard multi-Gaussian space using Lagrangian mechanics [van den Boogaart et al, 2017]. The two main input parameters that characterise this deformation are the starting $\sigma_0$ and the final $\sigma_1$ spreads of the kernel that control the smoothing. This means that $\sigma_0$ controls how strongly FA deforms the kernels, for which a smaller value results in more Gaussian factors. In contrast, $\sigma_1$ controls the ranges of the produced factors, and it is recommended that $\sigma_1=\sigma_0+1$ so that the marginal distributions of the transformed data have standard deviations close to 1 [Talebi et al, 2019]. For the details of the methodology of the FA algorithm, readers are referred to the original literature [van den Boogaart et al, 2017, Tolosana-Delgado and Mueller, 2021].

In this study, the “gmGeostats” CRAN package [Tolosana-Delgado and Mueller, 2021] was used for FA and the transformation was chained twice to achieve multivariate normality.

§.§.§ Chained multi-Gaussian transform and spatial decorrelation

MGT methods, similar to PCA and ICA linear transforms, can only guarantee a decorrelation at lag zero. A practical solution to ensure spatial decorrelation is to chain those methods with MAF [Desbarats and Dimitrakopoulos, 2000]. The original geostatistical application of MAF uses PCA twice, once on a covariance matrix at lag 0 and once on a single non-zero lag of the cross-variogram function. However, as RBIG, PPMT or FA already produce independent factors, only the second MAF is used. Chained MGT and MAF comprise the following steps:

* Transform the original data $Z$ into independent standard multi-Gaussian variables $Y^{MG}$.

* Calculate a sphering matrix $S^{-1/2}$ for multi-Gaussian factors:

\begin{equation}
\end{equation}

where $Q$ is an eigenvector matrix and $\Lambda$ is the corresponding diagonal matrix.

* Compute an eigenvector matrix $Q_h$ from the spectral decomposition of the cross-variogram matrix at lag $h$.

* Multiply the multi-Gaussian factors by the sphering and eigenvector matrices

\begin{equation}
\end{equation}

The resulting factors will also be spatially independent. Back transformation is performed by using transposed sphering and eigenvector matrices. A combination of MAF with MGT has been applied in various geostatistical case studies [Hosseini and Asghari, 2019, Erten and Deutsch, 2021, Tolosana-Delgado and Mueller, 2021]. However, MAF is applied on a single non-zero lag distance, so its spatial decorrelation is not perfect at all lags.

§.§ Metrics for comparison

In this study, RBIGP, RBIGI, PPMT and FA were compared using six metrics: multivariate normality, spatial decorrelation, execution times, qualitative assessment of the reproduction of multivariate distributions, and histogram and variogram validation. Multivariate normality (MVN) can be evaluated using different multivariate normality tests from the MVN R package [Korkmaz et al, 2014]. The MVN tests in this package include Mardia’s measures of multivariate skewness and kurtosis [Mardia, 1970], Royston’s techniques for assessing multivariate normality [Royston, 1983], the Henze-Zirkler invariant consistent tests for multivariate normality [Henze and Zirkler, 1990] and Energy statistics [Székely and Rizzo, 2013]. Many studies in geostatistics have used MVN tests to assess the performance of MGT e,g., [van den Boogaart et al, 2017, Tolosana-Delgado et al, 2019, Tolosana-Delgado and Mueller, 2021]. In this study, only Henze-Zirkler’s and Energy tests were used based on their consistency and robustness, confirmed in some comparison and review articles [Joenssen and Vogel, 2014, Ebner and Henze, 2020].

A common problem with MGT methods is that they cannot ensure spatial decorrelation. In a geostatistical context, spatial decorrelation is sometimes more crucial than multivariate normality. It is, therefore, important to assess spatial decorrelation even after applying MAF. For this purpose, the quality of the spatial decorrelation was evaluated through experimental cross-variograms. In addition, a quantitative assessment of spatial decorrelation was conducted using the relative deviation from diagonality $\tau(h)$ together with the spatial diagonalisation efficiency $\kappa(h)$ measures suggested by [Tercan, 1999]. The first measure compares the absolute sum of off-diagonal elements in the factor variogram matrix to the corresponding diagonal elements. The second one compares the sum of squares of off-diagonal elements in the factor variogram matrix to the sum of squares of off-diagonal elements in the sample variogram matrix. Perfect spatial decorrelation will result in zero for $\tau(h)$ and one for $\kappa(h)$.

The execution time for both forward and inverse multi-Gaussian transformations is another important factor. RBIGI, PPMT and FA were applied in a Jupyter Notebook environment, where their CPU times could be recorded. FA in the gmGeostats R package [Tolosana-Delgado and Mueller, 2021] and the Python implementations of RBIGI and PPMT were used for the work reported in this paper. Although the original Fortran program for PPMT is much faster, the Python code may provide a fairer comparison of PPMT and RBIGI. For RBIGP, the original MATLAB code [Laparra et al, 2011] was used, in which CPU time can also be measured.

Finally, reproduction of histograms and variograms are the most critical geostatistical properties. Even though their reproduction can be assessed qualitatively, the root mean square error (RMSE) was used to obtain a quantitative comparison of results. For example, the RMSE between a hundred percentiles of original and simulated cumulative distribution functions (CDFs) were calculated for histogram validation. Similarly, RMSE measures were calculated between experimental and simulated direct and cross-variograms at multiple lag distances. Metrics such as these can be found in other geostatistical studies [Mueller and Ferreira, 2012, Erten and Deutsch, 2021]. However, it is more difficult to provide quantitative comparisons of the reproduction of multivariate complexities, which is the primary objective of this study. For this purpose, cross-plots between simulated variables were qualitatively assessed and compared with the original plots from Figure <ref>.

§ RESULTS

§.§ Multi-Gaussian transformation

RBIGP, RBIGI, PPMT and FA were applied to three case studies with multivariate complexities. 150 iterations were used for RBIGP, RBIGI and PPMT. For FA, input parameters were $\sigma_0=0.1$ and $\sigma_1=1.1$ and this method was chained twice to achieve a standard multi-Gaussian distribution. As a result, all transforms produced independent multi-Gaussian factors for total and soluble copper grades in case A (Figure <ref>). RBIGI required 34 seconds, while PPMT and FA required 24 and 18 seconds, respectively.

Cross plots of multi-Gaussian transformation results for case A.

In case B, the non-linear relationship between iron and magnesium oxide was transformed into a multi-Gaussian distribution by all three methods (Figure <ref>). However, because this case has 9,990 sample points, it took slightly more time to execute: 50 seconds for RBIGI, 48 seconds for PPMT and 51 seconds for FA.

Cross plots of multi-Gaussian transformation results for case B.

Case C was more challenging to transform and Figure <ref> shows the transformed factors from all three methods. Declustering weights were not used in this case because the sum of the weights was not equal to the number of points. Furthermore, as case C comprises 33,021 titanium and zirconium samples, it takes more time to transform. For example, RBIGI and PPMT required 2 minutes 46 seconds and 2 minutes 42 seconds, respectively. On the other hand, it took FA 9 minutes 11 seconds to transform the heteroscedasticity between variables into independent factors with a multi-Gaussian distribution. In all three cases, RBIGP completed the transformation in less than one second. The significant difference in time between RBIGP and RBIGI can be explained by the different programming languages and simpler marginal Gaussianisation.

Cross plots of multi-Gaussian transformation results for case C.

Visual inspection of kernel density estimates can assist in assessing the quality of results. However, a better way to compare the results is to check for multivariate normality. Table <ref> shows the results from MVN tests applied to the transformed factors. For cases A and C, all three transforms produced multi-Gaussian distributions with perfect p-values according to the Henze-Zirkler and Energy tests. PPMT, however, could not produce multi-Gaussian factors and failed both MVN tests (i.e., with p-values of less than 0.05). RBIGI very slightly outperformed the other transforms based on multi-Gaussianity, particularly in cases A and B. PPMT failed to produce an MVN distribution in case B and showed slightly worse results in case A. On the other hand, all three methods produced similar results for case C. It is also important to note that the FA factors do not have a unit covariance matrix and have a standard deviation of 1.09 in all three cases.

Multivariate normality test results for three case studies

MVN test     Case     RBIGP     RBIGI     PPMT     FA

3*Henze-Zirkler’s     A     0.0932 (p = 1)     0.0205 (p = 1)     0.3023 (p = 1)     0.0237 (p = 1)

B     0.0621 (p = 1)     0.0226 (p = 1)     1.2135 (p $<$ 0.05)     0.0270 (p = 1)

C     0.3951 (p = 1)     0.3362 (p = 1)     0.3319 (p = 1)     0.3198 (p = 1)

3*Energy     A     0.0969 (p = 1)     0.0394 (p = 1)     0.2382 (p = 1)     0.0860 (p = 1)

B     0.0845 (p = 1)     0.0337 (p = 1)     1.2746 (p $<$ 0.05)     0.0740 (p = 1)

C     0.2835 (p = 1)     0.2214 (p = 1)     0.2307 (p = 1)     0.2377 (p = 1)

§.§ Spatial decorrelation

One of the limitations of the MGT approaches is that they can only assume spatial decorrelation. MGT methods always guarantee the decorrelation of variables at lag zero but not at further lags. Figure <ref> (left column) shows that the omni-directional cross-variograms of the MGT factors deviate slightly from zero in all three case studies. A better spatial decorrelation was achieved by chaining RBIGP, RBIGI, PPMT and FA with MAF. A 75 m lag for case A and a 15 m lag for cases B and C were selected for the MAF transformation, and the resulting matrices are shown in Table <ref>. As a result, the decorrelated factors for cases B and C are more spatially independent, whereas case A did not show much improvement. Nevertheless, the resulting factors are significantly more decorrelated than the normal score variograms.

Computed MAF matrices for multi-Gaussian factors in the three case studies

Case     Lag     RBIGP     RBIGI     PPMT     FA

A     75     $\begin{bmatrix}-0.90&-0.45\\0.45&-0.90\\\end{bmatrix}$     $\begin{bmatrix}-0.86&0.53\\-0.53&-0.86\\\end{bmatrix}$     $\begin{bmatrix}0.24&-0.98\\0.98&0.24\\\end{bmatrix}$     $\begin{bmatrix}0.23&-0.98\\0.98&0.23\\\end{bmatrix}$

B     15     $\begin{bmatrix}-0.20&-0.99\\0.99&-0.20\\\end{bmatrix}$     $\begin{bmatrix}-0.34&-0.95\\0.95&-0.34\\\end{bmatrix}$     $\begin{bmatrix}-0.79&0.64\\-0.64&-0.79\\\end{bmatrix}$     $\begin{bmatrix}-0.78&0.65\\-0.65&-0.78\\\end{bmatrix}$

C     15     $\begin{bmatrix}0.47&-0.89\\0.89&0.47\\\end{bmatrix}$     $\begin{bmatrix}-0.84&0.55\\-0.55&-0.84\\\end{bmatrix}$     $\begin{bmatrix}-0.77&-0.66\\0.66&-0.77\\\end{bmatrix}$     $\begin{bmatrix}-0.78&-0.65\\0.65&-0.78\\\end{bmatrix}$

Another way to assess spatial decorrelation is by checking the relative deviation from diagonality $\tau(h)$ and spatial diagonalisation efficiency $\kappa(h)$ introduced by [Tercan, 1999]. Figure <ref> (middle and right columns) shows that chaining MGT with MAF improves spatial decorrelation. Even though case A does not appear to be much better, there is a clear decorrelation at lower lags up to 100 m after applying MAF. It should also be noted that MGT methods alone are not sufficient to ensure decorrelation in case B. This is evident by $\tau(h)$ and $\kappa(h)$ not being close to zero and one, respectively. Such poor results can be explained by the normal score cross-variograms being very close to the factor cross-variograms. The average $\tau(h)$ and $\kappa(h)$ results are shown in Table <ref>.

Cross-variograms, relative deviation from diagonality $\tau(h)$ and spatial diagonalisation efficiency $\kappa(h)$ for case A (top row), case B (middle row) and case C (bottom row) before and after applying MAF.

Average relative deviation from diagonality $\tau(h)$ and spatial diagonalisation efficiency $\kappa(h)$

2*Method     2cCase A     2cCase B     2cCase C

$\overline{\tau}$     $\overline{\kappa}$     $\overline{\tau}$     $\overline{\kappa}$     $\overline{\tau}$     $\overline{\kappa}$

RBIGP     0.0550     0.9939     0.0825     0.8755     0.0285     0.9967

RBIGI     0.0606     0.9924     0.1106     0.8046     0.0287     0.9963

PPMT     0.0394     0.9960     0.1595     0.6562     0.0275     0.9960

FA     0.0397     0.9940     0.1662     0.4299     0.0279     0.9942

RBIGP+MAF     0.0416     0.9956     0.0364     0.9780     0.0105     0.9998

RBIGI+MAF     0.0424     0.9952     0.0338     0.9852     0.0107     0.9998

PPMT+MAF     0.0418     0.9956     0.0363     0.9805     0.0103     0.9998

FA+MAF     0.0426     0.9933     0.0364     0.9668     0.0103     0.9997

§.§ Geostatistical conditional simulation

Direct variograms were automatically fitted before conditionally simulating the decorrelated variables. Various automated and semi-automated variogram fitting algorithms can be used for this purpose [Emery, 2010, Desassis and Renard, 2013] and are available in various commercial software packages. There was no significant directional anisotropy for case A, and omni-directional variograms were modelled for the RBIG, PPMT and FA factors (Figure <ref>). The spatial variability is very similar among the three transforms, which means that any differences between generated realisations will be due mainly to the back-transformations.

\begin{equation}
\begin{aligned}
\gamma_{{CuT}_{RBIGP}}=0.342Sph(15m)+0.658Sph(82m) \\
\gamma_{{CuS}_{RBIGP}}=0.505Sph(15m)+0.495Sph(267m) \\
\gamma_{{CuT}_{RBIGI}}=0.342Sph(15m)+0.658Sph(82m) \\
\gamma_{{CuS}_{RBIGI}}=0.506Sph(15m)+0.494Sph(265m) \\
\gamma_{{CuT}_{PPMT}}=0.347Sph(15m)+0.653Sph(81m) \\
\gamma_{{CuS}_{PPMT}}=0.507Sph(15m)+0.493Sph(269m) \\
\gamma_{{CuT}_{FA}}=0.340Sph(15m)+0.660Sph(82m) \\
\gamma_{{CuS}_{FA}}=0.503Sph(15m)+0.497Sph(268m)
\end{aligned}
\label{eq6}
\end{equation}

Direct variograms of decorrelated factors of total (top row) and soluble (bottom row) copper grades for case A. Experimental variograms are indicated with points and fitted theoretical variograms with solid lines.

On the other hand, the data in case B showed significantly different variabilities in horizontal and vertical directions. Experimental variograms of decorrelated iron and magnesium oxide factors were automatically fitted as shown in Figure <ref>. However, for MAF to provide a decorrelation at all lags, spatial variability must be represented by a two-structured LMC [Desbarats and Dimitrakopoulos, 2000]. In cases A and B, the co-regionalisation models appear to be more complex and thus experimental variograms could not be fitted to two nested structures. Nevertheless, chaining MAF and MGT still demonstrates better results, even when two-structured LMC cannot be produced [Hosseini and Asghari, 2019, Erten and Deutsch, 2021].

\begin{equation}
\begin{aligned}
\gamma_{{Fe}_{RBIGP}}=0.380Nug+0.620Sph(20m,75m) \\
\gamma_{{MgO}_{RBIGP}}=0.145Nug+0.855Sph(33m,222m) \\
\gamma_{{Fe}_{RBIGI}}=0.370Nug+0.630Sph(20m,78m) \\
\gamma_{{MgO}_{RBIGI}}=0.140Nug+0.860Sph(33m,215m) \\
\gamma_{{Fe}_{PPMT}}=0.390Nug+0.610Sph(21m,81m) \\
\gamma_{{MgO}_{PPMT}}=0.146Nug+0.854Sph(34,213m) \\
\gamma_{{Fe}_{FA}}=0.395Nug+0.605Sph(21m,79m) \\
\gamma_{{MgO}_{FA}}=0.148Nug+0.852Sph(34m,220m)
\end{aligned}
\label{eq7}
\end{equation}

Direct variograms of decorrelated factors for iron (top row) and magnesium oxide (bottom row) for case B. Experimental variograms are indicated by points and fitted theoretical variograms by solid lines (blue: horizontal, and red: vertical).

Finally, omni-directional experimental variograms were calculated for case C. Figure <ref> shows the variograms for the RBIG, PPMT and FA factors. Unlike the other two cases, variograms in case C were fitted to two nested structures.

\begin{equation}
\begin{aligned}
\gamma_{{Ti}_{RBIGP}}=0.672Sph(15m)+0.328Sph(275m) \\
\gamma_{{Zr}_{RBIGP}}=0.530Sph(15m)+0.470Sph(275m) \\
\gamma_{{Ti}_{RBIGI}}=0.670Sph(15m)+0.330Sph(274m) \\
\gamma_{{Zr}_{RBIGI}}=0.531Sph(15m)+0.469Sph(274m) \\
\gamma_{{Ti}_{PPMT}}=0.665Sph(15m)+0.335Sph(275m) \\
\gamma_{{Zr}_{PPMT}}=0.534Sph(15m)+0.466Sph(275m) \\
\gamma_{{Ti}_{FA}}=0.668Sph(15m)+0.332Sph(275m) \\
\gamma_{{Zr}_{FA}}=0.534Sph(15m)+0.466Sph(275m)
\end{aligned}
\label{eq8}
\end{equation}

Direct variograms of decorrelated factors of titanium (top row) and zirconium (bottom row) for case C. Experimental variograms are indicated by points and the fitted theoretical variograms by solid lines.

It is evident that the factors produced by each method have almost identical spatial properties, which is true for all three cases. There were only minor differences in the variogram parameters generated by automated fitting. Using these variogram models, conditional turning bands simulation [Desassis and Renard, 1973, Emery and Lantuéjoul, 2006] was used to generate 50 realisations. In each case, the produced multi-Gaussian factors were modelled using identical neighbourhood parameters. Thus, most of the differences between realisations will be due to the variable transformation, which ensures a fair comparison of the MGT methods.

For case A, simulations were run on a 61 × 115 × 38 block model with grid dimension of 30m × 30m × 10m. A moving neighbourhood of 1000m × 1000m × 1000m with 8 octants and 8 points per octant was used for turning bands simulation. The block model in case B consists of 64 × 69 × 39 grids with 3m × 3m × 2m size, for which a 200m × 200m × 200m moving neighbourhood with 8 octants and 16 points per octant was used. Finally, realisations for case C were produced on a 65 × 36 × 54 block model with a 15m × 15m × 15m grid size. A 450m × 450m × 450m moving neighbourhood was used with 8 octants and 40 points per octant.

§.§ Analysis and validation

The simulated realisations were back-transformed to the original scale. To do so, MGT+MAF factors were first multiplied by the transposed MAF matrices (see Table <ref>), Then the RBIGP, RBIGI, PPMT and FA inverse transformations were applied to the corresponding multi-Gaussian realisations. In case A, RBIGP and PPMT required only 3 seconds to back-transform a single realisation and RBIGI required 22 seconds, whereas FA took 13 minutes and 29 seconds. A similar difference was observed in case B, where FA required 14 minutes and 30 seconds, RBIGI required 13 seconds and both RBIGP and PPMT required 2 seconds to back-transform one realisation. Case C is much larger than the other two, but RBIGP, RBIGI and PPMT needed only 1, 7 and 3 seconds, respectively. This is because their back-transformation does not require the original data, only Gaussian tables and rotation matrices at each iteration. However, as the original data are used in the FA back transformation, it took 34 minutes and 59 seconds to back-transform a single realisation for case C.

§.§.§ Reproduction of bivariate relationships

Figure <ref> (top) shows the cross plots of the back-transformed variables from a single realisation of case A. The bivariate relationships of the simulated data are similar to those of the original distribution. Although there are some artefacts above an inequality constraint, this is not surprising given the skewed distributions of the copper grades. In fact, even acceptance-rejection methods that reject and re-simulate values to be within the bounds of inequality constraints do so at the cost of other statistical properties. For example, the acceptance-rejection approach performs well when marginal distributions are moderately skewed [Madani and Abulkhair, 2020], but poorly reproduces other properties when dealing with very skewed distributions [Abulkhair and Madani, 2021]. Finally, reproduction of the Pearson and Spearman correlation coefficients suggests that PPMT performed slightly better than the others, while RBIGI underestimated the correlations (Figure <ref> bottom).

Top: cross plots of back-transformed simulated grades from one realisation compared to the original data in case A. Bottom: box plots of Pearson $\rho$ and Spearman $\rho_S$ correlation coefficients from 50 realisations of case A compared to the original values (solid blue lines).

In case B, the back-transformed results show a non-linear bivariate relationship, similar to the original data (Figure <ref>). Moreover, the reproduction of correlation coefficients is also almost identical.

Top: cross plots of back-transformed simulated grades from one realisation compared to the original data in case B. Bottom: box plots of Pearson $\rho$ and Spearman $\rho_S$ correlation coefficients from 50 realisations of case B compared to the original values (solid blue lines).

Finally, the cross plots in Figure <ref> demonstrate that all the presented methods reproduced the heteroscedastic relationship between titanium and zirconium in case C. It can also be observed in the reproduction of the Pearson and Spearman correlation coefficients, where all three methods show similar results. However, despite similar correlations, a visual inspection suggests that FA produces better results. RBIGI produced significant outliers in the top left and bottom right corners of its cross plots in cases A and the top right corner in case C (Figures <ref> and <ref>). On the other hand, RBIGP and PPMT produced minor artefacts in all three cases, which was not observed in the FA realisations.

Top: cross plots of back-transformed simulated grades from one realisation compared to the original data in case C. Bottom: box plots of Pearson $\rho$ and Spearman $\rho_S$ correlation coefficients from 50 realisations of case B compared to the original values (solid blue lines).

§.§.§ Histogram and variogram validation

Histogram reproduction is another important part of geostatistical validation. Histograms of back-transformed variables of a single realization were well reproduced in all three case studies (see Figures <ref>, <ref> and <ref>). However, it is more appropriate to check all realisations, which can be accomplished by CDF or Q-Q plots. Figure <ref> shows the reproduction of the CDFs of the back-transformed realisations and the RMSE values between percentiles of the realisations and the original distributions. In case A, the skewness of the distributions makes it challenging to assess and compare, but RMSE results suggest that RBIGP and FA performed slightly better than the others. Similarly, RBIGP has the least RMSE for both variables in case C, but other methods performed better in case B. Nevertheless, the difference in RMSE is insignificant, and the overall results are similar.

Histogram validation for case A (a), case B (b) and case C (c). The CDFs of the simulated grades are shown as solid black lines and the solid red lines are the original functions.

Reproduction of spatial variability was checked on the normal scores (i.e., first normal score transformation before iterations of RBIGP, RBIGI, PPMT and FA). The reason for this is to avoid significant deviations due to the skewness of the original variables. In addition, experimental normal score variograms were used as theoretical variograms of normal scores were not modelled. Finally, RMSE between variogram values for simulated and experimental variograms were calculated to provide a quantitative comparison of result.

Figure <ref> shows the variogram validation for case A. As expected, direct variograms are well reproduced, and cross-variograms of the simulated realisations are also similar to the original points. As was observed in the histogram validation, all three methods show almost identical results from a visual standpoint. Nevertheless, it is difficult to make a comparison based on RMSE results. RBIGP and FA have smaller RMSE for total copper, but RBIGI shows better reproduction for soluble copper and cross-variability. At the same time, RBIGI also has a higher RMSE for total copper.

Variogram validation for case A with direct variograms in the top and bottom rows and cross-variograms in the middle row. Solid lines represent normal score variograms of simulated realisations and experimental variograms are shown as points.

In case B, omni-horizontal and vertical variograms were used during the modelling due to anisotropy between those directions. Figure <ref> shows the variogram validation in both the vertical and horizontal directions. Visually, the results are identical for direct and cross-variograms in both directions. However, based on RMSE, RBIGP performed the best, followed closely by other methods. There is also an underestimation of the magnesium oxide variogram for all three methods.

Variogram validation for case B with direct variograms in the top and bottom rows and cross-variograms in the middle. Solid lines represent normal score variograms of simulated realisations and experimental variograms are shown as points (blue: horizontal, red: vertical).

Finally, variogram reproduction is visually indistinguishable in case C (Figure <ref>). RBIGP has a smaller RMSE for direct variograms, but the difference with other methods is insignificant. Overall, geostatistical validation metrics suggest that the performance of all four methods is similar, especially in the reproduction of marginal distributions and spatial variability. Moreover, without calculating RMSE, it is impossible to tell the difference between RBIGP, RBIGI, PPMT and FA results.

Variogram validation for case C with direct variograms in the top and bottom rows and cross-variogram in the middle row. Solid lines represent normal score variograms of simulated realisations and experimental variograms are shown as points.

§ DISCUSSION AND CONCLUSIONS

This paper compares four MGT methods, PPMT, FA and two types of RBIG (i.e., RBIGP with PCA rotations and RBIGI with ICA rotations), using three case studies with complex multivariate relationships. All three cases were obtained from undisclosed datasets, and each represents a particular bivariate complexity with moderate-to-high correlation coefficients. Case A is a bivariate dataset with an inequality constraint between total and soluble copper grades, with 6,074 data points. It can also be called a fractional constraint since the soluble copper grade is a fraction of the total grade. Case B consists of 9,990 iron and magnesium oxide samples with a non-linear relationship between them. Finally, case C has 33,021 data points, with clear heteroscedasticity between titanium and zirconium grades.

The most apparent difference between MGT approaches in this study is the execution time for forward and back transformation. Although RBIGI, PPMT and FA showed comparable times for forward transformation in cases A and B, FA took almost thrice as much time as PPMT and RBIGI in case C. Furthermore, the difference in execution time is significantly higher during the back transformation. RBIGP, RBIGI and PPMT do not need original data for the inverse transform because they save all the rotation matrices and Gaussian tables at each iteration during the fitting of transforms. Since FA requires the original data, the number of samples of the original hard data significantly affects its execution time. For example, FA required around 14 minutes to back-transform a single realization in cases A and B, while RBIGP, RBIGI and PPMT required only seconds. Moreover, because case C is significantly larger than the other two datasets, the time required by FA to back-transform one realisation was about 35 minutes. Due to a simpler and faster marginal Gaussianisation and more optimised matrix operations in MATLAB, RBIGP was significantly faster than other methods in forward transformation and similar to PPMT in inverse transformation. For example, RBIGP needed less than one second to transform the original data in all three cases.

Nevertheless, FA shows better reproduction of bivariate relationships, which was also observed by other studies reported in the geostatistical literature. PPMT has minor artefacts generated during back-transformation. While RBIGI appears to have good reproduction, it can produce significant outliers from the distribution but still within the bounds of marginals. The reason for these outliers is the convergence issue of ICA rotation, which can be unstable in some iterations. On the other hand, RBIGP does not appear to have this problem, producing results of similar quality to PPMT. Despite the minor differences, all three case studies showed that MGT methods effectively model datasets with multivariate complexities. Furthermore, all three methods show identical results in histogram and variogram validation. Even the RMSE between realisations and the original data shows no significant differences in the RBIGP, RBIGI, PPMT and FA results. Interestingly, the multivariate normality of the multi-Gaussian factors did not significantly impact the results. For example, PPMT failed MVN tests in case B but still produced realisations with a similar quality to those of other methods. However, multivariate normality plays a much bigger role in high-dimensional cases, particularly when correlations between variables are weak. Nevertheless, spatial decorrelation is much more important in geostatistical applications, and PPMT had similar spatial diagonalisation results to those of the other three methods.

Overall, the four MGT approaches analysed in this paper produce similar results in terms of multivariate normality, spatial decorrelation, reproduction of bivariate relationships, histogram and variogram validation. The main difference is the time it takes them to transform the original data and back-transform geostatistical realisations. It can be concluded that FA is not suitable for tasks requiring faster transformation, particularly when working on large datasets. For example, the recently developed rapid resource model updating methods require Gaussian transformation and must be in near real-time [Benndorf, 2020]. Geostatistical decorrelation methods, such as MAF and FA, were applied to rapidly update multiple cross-correlated variables [Kumar et al, 2020, Prior et al, 2021]. However, we believe that RBIG and PPMT are more suitable for such tasks, especially when chained with spatial decorrelation. In future work, we will apply RBIG and PPMT in the rapid updating of multiple variables within iron oxide copper-gold deposits.

Author contributions

Conceptualization: SA, PAD and CX; Methodology: SA; Investigation: SA; Writing - original draft preparation: SA; Writing - review and editing: SA, PAD and CX; Supervision: PAD and CX; Funding acquisition: PAD and CX. All authors read and approved the final manuscript.

The research reported here was supported by the Australian Research Council Industrial Transformation Training Centre for Integrated Operations for Complex Resources (ARC ITTC IOCR - project number IC190100017) and funded by universities, industry and the Australian Government. The first author also acknowledges the International Association for Mathematical Geosciences for providing a Travel Grant to attend the IAMG 2022 conference and present this study. Finally, we acknowledge Geovariances for providing a demonstration version of Isatis.neo Mining Edition software.

Code availability

Python package “pymgt” from ARC TC IOCR was used to perform RBIG and PPMT transformations (available at <https://github.com/exepulveda/pymgt>). FA was performed using the “gmGeostats” CRAN package (available at <https://cran.r-project.org/web/packages/gmGeostats>). For RBIGP, the original MATLAB code was used (available at <https://github.com/IPL-UV/rbig_matlab>). MAF, variogram fitting and geostatistical modelling were performed using the Isatis.neo software.

Conflicts of Interest

The authors declare no conflict of interest.

[Abildin et al, 2019]
Abildin Y, Madani N, Topal E (2019) A Hybrid Approach for Joint Simulation of
Geometallurgical Variables with Inequality Constraint. Minerals 9(1):24.

[Abulkhair and Madani, 2021]
Abulkhair S, Madani N (2021) Assessing heterotopic searching strategy in
hierarchical cosimulation for modeling the variables with inequality
constraints. Comptes Rendus Géoscience 353(1):115–134.

[Almeida and Journel, 1994]
Almeida AS, Journel AG (1994) Joint simulation of multiple variables with a
Markov-type coregionalization model. Mathematical Geology 26(5):565–588.

[Barnett et al, 2014]
Barnett RM, Manchuk JG, Deutsch CV (2014) Projection Pursuit Multivariate
Transform. Mathematical Geosciences 46:337–359.

[Barnett et al, 2016]
Barnett RM, Manchuk JG, Deutsch CV (2016) The Projection-Pursuit Multivariate
Transform for Improved Continuous Variable Modeling. SPE Journal
21(06):2010–2026. 10.2118/184388-PA

[Bassani et al, 2018]
Bassani MAA, Coimbra Leite Costa JF, Deutsch CV (2018) Multivariate
geostatistical simulation with sum and fraction constraints. Applied Earth
Science 127(3):83–93. 10.1080/25726838.2018.1468145

[Benndorf, 2020]
Benndorf J (2020) Closed Loop Management in Mineral Resource Extraction.
Springer, Cham

[Davis and Greenes, 1983]
Davis BM, Greenes KA (1983) Estimation using spatially distributed
multivariate data: An example with coal quality. Journal of the
International Association for Mathematical Geology 15:287–300.

[de Figueiredo et al, 2021]
de Figueiredo LP, Schmitz T, Lunelli R, et al (2021) Direct Multivariate
Simulation - A stepwise conditional transformation for multivariate
geostatistical simulation. Computers & Geosciences 147:104,659.

[Desassis and Renard, 1973]
Desassis N, Renard D (1973) The intrinsic random functions and their
applications. Advances in Applied Probability 5(3):439–468.

[Desassis and Renard, 2013]
Desassis N, Renard D (2013) Automatic Variogram Modeling by Iterative Least
Squares: Univariate and Multivariate Cases. Mathematical Geosciences
45:453–470. 10.1007/s11004-012-9434-1

[Desbarats and Dimitrakopoulos, 2000]
Desbarats A, Dimitrakopoulos R (2000) Geostatistical Simulation of
Regionalized Pore-Size Distributions Using Min/Max Autocorrelation Factors.
Mathematical Geology 32:919–942. 10.1023/A:1007570402430

[Deutsch and Journel, 1992]
Deutsch CV, Journel AG (1992) GSLIB: Geostatistical Software Library and User's
Guide. Oxford University Press, New York

[Ebner and Henze, 2020]
Ebner B, Henze N (2020) Tests for multivariate normality—a critical review
with emphasis on weighted ${L}^2$-statistics. TEST 29:845–892.

[Egozcue et al, 2003]
Egozcue JJ, Pawlowsky-Glahn V, Mateu-Figueras G, et al (2003) Reproduction of
secondary data in projection pursuit transformation. Mathematical Geology
35:279–300. 10.1023/A:1023818214614

[Emery, 2008]
Emery X (2008) A turning bands program for conditional co-simulation of
cross-correlated Gaussian random fields. Computers & Geosciences
34(12):1850–1862. 10.1016/j.cageo.2007.10.007

[Emery, 2010]
Emery X (2010) Iterative algorithms for fitting a linear model of
coregionalization. Computers & Geosciences 36(9):1150–1160.

[Emery and Lantuéjoul, 2006]
Emery X, Lantuéjoul C (2006) TBSIM: A computer program for conditional
simulation of three-dimensional Gaussian random fields via the turning bands
method. Computers & Geosciences 32(10):1615–1628.

[Erten and Deutsch, 2021]
Erten O, Deutsch CV (2021) Assessment of variogram reproduction in the
simulation of decorrelated factors. Stochastic Environmental Research and
Risk Assessment 35:2583–2604. 10.1007/s00477-021-02005-0

[Friedman, 1987]
Friedman JH (1987) Exploratory Projection Pursuit. Journal of the American
Statistical Association 82(397):249–266.

[Henze and Zirkler, 1990]
Henze N, Zirkler B (1990) A class of invariant consistent tests for
multivariate normality. Communications in Statistics - Theory and Methods
19(10):3595–3617. 10.1080/03610929008830400

[Hosseini and Asghari, 2015]
Hosseini SA, Asghari O (2015) Simulation of geometallurgical variables through
stepwise conditional transformation in Sungun copper deposit, Iran. Arabian
Journal of Geosciences 8:3821–3831. 10.1007/s12517-014-1452-5

[Hosseini and Asghari, 2019]
Hosseini SA, Asghari O (2019) Multivariate Geostatistical Simulation on
Block-Support in the Presence of Complex Multivariate Relationships: Iron Ore
Deposit Case Study. Natural Resources Research 28:125–144.

[Hyvarinen, 1999]
Hyvarinen A (1999) Fast and robust fixed-point algorithms for independent
component analysis. IEEE transactions on Neural Networks 10(3):626–634.

[Joenssen and Vogel, 2014]
Joenssen DW, Vogel J (2014) A power study of goodness-of-fit tests for
multivariate normality implemented in R. Journal of Statistical Computation
and Simulation 84(5):1055–1078. 10.1080/00949655.2012.739620

[Journel, 1999]
Journel AG (1999) Markov models for cross-covariances. Mathematical Geology
31(8):955–964. 10.1023/A:1007553013388

[Journel and Huijbregts, 1978]
Journel AG, Huijbregts CJ (1978) Mining geostatistics. Academic Press, London

[Korkmaz et al, 2014]
Korkmaz S, Göksülük D, Zararsiz G (2014) MVN: An R Package for Assessing
Multivariate Normality. R JOURNAL 6(2):151–162

[Kumar et al, 2020]
Kumar A, Dimitrakopoulos R, Maulen M (2020) Adaptive self-learning mechanisms
for updating short-term production decisions in an industrial mining
complex. Journal of Intelligent Manufacturing 31:1795–1811.

[Laparra et al, 2011]
Laparra V, Camps-Valls G, Malo J (2011) Iterative Gaussianization: From ICA to
Random Rotations. IEEE Transactions on Neural Networks 22(4):537–549.

[Leuangthong and Deutsch, 2003]
Leuangthong O, Deutsch CV (2003) Stepwise Conditional Transformation for
Simulation of Multiple Variables. Mathematical Geology 35(2):155–173.

[Madani and Abulkhair, 2020]
Madani N, Abulkhair S (2020) A hierarchical cosimulation algorithm integrated
with an acceptance–rejection method for the geostatistical modeling of
variables with inequality constraints. Stochastic Environmental Research and
Risk Assessment 34:1559–1589. 10.1007/s00477-020-01838-5

[Manchuk et al, 2017]
Manchuk JG, Barnett RM, Deutsch CV (2017) Reproduction of secondary data in
projection pursuit transformation. Stochastic Environmental Research and
Risk Assessment 31:2585–2605. 10.1007/s00477-016-1363-y

[Mardia, 1970]
Mardia KV (1970) Measures of multivariate skewness and kurtosis with
applications. Biometrika 57(3):519–530. 10.1093/biomet/57.3.519

[Mueller and Ferreira, 2012]
Mueller UA, Ferreira J (2012) The U-WEDGE Transformation Method for
Multivariate Geostatistical Simulation. Mathematical Geosciences
44:427–448. 10.1007/s11004-012-9384-7

[Mueller et al, 2020]
Mueller UA, Tolosana-Delgado R, Grunsky EC, et al (2020) Biplots for
compositional data derived from generalized joint diagonalization methods.
Applied Computing and Geosciences 8:100,044.

[Pawlowsky-Glahn et al, 2015]
Pawlowsky-Glahn V, Egozcue JJ, Tolosana-Delgado R (2015) Modeling and analysis
of compositional data. John Wiley & Sons

[Pedregosa et al, 2011]
Pedregosa F, Varoquaux G, Gramfort A, et al (2011) Scikit-learn: Machine
Learning in Python. Journal of Machine Learning Research 12:2825–2830

[Prior et al, 2021]
Prior A, Tolosana-Delgado R, van den Boogaart KG, et al (2021) Resource
Model Updating For Compositional Geometallurgical Variables. Mathematical
Geosciences 53:945–968. 10.1007/s11004-020-09874-1

[Rossi and Deutsch, 2014]
Rossi ME, Deutsch CV (2014) Mineral resource estimation. Springer, Dordrecht

[Royston, 1983]
Royston J (1983) Some Techniques for Assessing Multivarate Normality Based on
the Shapiro-Wilk W. Journal of the Royal Statistical Society: Series C
(Applied Statistics) 32(2):121–133. 10.2307/2347291

[Székely and Rizzo, 2013]
Székely GJ, Rizzo ML (2013) Energy statistics: A class of statistics based on
distances. Journal of statistical planning and inference 143(8):1249–1272.

[Talebi et al, 2019]
Talebi H, Mueller U, Tolosana-Delgado R, et al (2019) Geostatistical
Simulation of Geochemical Compositions in the Presence of Multiple Geological
Units: Application to Mineral Resource Evaluation. Mathematical Geosciences
51:129–153. 10.1007/s11004-018-9763-9

[Tercan and Sohrabian, 2013]
Tercan A, Sohrabian B (2013) Multivariate geostatistical simulation of coal
quality data by independent components. International Journal of Coal
Geology 112:53–66. 10.1016/j.coal.2012.10.007

[Tercan, 1999]
Tercan AE (1999) Importance of Orthogonalization Algorithm in Modeling
Conditional Distributions by Orthogonal Transformed Indicator Methods.
Mathematical Geology 31:155–173. 10.1023/A:1007557701073

[Tolosana-Delgado and Mueller, 2021]
Tolosana-Delgado R, Mueller U (2021) Geostatistics for Compositional Data with
R. Springer, Cham

[Tolosana-Delgado et al, 2019]
Tolosana-Delgado R, Mueller U, van den Boogaart KG (2019) Geostatistics for
compositional data: an overview. Mathematical Geosciences 51(4):485–526.

[van den Boogaart et al, 2017]
van den Boogaart KG, Mueller U, Tolosana-Delgado R (2017) An Affine
Equivariant Multivariate Normal Score Transform for Compositional Data.
Mathematical Geosciences 49:231–251. 10.1007/s11004-016-9645-y

[Verly, 1993]
Verly G (1993) Sequential Gaussian cosimulation: a simulation method
integrating several types of information. In: Soares A (ed) Geostatistics
Tróia ’92. Quantitative Geology and Geostatistics, vol 5. Springer,
Dordrecht, p 543–554, 10.1007/978-94-011-1739-5_42

[Wackernagel, 2003]
Wackernagel H (2003) Multivariate Geostatistics: An Introduction with
Applications. Springer, Berlin