# Developing a simulation tool to investigate a novel trapped two-state Bose-
Einstein condensate Ramsey interferometer driven by dipole oscillations and
gravitational sag

Anushka Thenuwara  Centre for Optical Sciences, Faculty of Science,
Engineering and Technology, Swinburne University of Technology, John Street,
Hawthorn, Victoria, Australia 3122  School of Physics and Astronomy, Monash
University Clayton Campus, Wellington Rd, Clayton, Victoria, Australia 3800
<EMAIL_ADDRESS>these authors contributed equally to this work
Andrei Sidorov  Centre for Optical Sciences, Faculty of Science, Engineering
and Technology, Swinburne University of Technology, John Street, Hawthorn,
Victoria, Australia 3122  these authors contributed equally to this work

###### Abstract

We propose and explore the feasibility of a novel Ramsey interferometer
created by a trapped two-state Bose-Einstein condensate (BEC) driven by dipole
oscillations and gravitational sag. The BEC is formed in a pure cigar shaped
compressed magnetic trap (CMT) via a dilute atom cloud of ${}^{87}Rb$ atoms in
state $\ket{F=2,m_{F}=+2}$ $(\ket{+2})$ of the $5^{2}S_{\frac{1}{2}}$ ground
state. Here, Rmasey interferometry is performed with states
$\ket{F=2,m_{F}=+1}$ $(\ket{+1})$ and $\ket{+2}$. The proposed interferometer
utilises the response of atoms to the harmonic oscillator trapping potential
and the gravitational sag due to the variation in the $m_{F}$ state. Briefly,
the state $\ket{+1}$ experiences a shallower radial trap with a larger
gravitational sag; whereas, state $\ket{+2}$ experiences a tighter radial trap
with a gravitational sag which is half of state $\ket{+1}$. Due to this, a
superposition between the states $\ket{+1}$ and $\ket{+2}$ experiences
multipath propagation resulting in an interference pattern. This may be
utilised to measure local gravitational fields and measure inter-sate
scattering lengths. Here, a theoretical framework is reported which is
developed via the two-level system in combination with the Gross-Pitaevskii
equation (GPE). Further, the development of a simulation tool via GPELabs in
MATLAB that explores the prosed interferometer is reported along with key
insights and findings.

## Introduction

The $5^{2}S_{\frac{1}{2}}$ ground state of ${}^{87}Rb$ atoms have three
magnetically trappable states: $\ket{m_{F}=+2}(\ket{+2})$ and
$\ket{m_{F}=+1}(\ket{+1})$ in the $F=2$ hyperfine manifold and
$\ket{m_{F}=-1}(\ket{-1})$ in the $F=1$ hyperfine manifold. Two of those
states ($\ket{+1}$ and $\ket{-1}$) have the same magnetic moment, can be
coupled by a two-photon MW-RF transition and have been extensively studied [1,
2, 3]. Moreover, several examples of Ramsey interferometry in the
$\ket{-1}-\ket{+1}$ system are presented in [4, 5, 6, 7]. The work in [4]
shows the spatially dependent relative phase evolution of an elongated two-
component Bose-Einstein condensate between the states $\ket{-1}$ and
$\ket{+1}$. In it, the relative phase evolution is probed via Ramsey
interferometry and goes on to show that the mean-field formalism provides a
good description of the decay of the Ramsey signal. It also concludes that the
loss of Ramsey contrast is due to the inhomogeneity of the collective relative
phase across the cloud rather than to decoherence or phase diffusion. Ramsey
interferometry with spin-echo techniques in the $\ket{-1}-\ket{+1}$ system is
performed in [5] and reports coherence times of several seconds. The cause of
this extended coherence is due to mean-field-driven collective oscillations of
the two components leading to periodic de-phasing and re-phasing of the
condensate wavefunctions with a slow decay of the interference fringe
visibility. The work in [6] describes the Wigner calculations developed with
the experimental data in [5] that goes beyond classical GPE modelling.
Finally, the work in [7] utilises Ramsey interferometry in the
$\ket{-1}-\ket{+1}$ system to perform precise scattering length measurements.
It reports the inter-state scattering length of
$a_{\ket{-1}-\ket{+1}}=98.006(16)a_{0}$ with an uncertainty of
$1.6\text{\times}{10}^{-4}$ and the intra-state scattering length of
$a_{\ket{+1}-\ket{+1}}=95.44(7)a_{0}$.

This paper focusses on exploring the capabilities of performing Ramsey
interferometry in the $\ket{+1}-\ket{+2}$ system. The states $\ket{+2}$ and
$\ket{+1}$ can be coupled by a single-photon magnetic dipole transition.They
have different magnetic moments resulting in different radial trapping
frequencies $\omega_{r}$ based on

$\omega_{r}=\sqrt{\frac{\mu_{B}m_{F}g_{F}}{m}\frac{\partial^{2}B}{\partial
r^{2}}},$ (1)

where $\mu_{B}$ is the Bohr magneton, $m_{F}$ is the magnetic quantum number,
$g_{F}$ is the Landé factor of the $F$ hyperfine state, $m$ is the mass of the
atom and $\frac{\partial^{2}B}{\partial r^{2}}$ is the curvature of the
magnetic field in the radial direction.

Further, BECs are affected by the gravity and the vertical position is shifted
by the gravitational sag $z_{g}$

$z_{g}=\frac{g}{\omega_{r}^{2}},$ (2)

where $g$ is the gravitational acceleration.

We propose a Ramsey interferometer by combining the effects of the variation
in $m_{F}$ values on $\omega_{r}$ and $z_{g}$ in the $\ket{+1}-\ket{+2}$
system. Firstly, prepare a BEC in state $\ket{+2}$ and create a superposition
by transfer a certain fraction of atoms to state $\ket{+1}$. This will be the
first pulse of the Ramsey sequence ($\hat{U}_{\tau}$ in Equation 9) where the
experimental results in [8] show the ability of the $\ket{+1}-\ket{+2}$ BEC
system to create a variety of superpositions. Once the superpostion is
created, the state dependent radial trap frequency as per Equation 1 and state
dependent gravitational sag as per Equation 2 cause the $\ket{+1}-\ket{+2}$
superposition to undergo dipole oscillations. Here, the fraction of the
superposition in state $\ket{+1}$ is on the slope of its trapping potential
which slides down, travels towards and goes beyond its trap minimum along the
$z$-axis, and oscillate back to the original position. The dipole oscillation
time of $\ket{+1}$ is $T=\frac{2\pi}{\omega_{1r}}$ which is the free evolution
time $T$ of the Ramsey sequence ($\hat{U}_{T}$ in Equation 9). This free
evolution time can be varied by adjusting the radial trap frequency via
changing the curvature of the external magnetic field. Once the superposition
undergoes one dipole oscillation, the recombination pulse $\hat{U}_{\tau}$ is
applied with varying phase $\phi$ (Equation 9). This completes the Ramsey
interferometric sequence, where the resulting system wavefunction follows
Equation 11. Based on the above, the proposed Ramsey sequence is depicted in
Figure 1.

Figure 1: Ramsey interferometry of the trapped two-state $\ket{+1}-\ket{+2}$
BEC system subjected gravitational sagging and dipole oscillations. Here, the
graph next to the BEC is a schematic representation of the trapping potentials
for states $\ket{+1}$ (blue) and $\ket{+2}$ (red).

Primarily, we explore the plausibility of creating the aforementioned Ramsey
interferometer governed by the interesting characteristics of gravitational
sagging and dipole oscillations. This plausibility is theoretically explored
in [9] where analytical expressions are derived to fit experimental results.
Here, we develop and perform numerical simulations to get an insight on how
the BECs behave during experimental runs. It should be noted that these
simulations follow the two-level model and only consider the trapped
$\ket{+1}-\ket{+2}$ states of the ${}^{87}Rb$ BEC in the hyperfine $F=2$
manifold of the $5^{2}S_{\frac{1}{2}}$ ground state. This paper is structured
in three sections, which outlines the theoretical framework, the development
of the simulation tool and the results from the simulations. The section on
the theoretical framework presents three areas. Firstly, the Rabi model for a
two-level atom which describes the creation of a superposition of states via
an electromagnetic (EM) field through the matrix formalism. Secondly, it
presents the description of the Ramsey interferometric sequence for the two-
state model through the unitary evolution operator. Thirdly, it presents the
theoretical description of the Gross-Pitaevskii equation (GPE) and the general
background behind numerical simulations of GPE. The next section presents the
development of the simulation tool which introduces the GPELabs MATLAB
toolbox. Then goes on to breakdown simulated equation for the proposed
interferometer along with capabilities of the developed simulation tool.
Finally, the last section presents the results of the simulations.

## Theoretical framework

### Matrix description of Rabi oscillations in a two-level system

As first presented in [10] and carefully explained in [11], a collection of
two-level atoms undergoes population oscillations between the two quantum
states when coupled to an external oscillating electromagnetic (EM) field.
These oscillations are known as Rabi oscillations and [10] provides an elegant
description on how this effect is utilised for precision measurements of
magnetic moments of atoms. As a brief description of the underlying process,
two-level atoms create a coherent superposition of the two states evolving
with time when exposed to a quasi-resonant EM field. This evolution is
dependent on several major parameters such as the dipole moment of the
transition, the amplitude of the EM wave and the difference (detuning) between
the transition frequency and the EM frequency.

Figure 2: Two-level atom interacting with an oscillating EM field where the
energy separation between the two states is $\hbar\omega_{0}$.

Consider the two-level system of $\ket{1}$ and $\ket{2}$ quantum states shown
in Figure 2 where the energy separation between the two states is
$\hbar\omega_{0}$. We use this definition of the quantum states because it
will be later related to the two magnetically trappable states
$\ket{m_{F}=+1}$ and $\ket{m_{F}=+2}$ in the $F=2$ hyperfine manifold of
${}^{87}Rb$ atoms. Consider this system being coupled to an external
oscillating EM field with a frequency $\omega_{EM}=(\omega_{0}-\Delta)$ where
$\Delta$ is the detuning. We will consider the magnetic dipole interaction
between the atom and an oscillating magnetic field, so the expression for the
magnetic component of the EM field is
$B_{0}\cos(\omega_{EM}t-\phi)=\frac{B_{0}}{2}\left(e^{i(\omega_{EM}t-\phi)}+e^{-i(\omega_{EM}t-\phi)}\right)$,
where $B_{0}$ and $\phi$ are respectively the amplitude and phase. For the
state vector,

$\Psi_{sys}=\sum_{i=1}^{2}C_{i}\ket{i}=\begin{bmatrix}C_{1}\\\ C_{2}\\\
\end{bmatrix},$ (3)

we write the time-dependent Schrödinger equation in the rotating wave
approximation

$i\hbar\frac{\partial\ket{\Psi_{sys}}}{\partial
t}=\hat{H}\ket{\Psi_{sys}}=\left(\hat{H}_{A}+\hat{H}_{I}\right)\ket{\Psi_{sys}},$
(4)

where $\hat{H}_{A}$ and $\hat{H}_{I}$ are the bare-atom and interaction
Hamiltonians

$\hat{H}_{A}=\frac{\hbar}{2}\begin{bmatrix}\omega_{0}&0\\\
0&-\omega_{0}\end{bmatrix},\hskip
28.45274pt{}\hat{H}_{I}=\frac{\hbar}{2}\begin{bmatrix}\Delta&\Omega_{R}e^{-i\phi}\\\
\Omega_{R}e^{i\phi}&-\Delta\end{bmatrix},$ (5)

where $\omega_{0}=\omega_{2}-\omega_{1}$, $\Delta$ is the detuning of the EM
field from resonance and
$\Omega_{R}=\frac{\bra{1}\hat{\mu}\ket{2}B_{0}}{\hbar}$ is the Rabi frequency
describing the magnetic dipole coupling of the $\ket{1}$ and $\ket{2}$ states.

We can use a matrix evolution formalism where the solution of the differential
equation 4 can be replaced by a time-evolving Hamiltonian applied to
$\ket{\Psi_{sys}(t)}$ at $t=0$ to obtain the time evolution of the system.
Instead of using the bare state basis of Equation 3, we can move to the
interaction picture and use the basis of the interaction Hamiltonian
$\hat{H}_{I}$. We derive the eigenvalues and eigenvectors of $\hat{H}_{I}$ and
apply them to the principles of unitary time evolution to obtain
$\hat{H}_{I}(t)$. This keeps the basis states independent of time and obtains
$\ket{\Psi_{sys}(t)}$ by applying the matrix operator $\hat{H}_{I}(t)$ to the
initial conditions of the system. The eigenvalues of $\hat{H}_{I}$ are
$\lambda_{\pm}=\pm\frac{\hbar\Omega_{G}}{2}$ and the eigenvectors for these
values are

$\ket{V_{+}}=\begin{bmatrix}e^{\frac{-i\phi}{2}}\cos(\frac{\theta}{2})\\\
e^{\frac{i\phi}{2}}\cos(\frac{\theta}{2})\end{bmatrix},\hskip
28.45274pt{}\ket{V_{-}}=\begin{bmatrix}-e^{\frac{-i\phi}{2}}\sin(\frac{\theta}{2})\\\
e^{\frac{i\phi}{2}}\cos(\frac{\theta}{2})\end{bmatrix},$ (6)

where $\tan(\theta)=\frac{\Omega_{R}}{\Delta}$. The application of unitary
time evolution
$\hat{H}_{I}(t)=e^{\frac{-i\lambda_{+}t}{\hbar}}\ket{V_{+}}\bra{V_{+}}+e^{\frac{-i\lambda_{-}t}{\hbar}}\ket{V_{-}}\bra{V_{-}}$
leads to the unitary time evolution operator

$\hat{U}_{t}=\begin{bmatrix}\cos\left(\frac{\Omega_{G}t}{2}\right)-\frac{i\Delta}{\Omega_{G}}\sin\left(\frac{\Omega_{G}t}{2}\right)&-ie^{-i\phi}\frac{\Omega_{R}}{\Omega_{G}}\sin\left(\frac{\Omega_{G}t}{2}\right)\\\
-ie^{i\phi}\frac{\Omega_{R}}{\Omega_{G}}\sin\left(\frac{\Omega_{G}t}{2}\right)&\cos\left(\frac{\Omega_{G}t}{2}\right)+\frac{i\Delta}{\Omega_{G}}\sin\left(\frac{\Omega_{G}t}{2}\right)\end{bmatrix}.$
(7)

The evolution of the system wavefuction $\ket{\Psi_{sys}}$ when the system
starts in state $\ket{2}$ (e.g., the $\ket{m_{F}=+2}$ state of the ${}^{87}Rb$
atom) can be obtained by applying the $\hat{U}_{t}$ operator to the initial
state vector $\begin{bmatrix}0&1\end{bmatrix}^{T}$. This leads to

$\ket{\Psi_{sys}(t)}=\hat{U}_{t}\begin{bmatrix}0\\\
1\end{bmatrix}=\begin{bmatrix}-ie^{-i\phi}\frac{\Omega_{R}}{\Omega_{G}}\sin\left(\frac{\Omega_{G}t}{2}\right)\\\
\cos\left(\frac{\Omega_{G}t}{2}\right)+\frac{i\Delta}{\Omega_{G}}\sin\left(\frac{\Omega_{G}t}{2}\right)\end{bmatrix}.$
(8)

### Ramsey method of separated oscillatory fields

Following the work of I.I. Rabi, N.F. Ramsey [12] significantly improved the
Rabi method by using two oscillatory fields with a short pulse length $\tau$
separated by a long free evolution time $T$ to study molecular resonances and
demonstrated linewidths that are $0.6$ times narrower compared to that of the
single oscillatory field method. This work of N.F. Ramsey won him the Nobel
prize in physics in $1989$ and is commonly referred to as Ramsey
interferometry in the physics community. The Ramsey interferometric method
provides the basis of the exquisite time standards using $Cs$ fountain clocks.
As examples, the NIM5 $Cs$ fountain clock [13] based in China has an
uncertainity of $1.6\text{\times}{10}^{-15}$ and the NIST-F1 $Cs$ fountain
clock [14] based in the USA has an uncertainty of
$0.97\text{\times}{10}^{-15}$. These uncerainities translate to a
$1\text{\,}\mathrm{s}$ error in time keeping for every $20$ to $30$ million
years. Further, Ramsey interferometry is utilised to obtain extremely
sensitive measurements of local gravity as reported in [15] and a Ramsey-type
method with a spin-echo pulse (i.e. a $\pi$-pulse during the free evolution
time $T$) is utilised in [16] to obtain precise measurements of the Newtonian
gravitational constant
$G=$6.671\,91(99)\text{\times}{10}^{-11}\text{\,}\mathrm{m}^{3}\mathrm{k}\mathrm{g}^{-1}\mathrm{s}^{-2}$$.

A typical Ramsey sequence interrogates a single two-level atom with a single-
frequency quasi-resonant electromagnetic field. Briefly, Ramsey interferometry
synchronizes two independent oscillators (an electromagnetic field and a
coherent superposition of two quantum states of an atom) and compares their
evolution in time. The atom is initially in state $\ket{2}$. The first
$\frac{\pi}{2}$-pulse of duration $\tau$ (Rabi method) prepares a coherent
superposition of the two states with equal ($50$-$50$) populations which is
synchronized with the EM radiation. Then the atom-EM coupling is turned off
and the two oscillators evolve independently for the evolution time $T$. The
second (recombination) $\frac{\pi}{2}$-pulse interrogates the superposition
state and extracts the phase difference (the relative phase) between the two
quantum states. This relative phase is contained in the population difference
of the states [4]. Further, the dynamical behavior of the two-level system
wavefunction and its behaviour between the two states can be mathematically
represented via the Bloch vector formalism [17].

As derived in Equation 7, the unitary time operator $\hat{U}$ can be used to
determine the system wavefunction $\Psi_{sys}(t)$ at the end of each step of
the Ramsey interferometric sequence. One of the unitary time evolution
operators relevant to the sequence is the operation for the combination and
the recombination pulses of the superposition. A $$50$:$50$$ superposition in
a two-level system is created when $\Omega_{G}t=\frac{\pi}{2}$. Applying this
condition to $\hat{U}_{t}$ in Equation 7 leads to the following expression for
the unitary evolution operator $\hat{U}_{\tau}$ for the $\frac{\pi}{2}$-pulse.
The other important case is the free evolution operator $\hat{U}_{T}$ which is
when $\Omega_{R}=0$ for $\hat{U}_{t}$ in Equation 7 leading to

$\hat{U}_{\tau}=\frac{1}{\sqrt{2}}\begin{bmatrix}1&-ie^{-i\phi}\\\
-ie^{i\phi}&1\end{bmatrix},\hskip
28.45274pt{}\hat{U}_{T}=\begin{bmatrix}e^{-\frac{i\Delta T}{2}}&0\\\
0&e^{\frac{i\Delta T}{2}}\end{bmatrix}.$ (9)

The evolution operators above are sequentially applied to the initial state of
$\Psi_{sys}(0)$ as

$\Psi_{sys}(t)=\hat{U}_{\tau}.\hat{U}_{T}.\hat{U}_{\tau}.\Psi_{sys}(0),$ (10)

to obtain an expression for $\Psi_{sys}(t)$ at the end of the Ramsey sequence.
This leads to

$\Psi(\tau+T+\tau)=\begin{bmatrix}-ie^{\frac{-i\phi}{2}}\cos(\frac{\Delta
T-\phi}{2})\\\ ie^{\frac{i\phi}{2}}\sin(\frac{\Delta
T-\phi}{2})\end{bmatrix},$ (11)

where $T$ is the free evolution time, $\Delta$ and $\phi$ are respectively the
detuning of the EM field from resonance and the phase of the second
$\frac{\pi}{2}$-pulse.

The intermediate expressions for $\Psi_{sys}(t)$ at different stages of the
Ramsey sequence are shown in Figure 3. The equation for $\Psi_{sys}(t)$ at the
end of the Ramsey sequence depends on the free evolution time $T$, the
detuning of the EM field from resonance $\Delta$ and the phase $\phi$ of the
second $\frac{\pi}{2}$-pulse. These are the three domains (the time domain,
the frequency domain and the phase domain) in which Ramsey interferometry can
be preformed. It should be noted that Ramsey interferometry in the frequency
domain may be riddled by the shift of the detuning over the course of the
experiment and frequency instability (e.g., magnetic noise) in the equipment.
However, the scheme proposed in [18] is an elegant method to utilise Ramsey
interferometry in the frequency domain. Further, Ramsey interferometry is also
performed in the atom number $N$ domain [5]. Here, we utilises Ramsey
interferometry in the phase domain where $\phi$ of the second
$\frac{\pi}{2}$-pulse is manipulated to extract the interferometric data
between $\ket{1}$ and $\ket{2}$ for a fixed $T$. The fixed $T$ is one of the
major advantages of operating in the phase domain as it avoids the shot-to-
shot variation of atom numbers due to trap losses for varying hold times.

Figure 3: Evolution of the state vector $\Psi_{sys}(t)$ of a two-level atom in
a typical Ramsey interferometric sequence.

The interesting aspect in this work is the state dependent spatial separation
of the two states due to gravitational sag based on Equation 2. This causes
the superposition of states to spatially separate and recombine during $T$.
The overlap of the two components before the recombination pulse defines the
visibility of the Ramsey fringe as denoted by Equation 12 [5]. At the end of
the Ramsey sequence we introduce a variable $P_{z}(T,\Delta,\phi)$ where,

$P_{z}(T,\Delta,\phi)=\frac{2}{N}Im\left[e^{i(\phi+\Delta
T)}\int\Psi_{(\tau+T)}^{\ket{+2}*}\Psi_{(\tau+T)}^{\ket{+1}}\right]d^{3}r,$
(12)

which contains all dependencies and can be measured via an absorption image.

Next we shall investigate the theoretical framework that describes
interactions of the proposed two-level Ramsey interferometric sequence.

### The Gross-Pitaevskii equation and its numerical simulation

The dynamics of the $\ket{+1}-\ket{+2}$ BEC system during the Ramsey
interferometric sequence is analysed via the Gross-Pitaevskii equation (GPE).
This foundational theoretical model was introduced by E. P. Gross and L. P.
Pitaevskii in 1961 [19] for analysing the dynamical evolution of the zero
temperature condensate wavefunction. The GPE is an extension of the standard
Schrödinger equation which incorporates interactions between atoms in a dilute
cloud through a mean-field interaction potential $U_{int}(r)$. The GPE model
has well-defined many experimental data whilst reducing the calculation power
necessary as shown in Figure 4 in [4]. The GPE is also used complementarily
with experimental results to extract enhanced measurements as shown in [7, 5].

The time-dependent GPE takes the form of Equation 13 [20], which describes
both the temporal and spatial dynamics of the BEC system.

$i\hbar\frac{\partial\Psi(r,t)}{\partial
t}=-\frac{\hbar^{2}}{2m}\triangledown^{2}\Psi(r,t)+V_{HO}(r)\Psi(r,t)+U_{0}|\Psi(r,t)|^{2}\Psi(r,t),$
(13)

where $\Psi(r,t)$ is the time-dependent wavefunction,
$\triangledown^{2}=\frac{\partial^{2}}{\partial r^{2}}$ is Laplace operator,
$V_{HO}$ is the external harmonic oscillator trapping potential and $U_{0}$ is
the effective interaction parameter.

The first term in the equation encapsulates the kinetic energy of the system
and the second term captures the potential energy via the external harmonic
potential defined by
$V_{HO}=\frac{m}{2}\left(2\omega_{r}^{2}r^{2}+\omega_{y}^{2}y^{2}\right)$ for
a cigar-shaped trap, where $\omega_{r,y}$ are the radial and axial trap
frequencies, respectively. The last term in the equation describes the energy
variation due to interactions with the special parameter $U_{0}$ know as the
interaction constant defined as $U_{0}=\frac{4\pi\hbar^{2}a_{s}}{m}$, where
$a_{s}$ is the s-wave scattering length between the interacting atoms.

The GPE in the form of Equation 13 describes the temporal and spatial
evolution of a single BEC with internal collisional interactions. For the two-
component $\ket{+1}-\ket{+2}$ BEC system, the GPE requires to include intra-
state (within the states) and inter-state (between states) interactions
resulting in a set of equations that describes a two-component BEC. The intra-
state interactions are described by $U_{ii}=\frac{4\pi\hbar^{2}a_{ii}}{m}$,
where $a_{ii}$ is the intra-state scattering length for elastic collisions
between atoms of the same state $i$. The inter-state interactions are
described by $U_{ij}=\frac{4\pi\hbar^{2}a_{ij}}{m}$ , where $a_{ij}$ is the
inter-state scattering length for elastic collisions between atoms of adjacent
states $ij$. This lead to coupled equations for a two-component interacting
BEC system of the form

$\begin{split}i\hbar\frac{\partial\Psi_{1}(r,t)}{\partial
t}&=-\frac{\hbar^{2}}{2m}\triangledown^{2}\Psi_{1}(r,t)+V_{HO1}(r)\Psi_{1}(r,t)+U_{11}|\Psi_{1}(r,t)|^{2}\Psi_{1}(r,t)+U_{12}|\Psi_{2}(r,t)|^{2}\Psi_{1}(r,t),\\\
i\hbar\frac{\partial\Psi_{2}(r,t)}{\partial
t}&=-\frac{\hbar^{2}}{2m}\triangledown^{2}\Psi_{2}(r,t)+V_{HO2}(r)\Psi_{2}(r,t)+U_{22}|\Psi_{2}(r,t)|^{2}\Psi_{2}(r,t)+U_{12}|\Psi_{1}(r,t)|^{2}\Psi_{2}(r,t),\end{split}$
(14)

which can be solved via numerical methods. Further, each equation has
parameters unique to either state allowing to incorporate state-dependent
external trapping potentials, where effects of gravitational sag are
incorporated into the analysis. Also, additional aspects of interactions with
an external EM field and atom losses from either state will be added leading
to an expanded analysis.

Numerical simulations are of great importance to obtain an insight into the
$\ket{+1}-\ket{+2}$ BEC dynamics. As a precursor, the coupled Gross-Pitaevskii
equations in Equation 14 are useful to analyse the impact of the inter-state
interactions on the suggested Ramsey sequence as well as to probe optimal
conditions to maximise the Ramsey contrast at the end of one dipole
oscillation. Though the in-house developed xSPDE Matlab simulation toolbox is
capable of simulating these coupled partial differential GPEs [21], the
theoretical framework presented in [22, 23] is used as the cornerstone for the
numerical simulations in this project. The main advantage of [22, 23] is the
readily available GPELabs MATLAB tool box [24, 25] which does not require
developing the simulation environment compared with xSPDE. Further, xSPDE is
tailored towards studying stochastic effects which is not considered here.

The GPELabs MATLAB tool box utilises several finite difference and spectral
methods to calculate the ground-state and the dynamics of wavefunctions for
multi-component Bose-Einstein condensate systems. To calculate the ground-
state, GPELabs offers four schemes: backward Euler finite difference scheme
(BEFD), Crank-Nicolson finite difference scheme (CNFD), backward Euler
spectral discretization scheme (BESP) and Crank-Nicolson spectral
discretization scheme (CNSP). From these the first two are finite difference
schemes and the second two are spectral schemes. To calculate the dynamics,
GPELabs has six schemes with two additional pseudo spectral schemes to the
above which are time splitting pseudo spectral scheme (TSSP) and relaxation
pseudo spectral scheme (ReSP).

An overall summary of the numerical schemes relating to solving BEC systems is
presented in [26]. It reports a vast variety of schemes spanning the areas of
simulating BECs at zero temperature, BECs at finite temperatures with a
thermal fraction, as well as systems of fermions and boson-fermion mixtures.
Here we are focused on a pure BEC system without any thermal fraction. A
detailed study of the TSSP method is presented in [27] and an extensive
descriptions of both BEFD and TSSP are presented in [28] which reserves the
BEFD for ground-state calculations with proven energy diminishing condition
under normalized gradient flow and reserves the TSSP for dynamical simulations
of BECs. The simulated ground-states via BEFD are confirmed by comparing with
the experimental results for the two component BEC system in the immiscible
regime in [29] and goes on to simulate the dynamics of the same system via
TSSP showing the separation and recombination of the two components due to
miscibility. In [30] the numerical results from BEFD and TSSP are compared
with CNFD and forward Euler finite difference (FEFD) methods. It shows that
BEFD and TSSP are much better in preserving the energy diminishing property
under continuous normalized gradient flow and reports the importance of
picking the appropriate spatial mesh and time-step size to obtain correct
results via CNFD and FEFD schemes. The work presented in [31, 32, 33] outlines
all the calculation methods discussed above in several avenues; however, the
optimal method for ground-state and dynamic calculations of coupled GPEs is
considered. Regarding ground-state calculations, [32, 33] compares the BEFD
and BESP methods and concludes BESP is better as it requires fewer grid points
which saves computational memory and time especially for the 3D scenario.
Regarding the dynamic scenario, the framework for the ReSP scheme was first
introuduced in [34], which reports a vast reduction in computational time for
the relaxation scheme compared to that of CNFD and highlights its ability to
conserve the density and energy compared to the TSSP method. For these
reasons, the BESP scheme is selected to calculate the ground-state and the
ReSP scheme is selected to calculate the dynamics of the two component coupled
GPEs. Further, scheme dependent variations of the simulated outputs are
considered out of the current scope.

## The development of the simulation tool

MATLAB was used as the primary software tool to simulate the two coupled
Gross-Pitaevskii equations for the two component BEC system due to the
dedicated GEPLab toolbox developed by Xavier Antoine and Romani Duboscq [24,
25]. This toolbox has the capability to solve the GPE to find the wavefunction
of the ground-state and the time-dependent behaviour (including Rabi coupling
and dipole oscillations) for a single or multi component BEC system in 1D, 2D
or 3D.

The GPE is calculated in the dimensionless form in GPELabs which vastly
expands its applicability for BEC systems with non-linear Schrödinger
equations. The dimensionless form of the GPE as per [24] is

$\left\\{i\frac{\partial\Psi(r,t)}{\partial
t}=\left(\frac{-1}{2}\Delta+V+\beta|\Psi(r,t)|^{2}-\Omega
L_{z}\right)\Psi(r,t)\right\\},$ (15)

where mathematical expressions within curly braces ($\\{Expression\\}$) or
variables with a bar ($\bar{X}$) denotes the dimensionless form.

An important facet of the simulation tool is the conversion between the
dimensionless values and the real values. It is important to encode the
experimental scenario into the simulation and then to decode simulated results
into real values. The experimental trap frequencies provide the footing for
the unit of time for the dimensionless framework where
$\omega_{m}=\min(\omega_{x,y,z})$. The unit length is then defined via the
harmonic oscillator characteristic length
$a_{0}=\sqrt{\frac{\hbar}{m\omega_{m}}}$, [24]. These allow to define the
dimensionless unit of time by $\bar{t}=\omega_{m}t$ and the dimensionless unit
of length $\bar{r}=\frac{r}{a_{0}}$. Based on these the 3D wavefunction has
the conversion
$\bar{\Psi}_{3D}(\bar{r},\bar{t})=a_{0}^{\frac{3}{2}}\Psi_{3D}(r,t)$. It
should be noted that the unit of $\Psi(r,t)$ is dependent upon the
dimensionality such that $[\Psi_{nD}]=[L]^{\frac{-n}{2}}$. Therefore, care
must be taken when converting the GPE between dimensions. All frequencies are
relative to $\omega_{m}$ where $\bar{\Omega}=\frac{\Omega}{\omega_{m}}$.
Finally, the dimensionless energy takes the form as
$\bar{E}=\frac{E}{\hbar\omega_{m}}$. These relations are of great importance
to correlate the experimental scenario with simulated results. Table 1
provides a summary of the conversion factors that convert simulated results to
real values where the dimensionless values are denoted by $\bar{X}$.

Real value | Converting form simulations
---|---
$t$ |  $\frac{\bar{t}}{\omega_{m}}$
$r$ |  $a_{0}\bar{r}$
$\Psi_{3D}(r,t)$ |  $a_{0}^{\frac{-3}{2}}\bar{\Psi}_{3D}(\bar{r},\bar{t})$
$\Omega$ |  $\bar{\Omega}\omega_{m}$
$E$ |  $\bar{E}\hbar\omega_{m}$
Table 1: Conversions of parameters from simulation to real values where the
dimensionless values from the simulation are denoted by $\bar{X}$

Furthermore, the simulated system of wavefunctions is normalised to unity
within GPELabs where $\int\Psi_{sys}^{*}\Psi_{sys}=1$. However, most of the
general literature normalise a system of wavefunctions to the total population
$N$ where $\int\Psi_{sys}^{*}\Psi_{sys}=N$. This should be carefully
considered when adapting terms into GPELabs as some terms acquire a factor of
$N$ to the standard expression. Due to the complicated nature of converting a
3D experimental to a 1D approximation, the simulation was developed for the 3D
scenarios.

The terms relevant to this work in Equation 15 are the Laplace term or the
dispersion retaliation ($\bar{\Delta}$), the linear term ($\bar{V}$) and the
non-linear term ($\\{\beta|\Psi(r,t)|^{2}\\}$). The angular momentum term for
the BEC ($\\{\Omega L_{z}\\}$) is not considered, leading an equation with the
following structure

$\left\\{\begin{bmatrix}i\dot{\Psi}_{1}\\\
i\dot{\Psi}_{2}\end{bmatrix}=\left(\frac{-1}{2}\begin{bmatrix}\Delta_{11}&0\\\
0&\Delta_{22}\end{bmatrix}+\begin{bmatrix}L_{11}&L_{12}\\\
L_{21}&L_{22}\end{bmatrix}+\begin{bmatrix}NL_{11}&NL_{12}\\\
NL_{21}&NL_{22}\end{bmatrix}\right)\begin{bmatrix}\Psi_{1}\\\
\Psi_{2}\end{bmatrix}\right\\},$ (16)

where $\bar{\Psi}_{n}$ is the dimensionless form of the wavefunction for state
$\ket{+n}$, $\bar{\Delta}_{ii}$ are dispersion relations, $\bar{L}_{ij}$ are
the linear terms and $\bar{NL_{ij}}$ are non-linear terms. The advantage of
this interpretation is the ease with which the linear terms and the non-linear
terms can be compartmentalised to simulate the BEC system.

The conversion of dimensionless parameters is well described in [22] which
shows the adaptation of equations in 1D, 2D and 3D. The work in [22] shows a
dependency of selecting a characteristic length scale with the interaction
strength of the simulated scenario. In short, $a_{0}$ is a valid selection if
the BEC is in the weak interaction regime where $4\pi|a_{s}|N\ll a_{0}$ or in
the moderate interaction regime where $4\pi|a_{s}|N\approx a_{0}$, where
$a_{s}$ is the relevant scattering length and $N$ is the atom number. If the
simulation is in the strong interaction regime where $4\pi|a_{s}|N\gg a_{0}$,
a better characteristic length $x_{s}$ is suggested instead of $a_{0}$ where
$x_{s}=(4\pi|a_{s}|Na^{4}_{0})^{\frac{1}{5}}$ and goes on to explain that the
selection of the appropriate characteristic length depicts the visibility of
phenomena and what can be resolved by the discretization of specified spatial
and temporal grids. This does not appear in the GPELabs documentation;
however, GPELabs provides full control over this discretization over space and
time which leads to needing several trails to set-up the appropriate grid to
observe the dynamics of the system. Based on the conversion factors in Table
1, the simulated scenario can be extracted out to the real world.

The outputs of GPELabs are in both numerical and graphical forms for all
components. The numerical outputs are the square of the wavefunction at the
origin $|\bar{\Psi}_{3D}(0,\bar{t})|^{2}$, the directional rms-radii, the
energy and the chemical potential of each component. The graphical outputs are
the three-dimensional distribution of the square of the wavefunction
$|\bar{\Psi}_{3D}(\bar{r},\bar{t})|^{2}$ and the phase for each component. The
frequency of acquiring this information can be set at desired intervals of
iterations.

As suggested above, the separation of linear terms and non-linear terms can be
done when all simulated phenomena are combined and adapted to the framework of
GPELabs. The linear terms $\bar{L}$ take the form

$\left\\{\begin{bmatrix}L_{11}&L_{12}\\\
L_{21}&L_{22}\end{bmatrix}=\begin{bmatrix}V^{HO}_{11}+V^{EM}_{11}&V^{EM}_{12}\\\
V^{EM}_{21}&V^{HO}_{22}+V^{EM}_{22}\end{bmatrix}\right\\},$ (17)

which accounts for the state dependent harmonic oscillator (HO) trapping
potential ($\bar{V}^{HO}_{i,i}$) and the interaction with the external EM
field ($\bar{V}^{EM}_{ij}$) as both of these are multiplied linearly by the
wavefunction. The HO potential takes up the diagonal terms where the
interaction with the EM field takes up both diagonal and off-diagonal terms.
When the EM field is applied, the off-diagonal elements couple the two states
of the wavefunction into a superposition of states. The non-linear terms
$\bar{NL}$ take the form

$\left\\{\begin{bmatrix}NL_{11}&NL_{12}\\\
NL_{21}&NL_{22}\end{bmatrix}=\begin{bmatrix}V^{INT}_{11}+V^{Loss}_{11}&0\\\
0&V^{INT}_{22}+V^{Loss}_{22}\end{bmatrix}\right\\},$ (18)

where the off-diagonal elements are zero and the diagonal terms account for
both interactions ($\bar{V}^{INT}_{i,i}$: inter- and intra-state) and loss
($\bar{V}^{Loss}_{i,i}$: two- and three-body) terms as both of these are
multiplied by a higher-order of the wavefunction. The inclusion of the loss
terms expand the capabilities of the simulation tool, but are not considered
in the primary results. Within the considered theoretical framework, the non-
linear interactions are impacted by the populations of either state. Further,
higher-order couplings that lead to effects such as the Josephson effect [35]
are not considered.

### Adapting linear terms to the GPELabs framework

The linear terms of the simulated phenomena are the external HO trapping
potential with gravitational sagging and the interactions with the EM fields.
The external trapping potential has the form
$V^{HO}=\frac{m}{2}\left(\omega^{2}_{x}x^{2}+\omega^{2}_{y}y^{2}+\omega^{2}_{z}z^{2}\right)$.
Further, when the effects of gravitational sagging based on Equation 2 are
added,
$V^{HO}=\frac{m}{2}\left(\omega^{2}_{x}x^{2}+\omega^{2}_{y}y^{2}+\omega^{2}_{z}(z-\Delta
z_{sag})^{2}\right)$. The dimensionless form of $\bar{V}^{HO}$ is

$\bar{V}^{HO}=\left\\{\begin{bmatrix}\frac{1}{2}\left(\gamma^{2}_{1x}x^{2}+\gamma^{2}_{1y}y^{2}+\gamma^{2}_{1z}(z-\Delta
z_{1sag})^{2}\right)&0\\\
0&\frac{1}{2}\left(\gamma^{2}_{2x}x^{2}+\gamma^{2}_{2y}y^{2}+\gamma^{2}_{2z}(z-\Delta
z_{2sag})^{2}\right)\end{bmatrix}\right\\},$

(19)

where $\bar{\gamma}_{x,y,z}=\frac{\omega_{x,y,z}}{\omega_{m}}$ are
dimensionless directional trap frequencies and
$\Delta\bar{z}_{sag}=\frac{\Delta z_{sag}}{a_{0}}$ is the dimensionless
gravitational sag of each component

A simulation for a harmonic oscillator potential with
$\omega_{r}=2\pi\times$400\text{\,}\mathrm{Hz}$$ for state $\ket{+1}$ in blue
and state $\ket{+2}$ in red is shown in the Figure 4. The horizontal axis of
the graph is the dimensionless vertical distance $z/a_{0}$ from the atom chip.
The $y$-axis of the graph is energy in dimensionless radial frequency units
($\bar{\Omega}$). A key feature displayed in the figure is the spatial
separation of the trapping potentials due to gravitational sagging. Moreover,
an important insight is the energy gap $\bar{\Omega}_{12}$ felt by atoms in
state $\ket{+1}$ at its starting point which is the potential minimum $\Delta
z_{2sag}$ of state $\ket{+2}$. This energy gap should be accounted for when
applying the EM field as otherwise the atoms in state $\ket{+1}$ are untrapped
hindering the population transfer in creating a superposition.

Figure 4: External vertical trapping potential with gravitational sagging
experienced by $\ket{+1}$ (blue) and $\ket{+2}$ (red) with trap frequency
$\omega_{r}=2\pi\times$400\text{\,}\mathrm{Hz}$$ for the $\ket{+1}$ state.

The other linear term that influences this project is the interaction term
with the external EM field. This follows Equation 5 leading to
$V^{EM}=\hat{H}_{I}$. The dimensionless form for this can also be derived by
converting the standard GPE into the dimensionless form

$\bar{V}^{EM}=\frac{1}{2\omega_{m}}\begin{bmatrix}0&\Omega_{R}e^{-i\phi}\\\
\Omega_{R}e^{i\phi}&-2\Delta\end{bmatrix}=\frac{1}{2}\left\\{\begin{bmatrix}0&\bar{\Omega}_{R}e^{-i\phi}\\\
\bar{\Omega}_{R}e^{i\phi}&-2\bar{\Delta}\end{bmatrix}\right\\},$ (20)

where $\bar{\Delta}=\frac{\Delta}{\omega_{m}}$ is the dimensionless detuning
of EM field from the resonance between $\ket{+1}$ and $\ket{+2}$ and
$\bar{\Omega}_{R}=\frac{\Omega_{R}}{\omega_{m}}$ is the dimensionless Rabi
frequency for RF coupling between the two states.

An interesting point is the lack of a conversion factor for the relative phase
$\phi$ of the EM field when converting to its dimensionless form. Therefore,
the value of $\phi$ in the simulation is the value for the phase in the real
experiment. There is a subtle connection between the two linear terms: the
trapping potential and the applied EM field. As per Figure 4 gravitational
sagging causes a separation of the two potential minima for either state of
the $\ket{+1}-\ket{+2}$ system. This leads to an energy gap between the two
states at $\Delta\bar{z}_{2sag}$ where the $\ket{+2}$ BEC is created. This gap
needs to be filled to bring the two states on-resonance prior to applying the
mixing pulse via the EM field. If the EM field is applied off-resonance, then
the transferred atoms from $\ket{+2}\rightarrow\ket{+1}$ are immediately lost
as there is no trap at $\Delta\bar{z}_{2sag}$ for $\ket{+1}$ atoms. The
detuning $\bar{\Delta}$ of the EM field is utilized to bring the two states to
on-resonance, where $\bar{\Omega}_{12}=2\bar{\Delta}$ at which point the
transfer of atoms from $\ket{+2}\rightarrow\ket{+1}$ is most efficient. The
detuning $\Delta$ required for the EM field of the real experiment can be
extrapolated out from $\bar{\Omega}_{12}$ in the simulation as
$\bar{\Delta}=\frac{\bar{\Omega}_{12}}{2}$ and
$\Delta=\frac{\bar{\Omega}_{12}}{2}\omega_{m}$, where $\bar{\Omega}_{12}$ is
readily available via the simulation. This may be of great assistance to
provide an estimated starting point for the experiment.

### Adapting non-linear terms to the GPELabs framework

There are two non-linear terms incorporated into the simulation tool which are
the interaction term $V^{INT}_{ii}$ and the loss term $V^{Loss}_{ii}$. The
general form of the interaction term follows $V^{INT}=U_{ij}|\Psi_{k}|^{2}$,
where $U_{ij}=\frac{4\pi\hbar^{2}a_{ij}}{m}$ and $i,j,k\in\\{1,2\\}$. The
dimensionless form of $\bar{V}^{INT}$ when $\int\Psi_{sys}^{*}\Psi_{sys}=1$ is

$\bar{V}^{INT}=\left\\{4\pi
N\begin{bmatrix}a_{11}|\Psi_{1}|^{2}+a_{12}|\Psi_{2}|^{2}&0\\\
0&a_{22}|\Psi_{2}|^{2}+a_{12}|\Psi_{1}|^{2}\end{bmatrix}\right\\},$ (21)

where $a_{ij}$ are the scattering lengths between the interacting two states,
$\bar{a}_{ij}=\frac{a_{ij}}{a_{0}}$ and $N$ is the total atom number of the
$\ket{+1}-\ket{+2}$ system.

Further, careful consideration should be given to $\bar{V}^{INT}$ when
converting to 1D and 2D simulations. As described on pages 324 - 339 of [20],
there is a dimensional dependence on the interaction constant $U_{ij}$. For
the 1D case, $\bar{U}^{1D}=\frac{2\hbar^{2}}{m}\frac{a_{s}}{a^{2}_{\perp}}$,
where $a_{\perp}=\sqrt{\frac{\hbar}{m\omega_{\perp}}}$. This is to be combined
with the dimensionality of the wavefunction when converting $U^{1D}$ into the
dimensionless form where $\bar{\Psi}^{1D}=a_{\perp}^{\frac{1}{2}}\Psi^{1D}$.
Further, for the 2D case,
$\bar{U}^{2D}=\frac{\sqrt{8\pi}\hbar^{2}}{m}\frac{a_{s}}{a_{z}}$, where
$a_{z}=\sqrt{\frac{\hbar}{m\omega_{z}}}$. Also, similar to the 1D case, the
dimensionality of the wavefunction should be accounted for when converting
$U^{2D}$ into the dimensionless form where
$\bar{\Psi}^{2D}=a_{\perp}\Psi^{2D}$.

Moreover, the additional aspect of loss terms $V^{Loss}_{ii}$ was coded into
the simulation, where the two-body and three-body loss terms take the general
form of $V^{Loss}=-i\hbar\Gamma$. By compacting all the loss terms together
and converting the 3D GPE to its dimensionless form along with the
normalisation $\int\Psi_{sys}^{*}\Psi_{sys}=1$ results in

$\bar{V}^{Loss}_{ii}=\begin{bmatrix}\frac{-iN(\xi_{11}|\Psi_{1}|^{2}+\xi_{12}|\Psi_{2}|^{2})}{2\omega_{m}a^{3}_{0}}&0\\\
0&\frac{-iN\xi_{12}|\Psi_{1}|^{2}}{2\omega_{m}a^{3}_{0}}+\frac{-iN\xi_{222}|\Psi_{2}|^{4}}{2\omega_{m}a^{6}_{0}},\end{bmatrix},$
(22)

where $N$ is the total atom number, $\xi$ are two- and three-body loss terms,
$\omega_{m}$ is the minimum trap frequency and $a_{0}$ is the harmonic
oscillator characteristic length. The two-body loss term for $\ket{+2}$ with
$\xi_{12}=0$ for the $\ket{+1}-\ket{+2}$ interaction as this interaction can
only switch the states between the two atoms causing no loss. The two-body
loss rate for $\ket{+1}$ with loss coefficient
$\xi_{11}=$8.1(3)\text{\times}{10}^{-14}\text{\,}\mathrm{c}\mathrm{m}^{3}\mathrm{/}\mathrm{s}$$
[7] for the $\ket{+1}-\ket{+1}$ interactions. The three-body loss rate for
$\ket{+2}$ with the loss coefficient
$\xi_{222}=$1.8\text{\times}{10}^{-29}\text{\,}\mathrm{c}\mathrm{m}^{6}\mathrm{/}\mathrm{s}$$
[36] which is the only loss channel via three-body recombination of three
$\ket{+2}$ atoms.

## Simulations of Ramsey interferometry in the trapped $\ket{+1}$ \-
$\ket{+2}$ BEC system

The BEC in the experiment has a total atom number of about
$N=$3\text{\times}{10}^{4}\text{\,}$$ and the conditions were selected such
that estimated axial and radial trapping frequencies for state $\ket{+1}$ are
$\omega_{ax}\approx 2\pi\times$400\text{\,}\mathrm{Hz}$$ and
$\omega_{r}\approx 2\pi\times$12\text{\,}\mathrm{Hz}$$. Here, the notations
for the properties of states $\ket{+1}$ and $\ket{+2}$ are represented by
subscripts $X_{1},X_{2}$. The values for the scattering length were set as
$a_{11}=95.44a_{0}$ and $a_{22}=98.98a_{0}$ [37]. Initially, we set the
parameters $a_{12}$, two- and three-body loss terms to zero to obtain insights
into the dynamics of Ramsey interferometry.

For these conditions, first the ground-state for the $\ket{+2}$ BEC was
simulated. To create a $50$-$50$ superposition, the two trapping potentials
required an energy gap of $\Delta\bar{\Omega}_{12}=$64\text{\,}$$ to bring
into resonance. This relates to a detuning of $\Delta\approx
2\pi\times$2.4\text{\,}\mathrm{kHz}$$. Further, a high Rabi frequency of
$10\text{\,}\mathrm{kHz}$ was selected. This is to minimize the loss of atoms
via leaking to the $\ket{+1}$ potential during the Rabi oscillations. Based on
this, a $50$-$50$ superposition was created with a resonant Rabi pulse of
duration $\bar{t}_{R}=0.00189$ corresponding to an actual pulse duration of
$t_{R}=$25\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ for the chosen trap
frequencies. Figure 5 shows the two-component BEC once this stage was reached.
Thereafter, several separate scenarios were simulated with the first to verify
the Ramsey signal and others for the standard Ramsey sequence.

Figure 5: The spatially overlapping two-component BEC profiles after the first
$\frac{\pi}{2}$ Rabi pulse.

In the first scenario, a second recombination pulse with varying phase was
applied immediately after the superposition was created. Ramsey fringes
(Figure 6 b)) exhibit an interference contrast of $100\text{\,}\%$ and a full
inversion of the populations for the zero phase of the recombination pulse.
This is the expected result confirming the capability of the simulation to
perform Ramsey interferometry in the phase domain.

Figure 6: Ramsey interference of the $\ket{+1}-\ket{+2}$ BEC system when
$a_{12}=0$ with immediate second $\frac{\pi}{2}$-pulse. a) Population
variation for $\ket{+1}$ (blue) and $\ket{+2}$ (red) BECs and b) Ramsey signal
in the phase domain.

The second simulated scenario was to check on dipole oscillations and to
observe key physical phenomena that a $\ket{+1}-\ket{+2}$ BEC system may
exhibit during the free evolution time. The centre of mass oscillations are
shown in Figure 7 a) and as expected the $\ket{+2}$ BEC remained stationary at
its trap bottom while the $\ket{+1}$ BEC oscillated back and forth. The
oscillation period was $0.1882$ in dimensionless units which relates to a
dipole oscillation frequency of $\omega_{DP}\approx
2\pi\times$399.3\text{\,}\mathrm{Hz}$$. This is in very good agreement with
the $2\pi\times$400\text{\,}\mathrm{Hz}$$ value of the $\ket{+1}$ radial
frequency. Further, the expected idea is for the $\ket{+1}$ BEC to return to
the starting spatial point where it overlaps the $\ket{+2}$ BEC. At this point
the second $\frac{\pi}{2}$-pulse in the phase domain can be applied to obtain
the Ramsey signal. However, as shown in Figure 7 b), the overlap region of the
two components is reduced.

Figure 7: a) Centre of mass positions after the first $\frac{\pi}{2}$-pulse of
$\ket{+1}$(blue) and $-\ket{+2}$ (red) with $a_{12}=0$. b) Overlap of
condensate wavefunctions $|\Psi_{1}\rvert^{2}$ (blue)and
$\lvert\Psi_{2}\rvert^{2}$ (red) at the end of the first dipole oscillation in
the case of $a_{12}=0$.

This is due to the settling of either wavefunction to its new chemical
potential during the free evolution time causing dynamical variation of the
full-width at half maximum (FWHM). At the end of the first
$\frac{\pi}{2}$-pulse, the atom number in $\ket{+2}$ is reduced by one half.
This causes its chemical potential to reduce to $\approx$76\text{\,}\%$$ of
its original value and the theoretical TF radius in the radial direction
reduces to $\approx$87\text{\,}\%$$ of its original value. The atoms in
$\ket{+1}$ experience an expansion of their TF radius due to the relaxed
trapping potential. During the free evolution time atoms in $\ket{+1}$ acquire
a larger FWHM. Due to these two conditions, the overlap region reduces at the
maximum overlap point.

The reduction in the overlap region as shown in Figure 8 reduces the Ramsey
contrast as depicted in Equation 12 which prevents the full inversion of
populations. When the Ramsey sequence is done the $\ket{+1}-\ket{+2}$ BEC
system with $a_{12}=0$ after one dipole oscillation, the resulting Ramsey
contrast reduces to $\approx$50\text{\,}\%$$ of its full value as shown in
Figure 8. This is an important find as it shows a reduction in the Ramsey
contrast even when the inter-state scattering length $a_{12}=0$, where the
contrast reduction mechanism is via the reduced overlap between the two states
after one dipole oscillation of $\ket{+1}$.

Figure 8: Ramsey interference of the $\ket{+1}-\ket{+2}$ BEC system when
$a_{12}=0$ with second $\frac{\pi}{2}$-pulse applied after one dipole
oscillation. a) Population variation for $\ket{+1}$ (blue) and $\ket{+2}$
(red) BECs and b) Ramsey signal in the phase domain.

The next scenarios consider the impact of the inter-state scattering length on
the Ramsey signal. The two considered cases are when $a_{12}=98a_{0}$ and
$a_{12}=-50a_{0}$ where the focus is the positive and negative nature of
$a_{12}$. The initial half of the simulation was repeated for these new
configurations in both scenarios and the behaviour of the overlap region is
observed.

The main finding for $a_{12}=$98\text{\,}\mathrm{a}$_{0}$ is reflected in the
centre of mass motions of the two states during the dipole oscillations. As
shown in Figure 9, during the first separation the $\ket{+1}$ BEC exerts a
force that pushes the $\ket{+2}$ BEC away from its rest point. This force
originates from the non zero inter-state scattering length which makes the two
states susceptible to each other. Further, the two BECs repel each other away
when $\ket{+1}$ BEC returns back from one dipole oscillation which is clearly
shown at $0.18$ dimensionless time units in Figure 9. This leads to a further
reduction in the overlap region between the two states at the recombination
pulse. Due to this the population does not undergo full inversion at the
second $\frac{\pi}{2}$-pulse (Figure 10 a)) and the Ramsey signal decreases
(Figure 10 b)) to about $35\text{\,}\%$ of its full value.

Figure 9: Centre of mass positions after the first $\frac{\pi}{2}$-pulse of
$\ket{+1}$(blue) and $-\ket{+2}$ (red) with $a_{12}=98a_{0}$. Figure 10:
Ramsey interference of the $\ket{+1}-\ket{+2}$ BEC system when
$a_{12}=98a_{0}$ with second $\frac{\pi}{2}$-pulse applied after one dipole
oscillation. a) Population variation for $\ket{+1}$ (blue) and $\ket{+2}$
(red) BECs and b) Ramsey signal in the phase domain.

When $a_{12}=-50a_{0}$, the two-component system exhibits attractive
interactions which is clearly reflected in Figure 11. Here, the centre of mass
motion of $\ket{+1}$ pulls on $\ket{+2}$ during the initial stage of the
dipole oscillation and the effect is reversed when the overlap occurs at about
$0.18$ dimensionless time units. This has an impact on the overlap of the two
components at the end of one dipole oscillation which is reflected in the
Ramsey signal (Figure 12 b)) which is reduced to $\approx$45\text{\,}\%$$ of
its full value. This is better than the case of $a_{12}=98a_{0}$, but not as
good as the best possible contrast of $50\text{\,}\%$ for the zero inter-state
scattering length.

Figure 11: Centre of mass positions after the first $\frac{\pi}{2}$-pulse of
$\ket{+1}$(blue) and $-\ket{+2}$ (red) with $a_{12}=-50a_{0}$.

Figure 12: Ramsey interference of the $\ket{+1}-\ket{+2}$ BEC system when
$a_{12}=-50a_{0}$ with second $\frac{\pi}{2}$-pulse applied after one dipole
oscillation. a) Population variation for $\ket{+1}$ (blue) and $\ket{+2}$
(red) BECs and b) Ramsey signal in the phase domain.

Finally, Figure 13 compares the Ramsey signals between the various scenarios
considered above. The blue solid line with stars is when $a_{12}=0$, where the
recombination pulse is applied immediately after the first pulse in the Ramsey
sequence when both states do not move. This has a phase offset of
$4.7\text{\,}\mathrm{rad}$ with Ramsey signal at $100\text{\,}\%$. The red
dashed line with open circles is when $a_{12}=0$, where the recombination
pulse is applied after one dipole oscillation of $\ket{+1}$ in the Ramsey
sequence. This has a phase offset of $0.3\text{\,}\mathrm{rad}$ with Ramsey
signal at about $50\text{\,}\%$. The particular feature is that $\ket{+2}$
does not move during this sequence and only $\ket{+1}$ moves (Figure 7 a))
which contributes to the phase offset compared to the in-trap scenario. The
next case is the black dashed line when $a_{12}=98a_{0}$, where the
recombination pulse is applied after one dipole oscillation of $\ket{+1}$ in
the Ramsey sequence. This has a phase offset of $3.5\text{\,}\mathrm{rad}$
with the Ramsey signal at about $35\text{\,}\%$. Here, the $\ket{+2}$ also
moves during this sequence (Figure 9) which contributes to the phase offset of
about $3.2\text{\,}\mathrm{rad}$ compared to $a_{12}=0$ with one dipole
oscillation (red dashed line with open circles). Finally, the green solid line
with squares when $a_{12}=-50a_{0}$, where the recombination pulse is applied
after one dipole oscillation of $\ket{+1}$ in the Ramsey sequence. This has a
phase offset of $2.0\text{\,}\mathrm{rad}$ with Ramsey signal at about
$45\text{\,}\%$. Here, the $\ket{+2}$ moves the most out of all the cases
during the Ramsey sequence (Figure 11) which contributes to the phase offset
of about $1.7\text{\,}\mathrm{rad}$ compared to $a_{12}=0$ with one dipole
oscillation (red dashed line with open circles).

Figure 13: Different Ramsey fringes of the $\ket{+1}-\ket{+2}$ BEC system,
where the blue solid line with stars is for $a_{12}=0$ with immediate
recombination pulse, red dashed line with open circles is for $a_{12}=0$ with
recombination pulse at one dipole oscillation, black dashed line is for
$a_{12}=98a_{0}$ and green solid line with solid squares is for
$a_{12}=-50a_{0}$.

## Discussion

The above shows intriguing results for the $\ket{+1}-\ket{+2}$ BEC system and
several contrast reduction mechanisms for the Ramsey signal. It should be
noted that, though the amplitude of the Ramsey signal is small, the integrity
of the signal over several dipole oscillations may be explored depending on
experimental findings. All in all, the above displays the capabilities of the
developed tool which is a good platform to explore and interpret experimental
data.

## References

  * [1] Deutsch, C. _et al._ Spin self-rephasing and very long coherence times in a trapped atomic ensemble.  _Physical Review Letters_ 105 (2010).  URL https://doi.org/10.1103/physrevlett.105.020401.
  * [2] Bernon, S. _et al._ Manipulation and coherence of ultra-cold atoms on a superconducting atom chip.  _Nature Communications_ 4 (2013).  URL https://doi.org/10.1038/ncomms3380.
  * [3] Sárkány, L., Weiss, P., Hattermann, H. & Fortágh, J.  Controlling the magnetic-field sensitivity of atomic-clock states by microwave dressing.  _Physical Review A_ 90 (2014).  URL https://doi.org/10.1103/physreva.90.053416.
  * [4] Anderson, R. P., Ticknor, C., Sidorov, A. I. & Hall, B. V.  Spatially inhomogeneous phase evolution of a two-component bose-einstein condensate.  _Physical Review A_ 80 (2009).  URL https://doi.org/10.1103/physreva.80.023603.
  * [5] Egorov, M. _et al._ Long-lived periodic revivals of coherence in an interacting bose-einstein condensate.  _Physical Review A_ 84 (2011).  URL https://doi.org/10.1103/physreva.84.021605.
  * [6] Opanchuk, B., Egorov, M., Hoffmann, S., Sidorov, A. I. & Drummond, P. D.  Quantum noise in three-dimensional BEC interferometry.  _EPL (Europhysics Letters)_ 97, 50003 (2012).  URL https://doi.org/10.1209/0295-5075/97/50003.
  * [7] Egorov, M. _et al._ Measurement of s-wave scattering lengths in a two-component bose-einstein condensate.  _Physical Review A_ 87 (2013).  URL https://doi.org/10.1103/physreva.87.053614.
  * [8] Petrovic, J., Herrera, I., Lombardi, P., Schäfer, F. & Cataliotti, F. S.  A multi-state interferometer on an atom chip.  _New Journal of Physics_ 15, 043002 (2013).  URL https://doi.org/10.1088/1367-2630/15/4/043002.
  * [9] Thenuwara, A. & Sidorov, A.  Ramsey interferometry in three-level and five-level systems of ${}^{87}rb$ bose-einstein condensates (2023).  2308.11095.
  * [10] Rabi, I. I., Millman, S., Kusch, P. & Zacharias, J. R.  The molecular beam resonance method for measuring nuclear magnetic moments. the magnetic moments ofLi63, li73andf199.  _Physical Review_ 55, 526–535 (1939).  URL https://doi.org/10.1103/physrev.55.526.
  * [11] Griffiths, D. J.  _Inroduction to Quantum Mechanics, 2nd ed._ ((Pearson Prentice Hall, 2005).
  * [12] Ramsey, N. F.  Experiments with separated oscillatory fields and hydrogen masers.  _Science_ 248, 1612–1619 (1990).  URL https://doi.org/10.1126/science.248.4963.1612.
  * [13] Fang, F. _et al._ NIM5 cs fountain clock and its evaluation.  _Metrologia_ 52, 454–468 (2015).  URL https://doi.org/10.1088/0026-1394/52/4/454.
  * [14] Heavner, T. P., Jefferts, S. R., Donley, E. A., Shirley, J. H. & Parker, T. E.  NIST-f1: recent improvements and accuracy evaluations.  _Metrologia_ 42, 411–422 (2005).  URL https://doi.org/10.1088/0026-1394/42/5/012.
  * [15] Peters, A., Chung, K. Y. & Chu, S.  Measurement of gravitational acceleration by dropping atoms.  _Nature_ 400, 849–852 (1999).  URL https://doi.org/10.1038/23655.
  * [16] Rosi, G., Sorrentino, F., Cacciapuoti, L., Prevedelli, M. & Tino, G. M.  Precision measurement of the newtonian gravitational constant using cold atoms.  _Nature_ 510, 518–521 (2014).  URL https://doi.org/10.1038/nature13433.
  * [17] Feynman, R. P., Vernon, F. L. & Hellwarth, R. W.  Geometrical representation of the schrödinger equation for solving maser problems.  _Journal of Applied Physics_ 28, 49–52 (1957).  URL https://doi.org/10.1063/1.1722572.
  * [18] Vitanov, N. V. _et al._ Fault-tolerant hahn-ramsey interferometry with pulse sequences of alternating detuning.  _Physical Review A_ 91 (2015).  URL https://doi.org/10.1103/physreva.91.033406.
  * [19] Gross, E. P.  Structure of a quantized vortex in boson systems.  _Il Nuovo Cimento_ 20, 454–477 (1961).  URL https://doi.org/10.1007/bf02731494.
  * [20] _Bose-Einstein Condensation_ (Oxford University Press, 2003).
  * [21] Kiesewetter, S., Polkinghorne, R., Opanchuk, B. & Drummond, P. D.  xSPDE: Extensible software for stochastic equations.  _SoftwareX_ 5, 12–15 (2016).  URL https://doi.org/10.1016/j.softx.2016.02.001.
  * [22] Bao, W., Jaksch, D. & Markowich, P. A.  Numerical solution of the gross–pitaevskii equation for bose–einstein condensation.  _Journal of Computational Physics_ 187, 318–342 (2003).  URL https://doi.org/10.1016/s0021-9991(03)00102-5.
  * [23] Bao, W. & Tang, W.  Ground-state solution of bose–einstein condensate by directly minimizing the energy functional.  _Journal of Computational Physics_ 187, 230–254 (2003).  URL https://doi.org/10.1016/s0021-9991(03)00097-4.
  * [24] Antoine, X. & Duboscq, R.  GPELab, a matlab toolbox to solve gross–pitaevskii equations i: Computation of stationary solutions.  _Computer Physics Communications_ 185, 2969–2991 (2014).  URL https://doi.org/10.1016/j.cpc.2014.06.026.
  * [25] Antoine, X. & Duboscq, R.  GPELab, a matlab toolbox to solve gross–pitaevskii equations II: Dynamics and stochastic simulations.  _Computer Physics Communications_ 193, 95–117 (2015).  URL https://doi.org/10.1016/j.cpc.2015.03.012.
  * [26] Minguzzi, A.  Numerical methods for atomic quantum gases with applications to bose–einstein condensates and to ultracold fermions.  _Physics Reports_ 395, 223–355 (2004).  URL https://doi.org/10.1016/j.physrep.2004.02.001.
  * [27] Bao, W., Jin, S. & Markowich, P. A.  On time-splitting spectral approximations for the schrödinger equation in the semiclassical regime.  _Journal of Computational Physics_ 175, 487–524 (2002).  URL https://doi.org/10.1006/jcph.2001.6956.
  * [28] Bao, W.  Ground states and dynamics of multicomponent bose–einstein condensates.  _Multiscale Modeling & Simulation_ 2, 210–236 (2004).  URL https://doi.org/10.1137/030600209.
  * [29] Hall, D. S., Matthews, M. R., Ensher, J. R., Wieman, C. E. & Cornell, E. A.  Dynamics of component separation in a binary mixture of bose-einstein condensates.  _Physical Review Letters_ 81, 1539–1542 (1998).  URL https://doi.org/10.1103/physrevlett.81.1539.
  * [30] Bao, W. & Du, Q.  Computing the ground state solution of bose–einstein condensates by a normalized gradient flow.  _SIAM Journal on Scientific Computing_ 25, 1674–1697 (2004).  URL https://doi.org/10.1137/s1064827503422956.
  * [31] Antoine, X., Bao, W. & Besse, C.  Computational methods for the dynamics of the nonlinear schrödinger/gross–pitaevskii equations.  _Computer Physics Communications_ 184, 2621–2633 (2013).  URL https://doi.org/10.1016/j.cpc.2013.07.012.
  * [32] Bao, W., & and, Y. C.  Mathematical theory and numerical methods for bose-einstein condensation.  _Kinetic & Related Models_ 6, 1–135 (2013).  URL https://doi.org/10.3934/krm.2013.6.1.
  * [33] Bao, W.  Mathematical models and numerical methods for bose-einstein condensation.  _Proceedings of the International Congress for Mathematicians_ 4, 1–26 (2014).
  * [34] Besse, C.  A relaxation scheme for the nonlinear schrödinger equation.  _SIAM Journal on Numerical Analysis_ 42, 934–952 (2004).  URL https://doi.org/10.1137/s0036142901396521.
  * [35] Tacla, A. B., Boixo, S., Datta, A., Shaji, A. & Caves, C. M.  Nonlinear interferometry with bose-einstein condensates.  _Physical Review A_ 82, 053636–1–053636–9 (2010).  URL https://doi.org/10.1103/physreva.82.053636.
  * [36] Söding, J. _et al._ Three-body decay of a rubidium bose–einstein condensate.  _Applied Physics B_ 69, 257–261 (1999).  URL https://doi.org/10.1007/s003400050805.
  * [37] Widera, A. _et al._ Precision measurement of spin-dependent interaction strengths for spin-1 and spin-287rb atoms.  _New Journal of Physics_ 8, 152–152 (2006).  URL https://doi.org/10.1088/1367-2630/8/8/152.


# Speaker anonymisation using the McAdams coefficient

###### Abstract

Anonymisation has the goal of manipulating speech signals in order to degrade
the reliability of automatic approaches to speaker recognition, while
preserving other aspects of speech, such as those relating to intelligibility
and naturalness. This paper reports an approach to anonymisation that, unlike
other current approaches, requires no training data, is based upon well-known
signal processing techniques and is both efficient and effective. The proposed
solution uses the McAdams coefficient to transform the spectral envelope of
speech signals. Results derived using common VoicePrivacy 2020 databases and
protocols show that random, optimised transformations can outperform competing
solutions in terms of anonymisation while causing only modest, additional
degradations to intelligibility, even in the case of a semi-informed privacy
adversary.

Index Terms: anonymisation, pseudonymisation, privacy, de-identification,
automatic speaker recognition

##  1 Introduction

Recent years have seen an increase in privacy legislation. Much of it covers
what is referred to as _personally identifiable information_ (PII), e.g.
biometric data such as speech [1, 2]. The community-led VoicePrivacy
initiative [3] aims to foster progress in the development of anonymisation
techniques that can be employed to conceal PII contained within speech
signals.

The VoicePrivacy 2020 evaluation plan [4, Sec. 3.2] specifies four
requirements for successful anonymisation solutions. They should: (i) produce
a speech waveform; (ii) suppress speaker-specific information as much as
possible; (iii) preserve intelligibility and naturalness; (iv) protect voice
distinctiveness. Anonymised speech should degrade the reliability of
identification applications such as automatic _speaker_ recognition, but
should not interfere with the correct functioning of utility applications such
as automatic _speech_ recognition. Anonymised speech recordings originating
from the same source speaker should have the same voice, while the linking of
anonymised voices to source speakers should be prevented. Given that each
speaker should retain a distinctive, but anonymised voice, the task is perhaps
referred to more appropriately as pseudonymisation [5].

Best practice dictates that assessment is applied according to different
attack models or scenarios whereby a privacy adversary has knowledge of, or
access to, a similarly (but not identically configured) anonymisation system.
To support assessment according to the various VoicePrivacy scenarios,
anonymisation is required to be applied separately, with different
configurations, to two different data subsets, enrolment and test, such that
the anonymised voice for each speaker in each set is different. Further
details concerning the attack models are discussed in [4, Sec. 3.3].

Two baseline anonymisation systems111https://github.com/Voice-Privacy-
Challenge/Voice-Privacy-Challenge-2020 were made available to VoicePrivacy
2020 participants: (i) a primary baseline based upon state-of-the-art x-vector
embeddings and neural waveform techniques [6]; (ii) a secondary baseline,
inspired from the McAdams coefficient [7], consisting of well-known signal
processing techniques [8]. The primary baseline is comparatively complex and
requires substantial training data and computational resources. Based upon a
simple contraction or expansion of pole locations derived using linear
predictive coding (LPC), the secondary baseline requires no training data and
is comparatively straightforward and efficient. While the primary baseline
better suppresses speaker specific information [3] (requirement ii), the
secondary baseline better preserves voice distinctiveness [5] (requirement
iv). The secondary baseline was made available to provide additional
inspiration, to expose a broader potential solution space and to lower the
cost of entry to the VoicePrivacy initiative. The intention was to help
attract potential participants with a signal processing background that might
lack the background in speaker characterisation required to explore the
primary baseline. The secondary baseline was not optimised in any way (other
than through casual listening tests), was applied in a deterministic fashion
such that anonymisation is easily reversible and was adopted by others [9,
10].

The work reported in this paper aims (i) to explore more thoroughly the
potential of well-known signal processing techniques as a solution to the
anonymisation problem, and (ii) to assess the real need for and benefit of
more complex, demanding solutions. The paper reports our efforts to optimise
the original McAdams-based solution and its adaptation to a more stochastic
approach which affords better protection from reversibility. Longer-term, it
is hoped that this work might expose opportunities to improve performance by
combining the components or techniques used by each baseline solution.

The remainder of this paper is organised as follows. Related past work is
described in Section 2. The proposed approach to anonymisation is introduced
in Section 3. Experiments are reported in Section 4 while a discussion of
findings, conclusions and ideas for future work are presented in Section 5.

##  2 Previous work

A popular synthetic sound generation techniques in the field of music signal
processing is that of additive synthesis [11]. The technique is used to
generate timbre through the addition of multiple cosinusoidal oscillations:

$y(t)=\sum_{k=1}^{K}r_{k}(t)\operatorname{cos}(2\pi(kf_{0})^{\alpha}t+\phi_{k})$
(1)

where $k$ is the harmonic index, $r_{k}(t)$ is amplitude, $\phi_{k}$ is the
phase, $t$ is time and $\alpha$ is the so-called McAdams coefficient [7].
Equation 1 resembles the synthesis of a periodic signal with an inverse
Fourier series which combines harmonic cosinusoidal oscillations, each with
some magnitude and some phase shift.

The McAdams coefficient is used to adjust the frequency of each harmonic.
Adjustments to the distribution of harmonics acts to modify the resulting
timbre. When applied to additive synthesis as in (1), the McAdams coefficient
transforms harmonics to non-harmonic components which are also referred to as
_partials_ or overtones.

This technique bears similarities to other, well-known techniques in speech
processing [12, 13, 14, 15, 16, 17]. One example is vocal tract length
normalization (VTLN) [18, 19] which has been explored previously in
anonymisation-related work [20, 21]. Warping functions are applied to the
spectral log amplitude of each speech frame to modify not just the envelope,
as in our approach described in Section 3, but also the pitch. It can hence be
argued that these approaches offer greater anonymisation strength and that our
approach may then be inferior where identification is performed by human
listeners. However, conversion of the pitch may be somewhat redundant if
anonymisation is to be performed purely for the purposes of deceiving
_automatic_ approaches to identification which rarely use estimates of pitch
or source/residual signals.

##  3 Anonymisation using
the McAdams coefficient

Figure 1: Pipeline of the application of the proposed McAdams coefficient-
based approach to anonymisation on a speech frame basis. The angle $\phi$ of
poles with a non-zero imaginary part are raised to the power of the McAdams
coefficient $\alpha$ to provoke a expansion/contraction in frequency in its
associated formant.

Our approach to pseudonymisation is explained in this section. It is based
upon the idea of applying a shift to the formant positions in a speech
utterance, thereby adjusting the timbre or spectral envelope. Representations
of the latter are often used in the front-end of an automatic speaker
verification (ASV) system. Consequently, manipulations to the formant
positions should degrade ASV reliability and hence provide some level of
anonymisation. In our approach, the degree of formant shifting is controlled
by the McAdams coefficient [7] as in (1). The principal ideas are outlined
first, before we show the resulting system is applied to the anonymisation of
speech signals.

###  3.1 Application to speech signals

We use the McAdams coefficient to manipulate the formant positions of speech
signals at the frame level. The system described below is that which we
proposed as the secondary baseline for the VoicePrivacy 2020 challenge. The
process is illustrated in Figure 1. First, source-filter analysis is applied
to each frame using linear predictive coding (LPC). The source or residual is
set aside for later resynthesis, whereas filter coefficients are used to
derive the set of pole positions. Real-valued pole positions (with zero-valued
imaginary terms) are left unmodified, whereas complex-valued poles (with non-
zero imaginary terms) are shifted according to the higher branch in Figure 1.

The shift, applied via the same transformation in (1), operates on the angle
between the positive real axis and the vector extending from the origin in the
z-plane to the complex pole position. This angle corresponds to frequency,
with the upper half of the unit circle (angle of $\pi$ radians) corresponding
to the sampling frequency. The full set of complex-valued pole positions is
shifted according to $\phi^{\alpha}$ resulting in either clockwise or counter-
clockwise shifts in pole locations.

This result of this process is illustrated in the z-plane of Figure 2 where
the blue circles indicate original pole locations and where other points
indicate induced shifts. As shown, the shift is different for each pole in
terms of both direction and the scale of the shift. For poles with an angle of
$\phi<1$ radian, the shift $\phi^{\alpha}$ is counter-clockwise for values of
$\alpha<1$. While not explored in this work, values of $\alpha>1$ would induce
shifts in a clockwise direction for values of $\phi<1$ radian. Poles with an
angle of $\phi>1$ radians undergo shifts in the opposite directions. Lastly,
the scale of the shift is dependent upon the gap between $\phi$ and $\phi=1$
radian; the greater the gap, the greater the shift. Example shifts in pole
positions are illustrated in Figure 2 for values of $\alpha<1$ for values
$\alpha\in\\{0.9,0.7,0.5\\}$.

An illustration of the result in terms of spectral envelope is shown in Figure
3. Four spectra are illustrated. The spectra illustrated by the solid blue
line is the original speech spectrum whereas the others show the spectra after
anonymisation for values of $\alpha\in\\{0.9,0.7,0.5\\}$. The nearer a pole to
the arc of the unit circle (Figure 1), the more apparent the corresponding
peak in the spectral envelope, with each pole and peak loosely corresponding
to a speech formant. Given a sampling rate of 16 kHz (as for VoicePrivacy 2020
data), the value of $\phi=1$ radian corresponds to a frequency of
approximately 2.5 kHz. The shift in pole positions is hence in opposite
directions for frequencies either side of this threshold and is either an
expansion away from 2.5 kHz or a contraction towards it.

Returning to Figure 1, the set of new poles, including original real-valued
poles and those modified according to the format shifting procedure, are then
converted back to LPC coefficients. The latter are combined with the residual
from the original speech signal and then used to resynthesise an anonymised
speech frame in the time domain. The stream of speech frames, each treated
according to the procedure described above, is then combined by means of a
standard overlap and add (OLA) technique to produce the final anonymised
speech signal.

Figure 2: An illustrative example of original pole positions and pole shifts
for values of $\alpha\in\\{0.9,0.7,0.5\\}$. Only the top-half of the pole-zero
map is shown. Complex conjugate poles undergo the same shift, but in the
opposite direction.

###  3.2 Stochastic anonymisation

The original VoicePrivacy baseline used an arbitrary, fixed McAdams
coefficient of $\alpha=0.8$. Clearly, there is scope for optimisation since
different values further away from $\alpha=1$ will potentially induce greater
shifts in the pole positions and should hence provide greater levels of
anonymisation (higher EERs). However, better anonymisation will likely come at
the expense of speech distortion, with the latter likely also causing
degradation to intelligibility and naturalness. Accordingly, we report new
work which aims to explore the trade off between these two competing
objectives. We restrict the study to values of $\alpha<1$. Casual listening
tests showed that the use of values of $\alpha>1$ resulted in less acceptable
degradations to intelligibility.

Next, rather than using a fixed, deterministic value of $\alpha$, which
implies that anonymisation is reversible, we sought to confirm that different
values of $\alpha$ can produce different pseudo-voices for the same speaker.
The goal here is to support experimentation with different attack models and
scenarios by ensuring that the anonymisation of enrolment and test data
results in _different_ , or sufficiently dissimilar pseudo-voices (see Section
1 and [4, Sec. 3.2]). The use of a fixed $\alpha$ in the original baseline
meant that the anonymisation of enrolment and test utterances previously
resulted in the _same_ , or similar pseudo-voices.

The result is a stochastic approach to anonymisation whereby the McAdams
coefficient is sampled from within a range of the uniform distribution, i.e.
$\alpha\in\text{U}(\alpha_{\text{min}},\alpha_{\text{max}})$. The application
of anonymisation to enrolment and test data with different, random McAdams
coefficients drawn from different ranges of the same uniform distribution
should result in different pseudo-voices for each speaker’s enrolment and test
utterances, while also providing some protection from reversibility. We note
that anonymisation is applied on a speaker-dependent basis, with each speaker
having a randomly chosen McAdams coefficient. A privacy adversary would then
need to know the exact McAdams coefficient used to anonymise the speech of any
particularly speaker in order to revert the transform.

##  4 Experiments

This section describes the VoicePrivacy 2020 database used in this work, the
metrics used for assessment, the different attack models or scenarios and our
results.

###  4.1 Data and metrics

We used the VoicePrivacy 2020 database and protocols described in [3, 4].
However, in contrast to the primary baseline [6, 22], our anonymisation system
has no requirement for training data. Consequently, our solution makes no use
of the training partition. We report results for evaluation data only, which
is drawn from the LibriSpeech [23] and VCTK [24] source datasets. While we
have observed differences in performance for each dataset, as well as gender
dependencies, for reasons of space limitations we report only average
performance here.222Full results available at
https://github.com/josepatino/Voice-Privacy-
Challenge-2020/blob/master/results/ Database, gender dependencies and other
potential biases, obviously warrant investigation in further work. Assessment
is performed using the standard VoicePrivacy 2020 ASV and ASR systems [3] both
trained using the train-clean-360 partition of the LibriSpeech database.

Degradation to intelligibility is measured through the ASR word error rate
(WER) whereas anonymisation performance is measured through the ASV equal
error rate (EER). Also reported here are estimates of privacy using the zero
evidence biometric recognition assessment (ZEBRA) framework [25], a recently
proposed, adversary-agnostic, metric inspired by the work of the forensic
sciences community.

Figure 3: Spectral envelopes for an original speech frame and pole positions,
and for anonymised versions for $\alpha\in\\{0.9,0.7,0.5\\}$. Anonymisation
results in an expansion or contraction of the spectrum either side of $\phi=1$
radiant or 2.5 kHz.

###  4.2 Scenarios

In adopting the terminology of [21], we investigated anonymisation performance
according to two different attack models or scenarios. Both involve a privacy
adversary who seeks to determine whether two utterances belong, or not, to the
same speaker. In both cases there is a test utterance which is anonymised, and
an enrolment utterance which is either an original utterance, or a similarly-
anonymised utterance.

The first scenario assumes an  _ignorant_ privacy adversary who has no
knowledge of the test utterance being anonymised. In this case the adversary
attempts to determine the match between an original enrolment utterance and an
anonymised test utterance. This is referred to as the o-a scenario. The second
scenario assumes a _semi-informed_ privacy adversary who, while knowing that
the test utterance has been anonymised and having access to the systems used,
does not know the configuration. In this case, the adversary will likely
determine more reliably the match between the two utterances by anonymising
the enrolment utterance such that the comparison is now between two similarly
anonymised utterances, even though each is anonymised with differently
configured systems. It is referred to as the a-a scenario. Like [21], we do
not consider the unlikely scenario of an _informed_ adversary that has
knowledge of both the algorithm and the configuration.

We also report ASV and ASR results derived from systems that are retrained on
either original or anonymised data. These experiments aim to determine
robustness in the case that the adversary is able to harness knowledge of the
anonymisation system to break the protection it affords. Similarly, we seek to
determine the improvements one can obtain in intelligibility, if the ASR
system is adapted to anonymised data. See [4] for further details.

###  4.3 Results

Figure 4: ASV performance in terms of the equal error rate (EER, %) for
various systems and for the VoicePrivacy 2020 test set. B1 is the primary
baseline as in [3, 6, 22], while B2 is the secondary baseline as in [8] with a
fixed $\alpha=0.8$. Other cases are for the proposed solution with
$\alpha\in\text{U}(\alpha_{min},\alpha_{max}$).

Anonymisation results in terms of ASV EER are illustrated in Figure 4. Results
are shown for an original setup without anonymisation, for the primary and
secondary VoicePrivacy 2020 baselines (B1 and B2) and for B2 stochastic
variants where the McAdams coefficient is drawn from different intervals in
$\alpha\in\mbox{U}(\alpha_{min},\alpha_{max})$. All systems based on the
McAdams coefficient have an LPC order of 20.

It is clear that all systems increase the EER from the original 3.3%. Blue,
left-most bars show EERs for ignorant privacy adversaries (o-a scenario) when
test utterances are anonymised but enrolment utterances are not. Results show
that B1 out-performs B2 and all related systems by a substantial margin.
Yellow, middle bars show performance for semi-informed privacy adversaries
(a-a scenario) when both enrolment and test utterances are anonymised (with
the same, but differently configured algorithm). Now, B2 and related systems
are more competitive, e.g. 33.8% for U(0.7,0.9) _cf._ 32.3% for B1.

While lower values of $\alpha_{min}$ give better anonymisation, improvements
are accompanied by degradations to intelligibility. ASR results for the same
setups are illustrated in Figure 5. An original, baseline WER of 8.45%
increases only to 10.95% for B1, but to 13.7% for U(0.8,0.9) and worse for
other configurations. Thus, while the McAdams approach to anonymisation can
better protect privacy, it does so at the cost of reduced intelligibility.

The question then is what level of improvement to intelligibility can be
delivered simply by retraining the ASR systems with anonymised speech and what
would be the corresponding impacts to ASV. Answers are provided by the green
bars in Figures 4 and 5 for ASV and ASR systems respectively. For ASV, the
McAdams approach to anonymisation performs best, e.g. 22.6% for U(0.8,0.9)
_cf._ 13.7% for B1. Now, though, while better anonymisation is achieved for
$\alpha_{min}<0.8$, the corresponding WER is much more competitive with the
degradation for B1. Whereas for the latter, the WER decreases from 10.95% to
8.25%, it decreases from 45% to 9.6% for U(0.5,0.9), while the EER is almost
twice as high than for B1, e.g. up to 37.5% for U(0.5,0.9).

Table 1 shows ZEBRA results which reflect the evidence remaining to a privacy
adversary after anonymisation. They correspond to the semi-informed scenario
where the ASV system is retrained on anonymised data, i.e. the green bars in
Figures 4. While the scales are inverted (less bits of disclosure infer better
privacy), ZEBRA results shown in column 2 confirm much the same trend shown by
EER results. Noting that anonymisation affords different levels of privacy to
different individuals, the ZEBRA framework also allows one to determine the
_worst case_ level of privacy disclosure. According to [25], this is expressed
according to categorical tags where a tag of ‘0’ reflects no privacy
disclosure, ’A’ is the next best case and where tag ‘F’ reflects the worst.
Identical categorical tags of ‘C’ illustrated in column 3 for each system
suggest that the differences in performance discussed above are not so great
from the perspectives of a worst case scenario.

Figure 5: As for Fig. 4 except for the ASR performance in terms of the word
error rate (WER, %).  Table 1: ZEBRA results for the semi-informed (a-a)
scenario and where ASV systems are trained using anonymised data.

| System
---
| Expected privacy
---
disclosure (bits)
| Worst case
---
categorical tag
Original | 0.647 | C
B1 | 0.403 | C
B2 | 0.598 | C
$\alpha\in(0.8,0.9)$ | 0.366 | C
$\alpha\in(0.7,0.9)$ | 0.261 | C
$\alpha\in(0.6,0.9)$ | 0.199 | C
$\alpha\in(0.5,0.9)$ | 0.155 | C

##  5 Conclusions

This paper shows that well-known signal processing techniques can provide
efficient and effective solutions to anonymisation for speech signals. The
approach reported in this paper, based upon transformations to the spectral
envelope using the McAdams coefficient, can increase by ten-fold the equal
error rate of a baseline automatic speaker recognition system, with only
modest degradations to intelligibility, even in a semi-informed privacy
adversary scenario. This is achieved with a randomised transformation which
provides some level of protection from reversibility.

Despite this encouraging result, we must acknowledge that, until error rates
increase to the equivalent of random performance, or a ZEBRA category of ‘0’,
we remain far from meeting the goal of _true_ anonymisation. Database, gender
and other biases mean that an ‘averaged’ view of privacy tells an incomplete
picture, especially when we know that even the stronger anonymisation
solutions still leave some with relatively weaker protection. Encouragingly,
there is plenty of scope to extend this work. Thus far we have considered only
adjustments to the formant positions in terms of frequency. Future work should
explore more elaborate transformations to the spectral envelope which could be
investigated through random pole perturbations. There is also scope to explore
the resilience of the approach to brute-force attacks, and new approaches
which combine the merits of our approach with those of the more sophisticated
VoicePrivacy primary baseline based upon state-of-the-art x-vector embeddings
and neural waveform techniques.

Last, all results reported in this paper are reproducible with open source
code and scripts available online333https://github.com/josepatino/Voice-
Privacy-Challenge-2020/.

##  6 Acknowledgements

This work is partly funded by the VoicePersonae project which is supported by
the French Agence Nationale de la Recherche (ANR) and the Japan Science and
Technology Agency (JST). It is also linked to the VoicePrivacy initiative and
the Harpocrates project also funded by the ANR.

## References

  * [1] A. Nautsch, C. Jasserand, E. Kindt, M. Todisco, I. Trancoso, and N. Evans, “The GDPR & speech data: Reflections of legal and technology communities, first steps towards a common understanding,” in _Proc. INTERSPEECH_ , 2019\.
  * [2] A. Nautsch, A. Jiménez, A. Treiber, J. Kolberg, C. Jasserand, E. Kindt, H. Delgado, M. Todisco, M. A. Hmani, A. Mtibaa _et al._ , “Preserving privacy in speaker and speech characterisation,” _Computer Speech & Language_, vol. 58, pp. 441–480, 2019.
  * [3] N. Tomashenko, B. M. L. Srivastava, X. Wang, E. Vincent, A. Nautsch, J. Yamagishi, N. Evans, J. Patino, J.-F. Bonastre, P.-G. Noé _et al._ , “Introducing the VoicePrivacy initiative,” in _Proc. INTERSPEECH_ , 2020.
  * [4] ——, “The VoicePrivacy 2020 Challenge evaluation plan,” 2020.
  * [5] P.-G. Noé, J.-F. Bonastre, D. Matrouf, N. Tomashenko, A. Nautsch, and N. Evans, “Speech Pseudonymisation Assessment Using Voice Similarity Matrices,” in _Proc. INTERSPEECH_ , 2020.
  * [6] F. Fang, X. Wang, J. Yamagishi, I. Echizen, M. Todisco, N. Evans, and J.-F. Bonastre, “Speaker Anonymization Using X-vector and Neural Waveform Models,” in _Proc. 10th ISCA Speech Synthesis Workshop_ , 2018, pp. 155–160.
  * [7] S. McAdams, “Spectral fusion, spectral parsing and the formation of the auditory image,” _Ph. D. Thesis, Stanford_ , 1984.
  * [8] J. Patino, M. Todisco, A. Nautsch, and N. Evans, “Speaker anonymisation using the McAdams coefficient,” Eurecom, Tech. Rep. RR-20-343. 2020 [Online]. Available: http://www.eurecom.fr/publication/6190, Tech. Rep., 2020.
  * [9] P. Gupta, G. P. Prajapati, S. Singh, M. R. Kamble, and H. A. Patil, “Design of voice privacy system using linear prediction,” in _Proc. APSIPA_. IEEE, 2020, pp. 543–549.
  * [10] H. Kai, S. Takamichi, S. Shiota, and H. Kiya, “Lightweight voice anonymization based on data-driven optimization of cascaded voice modification modules,” in _Proc. IEEE SLT_ , 2021.
  * [11] C. Dodge and T. A. Jerse, _Computer Music: Synthesis, Composition and Performance_ , 2nd ed. Macmillan Library Reference, 1997.
  * [12] Q. Jin, A. R. Toth, T. Schultz, and A. W. Black, “Speaker de-identification via voice transformation,” in _Proc. ASRU_. IEEE, 2009, pp. 529–533.
  * [13] M. Pobar and I. Ipšić, “Online speaker de-identification using voice transformation,” in _2014 37th International convention on information and communication technology, electronics and microelectronics (mipro)_. IEEE, 2014, pp. 1264–1267.
  * [14] F. Bahmaninezhad, C. Zhang, and J. H. Hansen, “Convolutional neural network based speaker de-identification.” in _Proc. Odyssey_ , 2018, pp. 255–260.
  * [15] J. Qian, H. Du, J. Hou, L. Chen, T. Jung, X.-Y. Li, Y. Wang, and Y. Deng, “Voicemask: Anonymize and sanitize voice input on mobile devices,” _arXiv preprint arXiv:1711.11460_ , 2017.
  * [16] C. Magariños, P. Lopez-Otero, L. Docio-Fernandez, E. Rodriguez-Banga, D. Erro, and C. Garcia-Mateo, “Reversible speaker de-identification using pre-trained transformation functions,” _Computer Speech & Language_, vol. 46, pp. 36–52, 2017.
  * [17] B. M. L. Srivastava, A. Bellet, M. Tommasi, and E. Vincent, “Privacy-preserving adversarial representation learning in ASR: Reality or illusion?” _arXiv preprint arXiv:1911.04913_ , 2019.
  * [18] J. Cohen, T. Kamm, and A. G. Andreou, “Vocal tract normalization in speech recognition: Compensating for systematic speaker variability,” _The Journal of the Acoustical Society of America_ , vol. 97, no. 5, pp. 3246–3247, 1995.
  * [19] L. Lee and R. Rose, “A frequency warping approach to speaker normalization,” _IEEE Transactions on speech and audio processing_ , vol. 6, no. 1, pp. 49–60, 1998.
  * [20] J. Qian, H. Du, J. Hou, L. Chen, T. Jung, and X.-Y. Li, “Hidebehind: Enjoy Voice Input with Voiceprint Unclonability and Anonymity,” in _Proc. ACM Conference on Embedded Networked Sensor Systems_ , 2018, pp. 82–94.
  * [21] B. M. L. Srivastava, N. Vauquier, M. Sahidullah, A. Bellet, M. Tommasi, and E. Vincent, “Evaluating voice conversion-based privacy protection against informed attackers,” in _Proc. ICASSP_. IEEE, 2020, pp. 2802–2806.
  * [22] B. M. L. Srivastava, N. Tomashenko, X. Wang, E. Vincent, J. Yamagishi, M. Maouche, A. Bellet, and M. Tommasi, “Design choices for x-vector based speaker anonymization,” in _Proc. INTERSPEECH_ , 2020.
  * [23] V. Panayotov, G. Chen, D. Povey, and S. Khudanpur, “Librispeech: an ASR corpus based on public domain audio books,” in _Proc. ICASSP_. IEEE, 2015, pp. 5206–5210.
  * [24] C. Veaux, J. Yamagishi, K. MacDonald _et al._ , “CSTR VCTK corpus: English multi-speaker corpus for CSTR Voice Cloning Toolkit,” 2016.
  * [25] A. Nautsch, J. Patino, N. Tomashenko, J. Yamagishi, P.-G. Noe, J.-F. Bonastre, M. Todisco, and N. Evans, “The Privacy ZEBRA: Zero Evidence Biometric Recognition Assessment,” in _Proc. INTERSPEECH_ , 2020.