###### Contents

  1. 1 Introduction
  2. 2 Particle distributions
    1. 2.1 Particles as Dirac delta
    2. 2.2 Isotropic distributions
  3. 3 Particles in two regions
    1. 3.1 Isotropic distributions
  4. 4 The discrete form for isotropic pair-correlations
    1. 4.1 Particles in one region
    2. 4.2 Particles in two regions
  5. 5 The structure factor
  6. 6 Particle configurations from the structure factor
    1. 6.1 Restrictions
    2. 6.2 Gradient optimisation
    3. 6.3 Preliminary numerical results
  7. 7 Discussion

# Calculating Pair Correlations from Random Particle Configurations

Aristeidis Karnezis  Art L. Gower Website: https://arturgower.github.io

(28th August 2024)

## Abstract

Particle pair-correlations are broadly used to describe particle distributions
in chemistry, physics, and material science. Many theoretical methods require
the pair-correlation to predict material properties such as fluid flow,
thermal properties, or wave propagation. In all these applications it is
either important to calculate a pair-correlation from specific particle
configurations, or vice-versa: determine the likely particle configurations
from a pair-correlation which is needed to fabricate a particulate material.
Most available methods to calculate the pair-correlation from a particle
configuration require that the configuration be very large to avoid effects
from the boundary. Here we show how to avoid boundary effects even for small
particle configurations. Having small particle configurations leads to far
more efficient numerical methods. We also demonstrate how to use techniques
from smooth nonlinear optimisation to quickly recover a particle configuration
from a pair-correlation.

##  1 Introduction

#### Background.

In the field of material science, the investigation of structural properties
has witnessed significant progress. Molecular dynamics simulations, coupled
with neutron and X-ray scattering experiments, have provided valuable insights
into the structural characteristics of materials. For disordered materials,
most techniques to probe and analyse these structures focus on the pair-
correlation function and the structure factor [35, 34, 26, 25, 39]. The
structure factor appears naturally from scattering experiments [8, 11], for
example with small-angle neutron scattering [38, 7]. Pair-correlations and
structure factors are also used when calculating shear viscosity, electrical
conductivity, and thermal conductivity from a molecular perspective [28].

From the authors background, the pair-correlation appears when taking an
ensemble average of waves in disordered materials [14, 15, 6, 20]. That is,
the pair-correlation is the only way the material structure affects wave
propagation. This suggestions a route to design materials to control waves:
1) to determine pair-correlations that lead to band-gaps, frequency filters,
or enhanced transmission, and then
2) determine configurations of particles that match the desired pair-
correlation. Step 2) is known as the realizability problem [52, 21, 10, 11,
46], and there are open questions about when it is possible to solve [33, 12,
13, 9, 19, 26].

#### Realizability problem.

Determining a likely configuration of particles from a specific pair-
correlation is called the realizability problem. It appears in the study of
many-body systems, such as liquids, and disordered materials [39, 46, 21, 52,
22, 10, 40, 11, 50]. One of the challenges is that it is not known when it is
possible to find a configuration of particles for any given pair-correlation
[11, 41]. While several necessary conditions have been identified for
potential pair-correlation functions, including non-negativity and
restrictions on their associated structure factors [39], achieving a set of
sufficient conditions remains an open challenge. We note that the problem is
not completely resolved even in just one spatial dimension [39].

#### Further Applications.

There is a wide range of applications of determining particle configurations
from structure factor. One example is to better understand the enigmatic
"nuclear pasta" phenomenon found in the extreme conditions of neutron star
crusts [32, 17, 16]. Nuclear pasta refers to a non-uniform arrangements of
subatomic particles, in particular neutrons and protons, within the intense
gravitational fields of neutron stars. Neutron and X-ray scattering provide a
measure of the structure factor, from which researchers attempt to learn about
the internal structure of the nuclear pasta. So far, pasta-like formations
have been uncovered [32, 17, 16].

#### Reverse Monte-Carlo.

Reverse Monte-Carlo structural modeling is one technique that has been used to
calculate particle configurations that match a measured structure factors or
pair-correlations [33, 12, 13, 9, 19, 26, 44]. Typically these Monte-Carlo
simulations are guided by a genetic algorithm, or similar random searches,
which are computationally intensive [26, 25, 24, 36]. These methods are brute
force, and typically use non-gradient based methods. As the gradients can be
easily calculated, it is likely that gradient based methods [37] will out
perform non-gradient based methods.

#### Paper summary.

We start by deducing in a simple self-contained way how to describe
probability distributions and the pair-correlation in terms of a set of
particles Section 2. We then show how to avoid the affects of boundaries when
calculating pair-correlations in Section 3. The affects of boundaries are
usually undesirable, and the details we present seem to be missing from most
references. In Section 4 we show the same calculations again, but without the
use of Dirac deltas for didactic purposes.

In Section 5 we deduce the structure factor for isotropic distributions both
from any given pair-correlation and a set of particles, these results are need
in Section 6, where we present a method to calculate a configuration of
particles that matches a given structure factor. Developing more efficient
methods to reconstruct particle configurations from pair-correlations remains
an ongoing challenge in material science and computational chemistry [29]. The
method we propose uses techniques from smooth nonlinear optimisation to
improve the efficiency, which we are able to do because the structure-factor
is a smooth function of the particle positions. We also present some
preliminary numerical results. Finally, in Section 7 we summarise what the
paper achieved and possible future directions.

##  2 Particle distributions

Consider there are $J$ particles placed within some region $\mathcal{R}$. That
is, the centre of every particle $\bm{r}_{j}\in\mathcal{R}$. We represent one
possible position of the particles, or ensemble, by the set $\mathcal{X}^{s}$,
where every $\bm{r}\in\mathcal{X}^{s}$ is the centre of a particle in the
ensemble.

The function $\mathrm{p}(\bm{r})$ is probability density of finding a particle
centered at $\bm{r}$. We can approximate $\mathrm{p}(\bm{r})$ by defining a
mesh of volume elements $V(\bm{v}_{i})$, where the vector $\bm{v}_{i}$ is the
centre of the volume element, and then counting the number of particles in
each $V(\bm{v}_{i})$ divided by the total number of particles. This leads us
to define $\mathrm{p}_{V}(\bm{v}_{i})$ below to approximate
$\mathrm{p}(\bm{v}_{i})\approx\mathrm{p}_{V}(\bm{v}_{i})$ in the form

$\mathrm{p}_{V}(\bm{v}_{i})=\frac{1}{S}\sum_{s}^{S}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{i})\right]}{\\#\mathcal{X}^{s}}\frac{1}{|V(\bm{v}_{i})|},$ (1)

where $|V(\bm{v}_{i})|$ is the volume of $V(\bm{v}_{i})$, $\\#\mathcal{A}$ is
the number of elements in the set $\mathcal{A}$, and $S$ is the number of
ensembles considered.

We require that the mesh

$V(\bm{v}_{i})\cap V(\bm{v}_{j})=\varnothing,\;\;\text{for}\;\;j\not=i.$

In a similar way, we can approximate $\mathrm{p}(\bm{x}_{1},\bm{x}_{2})$,
which is the joint probability density of finding one particle centred at
$\bm{x}_{1}$ and another centred at $\bm{x}_{2}$, while averaging over all
over particle positions. If we assume that both $\bm{x}_{1}$ and $\bm{x}_{2}$
are distributed within $\mathcal{R}$, then we can approximate the function
$\mathrm{p}(\bm{x}_{1},\bm{x}_{2})$ with the formula:

$\mathrm{p}_{V}(\bm{v}_{i},\bm{v}_{j})=\frac{1}{S}\sum_{s}^{S}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{i})\right]}{\\#\mathcal{X}^{s}}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{j})\right]}{\\#\left[\mathcal{X}^{s}\backslash
V(\bm{v}_{i})\right]}\frac{1}{|V(\bm{v}_{i})||V(\bm{v}_{j})|}\quad\text{for}\;\;i\not=j,$
(2)

where $\mathrm{p}_{V}(\bm{v}_{i},\bm{v}_{i})=0$ for every $i$, and
$\mathcal{X}^{s}\backslash V(\bm{v}_{i})$ is defined as the set
$\mathcal{X}^{s}$ without $V(\bm{v}_{i})$. Note that

$\mathrm{p}(\bm{v}_{i},\bm{v}_{j})\approx\mathrm{p}_{V}(\bm{v}_{i},\bm{v}_{j}),$

for every $\bm{v}_{i}$ and $\bm{v}_{j}$.

As expected, using this approximation, the integral of
$\mathrm{p}(\bm{x}_{1},\bm{x})$ for $\bm{x}_{1},\bm{x}_{2}\in\mathcal{R}$
gives one:

$\int\mathrm{p}(\bm{x}_{1},{\bm{x}_{2}})\mathrm{d}\bm{x}_{1}\mathrm{d}\bm{x}_{2}\approx\sum_{i,j}\mathrm{p}_{V}(\bm{v}_{i},\bm{v}_{j})|V(\bm{v}_{i})||V(\bm{v}_{j})|\\\
=\frac{1}{S}\sum_{s}^{S}\sum_{i}\sum_{j\not=i}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{i})\right]}{\\#\mathcal{X}^{s}}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{j})\right]}{\\#\left[\mathcal{X}^{s}\backslash
V(\bm{v}_{i})\right]}\\\
=\frac{1}{S}\sum_{s}^{S}\sum_{i}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{i})\right]}{\\#\mathcal{X}^{s}}\frac{\\#\left[\mathcal{X}^{s}\cap\cup_{j\not=i}V(\bm{v}_{j})\right]}{\\#\left[\mathcal{X}^{s}\backslash
V(\bm{v}_{i})\right]}=\frac{1}{S}\sum_{s}^{S}\frac{\\#\left[\mathcal{X}^{s}\cap\cup_{i}V(\bm{v}_{i})\right]}{\\#\mathcal{X}^{s}}=1,$
(3)

where we used that $\cup_{i}V(\bm{v}_{i})=\mathcal{R}$,
$\mathrm{d}\bm{x}_{1}\approx|V(\bm{v}_{i})|$,$\mathrm{d}\bm{x}_{2}\approx|V(\bm{v}_{j})|$,
and $\mathcal{R}\cap\cup_{j\not=i}V(\bm{v}_{j})=\mathcal{R}\backslash
V(\bm{v}_{i})$ which implies that
$\mathcal{X}^{s}\cap\cup_{j\not=i}V(\bm{v}_{j})=\mathcal{X}^{s}\backslash
V(\bm{v}_{i})$.

The term that often appears in methods that use ensemble averaging of
particulates [4, 20, 43] is the particle pair-correlation, which is defined as

$\mathrm{g}(\bm{x}_{1},\bm{x}_{2}):=\frac{J-1}{J}\frac{\mathrm{p}(\bm{x}_{1},\bm{x}_{2})}{\mathrm{p}(\bm{x}_{1})\mathrm{p}(\bm{x}_{2})},$
(4)

where $J=\\#\mathcal{X}_{s}$ is the total number of particle, where we have
assumed that every configuration $\\#\mathcal{X}_{s}$ has the same number of
particles. The factor ${J-1}/{J}$ is correction for a finite number of
particles so that $\mathrm{g}(\bm{x}_{1},\bm{x}_{2})\to 1$ when particles
become uncorrelated, as confirmed by [20, Equation (8.1.2)]. We demonstrate
this in the next section.

###  2.1 Particles as Dirac delta

For a finite, but very large number, of particles we can rewrite the pair-
correlation in terms of Dirac deltas. For example, turning to (1) we assume
there is a finite number of ensembles $S$, each with a finite number of
particles $J$. Then we can make the volume elements $V(\bm{v}_{i})$ small
enough so that there is at most one particle in each volume element, which
implies that $\\#\mathcal{X}^{s}=J$ and

$\\#\left[\mathcal{X}^{s}\cap V(\bm{v}_{i})\right]=1,$

if there is a particle in $V(\bm{v}_{i})$ or otherwise the above is $0$. For
notational convenience, we choose the mesh of volume elements such that each
element is centred at a particle $\bm{r}_{i}\in\mathcal{X}^{s}$, if it
contains a particle. The above allows us to rewrite (1) in the form:

$\mathrm{p}_{V}(\bm{r}_{i})=\frac{1}{S}\sum_{s}^{S}\frac{1}{|V(\bm{r}_{i})|}\frac{1}{J},$
(5)

for every $\bm{r}_{i}\in\mathcal{X}^{s}$.

As we are taking the limit of all the volume elements going to zero,
$|V(\bm{r}_{i})|\to 0$, we can approximate:

$\mathrm{p}(\bm{x})=\begin{cases}\mathrm{p}_{V}(\bm{r}_{i})&\text{if}\;\;\bm{x}\in
V(\bm{r}_{i}),\\\ 0&\text{else},\end{cases}$

where $\mathrm{p}(\bm{x})=0$ if $\bm{x}$ is in a volume element that does not
contain a particle. With this definition, and taking the limit
$|V(\bm{v}_{i})|\to 0$ for every $i$, we can rewrite (1) in the form

$\mathrm{p}(\bm{x})=\frac{1}{JS}\sum_{s}^{S}\sum_{\bm{r}_{i}\in\mathcal{X}^{s}}\delta(\bm{x}-\bm{r}_{i}),$
(6)

where $\delta(\bm{x}-\bm{r}_{i})$ is the Dirac delta function. For details on
the Dirac delta see [2].

Repeating analogous steps for the joint probability density (2) we obtain

$\mathrm{p}(\bm{x}_{1},\bm{x}_{2})=\frac{1}{S}\frac{1}{J(J-1)}\sum_{s}^{S}\sum_{\bm{r}_{i}\in\mathcal{X}^{s}}\sum_{\bm{r}_{j}\not=\bm{r}_{i},\bm{r}_{j}\in\mathcal{X}^{s}}\delta(\bm{x}_{1}-\bm{r}_{i})\delta(\bm{x}_{2}-\bm{r}_{j}).$
(7)

###  2.2 Isotropic distributions

For most of this paper we focus on homogeneous and isotropic distributions.
The homogeneous assumptions leads to the particles are equally likely to be
placed anywhere in space, which leads to $\mathrm{p}(\bm{r})$ being a
constant. If we choose all volume elements to have the same volume
$|V(\bm{v}_{i})|=|V|$, then using that $\mathrm{p}_{V}(\bm{v}_{i})$ is a
constant together with (1) leads to

$\mathrm{p}_{V}(\bm{v}_{j})=\frac{1}{SI}\sum_{i}\sum_{s}^{S}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{i})\right]}{\\#\mathcal{X}^{s}}\frac{1}{|V|}=\frac{1}{SI}\frac{1}{|V|}\sum_{s}^{S}\frac{\\#\left[\mathcal{X}^{s}\cap\cup_{i}V(\bm{v}_{i})\right]}{\\#\mathcal{X}^{s}}=\frac{1}{I|V|}=\frac{1}{|\mathcal{R}|},$
(8)

where $I$ is the total number of volume elements $V(\bm{v}_{i})$,
$|\mathcal{R}|$ is the volume of $\mathcal{R}$, and we used
$\cup_{i}V(\bm{v}_{i})=\mathcal{R}$ and similarity $I|V|=|\mathcal{R}|$.

For an isotropic distribution we have that only the inter-particle distance is
needed for their joint probability [30, 31]. That is
$\mathrm{p}(\bm{x}_{1},\bm{x}_{2})=\mathrm{p}(|\bm{x}_{1}-\bm{x}_{2}|)$, where
$|\bm{r}|$ is the magnitude of the vector $\bm{r}$. To use this property to
simplify the formula for the pair-correlation (4), we start by writing

$\mathrm{g}(|\bm{x}_{1}-\bm{x}_{2}|)=\mathrm{g}(\bm{x}_{1},\bm{x}_{2}),$ (9)

then integrate both sides over all values such that
$|\bm{x}_{1}-\bm{x}_{2}|=z$ for fixed $z$ relative to the variables of
integration. To make this clearer we introduce the ball region using standard
set-builder notation:

$\mathcal{B}(\bm{x};r)=\\{\bm{y}\in\mathbb{R}^{3}:\;|\bm{x}-\bm{y}|\leq r\\}.$
(10)

Using this notation, we will integrate both sides of (9) over all
$\bm{x}_{1}\in\mathcal{R}$, followed by over every
$\bm{x}_{2}\in\partial\mathcal{B}(\bm{x}_{1};z)$ such that
$\bm{x}_{2}\in\mathcal{R}$ which results in

$\int_{\mathcal{R}}\int_{\partial\mathcal{B}(\bm{x}_{1};z)\cap\mathcal{R}}\mathrm{g}(|\bm{x}_{1}-\bm{x}_{2}|)\mathrm{d}S_{2}\mathrm{d}\bm{x}_{1}=\int_{\mathcal{R}}\int_{\partial\mathcal{B}(\bm{x}_{1};z)\cap\mathcal{R}}\mathrm{g}(\bm{x}_{1},\bm{x}_{2})\mathrm{d}S_{2}\mathrm{d}\bm{x}_{1}.$
(11)

The left side of (11) can move out of the integral, resulting in

$\int_{\mathcal{R}}\int_{\partial\mathcal{B}(\bm{x}_{1};z)\cap\mathcal{R}}\mathrm{g}(z)\mathrm{d}S_{2}\mathrm{d}\bm{x}_{1}=L(z)\mathrm{g}(z),\;\;\text{with}\;\;L(z)=\int_{\mathcal{R}}\int_{\partial\mathcal{B}(\bm{x}_{1};z)\cap\mathcal{R}}\mathrm{d}S_{2}\mathrm{d}\bm{x}_{1}.$
(12)

For the right side of (11) we use the definition (4), and then rearrange (11)
to reach:

$\displaystyle\mathrm{g}(z)=\frac{J-1}{J}\frac{1}{L(z)}\int_{\mathcal{R}}\int_{\partial\mathcal{B}(\bm{x}_{1};z)\cap\mathcal{R}}\frac{\mathrm{p}(\bm{x}_{1},\bm{x}_{2})}{\mathrm{p}(\bm{x}_{1})\mathrm{p}(\bm{x}_{2})}\mathrm{d}\bm{x}_{1}\mathrm{d}S_{2}.$
(13)

To reach a simple formula to calculate $\mathrm{g}(z)$, we need to calculate
the integral $L(z)$. We show an example of this for a sphere in the Optional
box 2.2. In general, calculating $L(z)$ can be awkward and we show how to
avoid this in the next section.

1\. A sphere filled with spherical particles Using (12) we can calculate an
analytic expression for $L(z)$ be using the formula for the area of a
spherical cap [49], followed by integrating over $\bm{x}_{1}$ to obtain
$L(z)=\frac{\pi^{2}}{3}(z-2\bar{R})^{2}(z+4\bar{R}).$ (14) One way to check
this formula is to verify that
$\int_{0}^{2R}L(z)z^{2}\mathrm{d}z=\pi^{2}\bar{R}^{6}16/9=|\mathcal{R}|^{2}$.

To simplify the radial pair-correlation (13) we combine (7) and (8) to write
the pair-correlation (4) in terms of Dirac delta functions, which we then
substituted into (13), and integrate over $\bm{x}_{1}$ to obtain:

$\displaystyle\mathrm{g}(z)$
$\displaystyle=\frac{1}{S}\frac{1}{{\mathfrak{n}}^{2}}\sum_{s}^{S}\sum_{\bm{r}_{i}\in\mathcal{X}^{s}}\sum_{\bm{r}_{j}\not=\bm{r}_{i},\bm{r}_{j}\in\mathcal{X}^{s}}\frac{1}{L(z)}\int_{\partial\mathcal{B}(\bm{r}_{i};z)\cap\mathcal{R}}\delta(\bm{x}_{2}-\bm{r}_{j})\mathrm{d}S_{2}.$
(15)

To simplify this note that because $\bm{r}_{j}\in\mathcal{R}$ we have that

$\int_{\partial\mathcal{B}(\bm{r}_{i};z)\cap\mathcal{R}}\delta(\bm{x}_{2}-\bm{r}_{j})\mathrm{d}S_{2}=\int_{\partial\mathcal{B}(\bm{r}_{i};z)}\delta(\bm{x}_{2}-\bm{r}_{j})\mathrm{d}S_{2}.$

In more detail, this is a result of the integral over
$\bm{x}_{2}\not\in\mathcal{R}$ having an integrand which is zero,
$\delta(\bm{x}_{2}-\bm{r}_{j})=0$, because $\bm{r}_{j}\in\mathcal{R}$. The
above allows us to rewrite (15) in the form

$\displaystyle\mathrm{g}(z)$
$\displaystyle=\frac{1}{S}\frac{1}{{\mathfrak{n}}^{2}}\sum_{s}^{S}\sum_{\bm{r}_{i}\in\mathcal{X}^{s}}\sum_{\bm{r}_{j}\not=\bm{r}_{i},\bm{r}_{j}\in\mathcal{X}^{s}}\frac{1}{L(z)}\int_{\partial\mathcal{B}(\bm{0};z)}\delta(\bm{r}_{i}-\bm{r}_{j}+\bm{x})\mathrm{d}S_{x},$
(16)

where we used the change of variables $\bm{x}=\bm{x}_{2}-\bm{r}_{i}$. We can
now show that:

$\int_{\partial\mathcal{B}(\bm{0};z)}\delta(\bm{r}+\bm{x})\mathrm{d}S_{x}=\delta(|\bm{r}|-z),$
(17)

by noting that for any $z_{1}<z_{2}$ we have that

$\int_{z_{1}}^{z_{2}}\int_{\partial\mathcal{B}(\bm{0};z)}\delta(\bm{r}+\bm{x})\mathrm{d}S_{x}\mathrm{d}z=\int_{\mathcal{B}(\bm{0};z_{2})\setminus\mathcal{B}(\bm{0};z_{1})}\delta(\bm{r}+\bm{x})\mathrm{d}V_{x}=\begin{cases}1&\text{if}\;\;z_{1}<|\bm{r}|<z_{2},\\\
0&\text{else},\end{cases}$

where we note that $\mathrm{d}S_{x}\mathrm{d}z=\mathrm{d}V_{x}$ is a volume
element. Considering the properties of the Dirac delta [2], the above can be
used to deduce (17).

Substituting these results into (18) leads to

$\displaystyle\mathrm{g}(z)$
$\displaystyle=\frac{1}{S}\frac{1}{{\mathfrak{n}}^{2}L(z)}\sum_{s}^{S}\sum_{\bm{r}_{i}\in\mathcal{X}^{s}}\sum_{\bm{r}_{j}\not=\bm{r}_{i},\bm{r}_{j}\in\mathcal{X}^{s}}\delta(|\bm{r}_{i}-\bm{r}_{j}|-z).$
(18)

The only change in the above when changing from three spatial dimensions to
two spatial dimensions is that $L(z)$ is given by (12) but with $\mathcal{R}$
being two dimensional and the integral over $S_{2}$ being a line integral.

Figure 1: The left image shows a finite set of particles in a region
$\mathcal{R}_{1}$ taken from a larger set of disordered particles in the
region $\mathcal{R}_{2}$. Note that $\mathcal{R}_{1}$ is contained within
$\mathcal{R}_{2}$. On the right a unit cell of random particles in a region
$\mathcal{R}_{1}$ that is periodically tilled. On the right, the region
$\mathcal{R}_{2}$ is a cut out from the periodic tilling of the particles in
$\mathcal{R}_{1}$. We use $\ell$ to indicate the minimum length of
periodicity, which is the height of the unit cell shown in the image on the
right.

##  3 Particles in two regions

To avoid the influence of the boundary that encloses the particles, and avoid
calculating $L(z)$ appearing in (32), we consider two different regions where
particles can be placed $\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ as shown in
Figure 1. We consider that $\bm{x}_{1}$ and $\bm{x}_{2}$ are two different
random variables with $\bm{x}_{1}\in\mathcal{R}_{1}$ and
$\bm{x}_{2}\in\mathcal{R}_{2}$ and noting that
$\mathcal{R}_{1}\subset\mathcal{R}_{2}$.

To calculate the joint probability distribution
$\mathrm{p}(\bm{x}_{1},\bm{x}_{2})$, let $\mathcal{X}^{s}$ be such that every
$\bm{x}_{2}\in\mathcal{X}^{s}$ is also in $\bm{x}_{2}\in\mathcal{R}_{2}$, then
considering that $\bm{x}_{1}\in\mathcal{R}_{1}$ and
$\bm{x}_{2}\in\mathcal{R}_{2}$ we have that

$\mathrm{p}_{V}(\bm{v}_{i},\bm{v}_{j})=\frac{1}{S}\sum_{s}^{S}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{i})\right]}{\\#[\mathcal{X}^{s}\cap\mathcal{R}_{1}]}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{j})\right]}{\\#[\mathcal{X}^{s}\backslash
V(\bm{v}_{i})]}\frac{1}{|V(\bm{v}_{i})||V(\bm{v}_{j})|}\quad\text{for}\;\;i\not=j,$
(19)

such that $\bm{v}_{i}\in\mathcal{R}_{1}$, $\bm{v}_{j}\in\mathcal{R}_{2}$, and
$V(\bm{v}_{i})\cap V(\bm{v}_{j})=\varnothing$ for every $i,j$, and

$\cup_{i=1}^{J_{1}}V(\bm{v}_{i})=\mathcal{R}_{1}\quad\text{and}\quad\cup_{j=1}^{J_{2}}V(\bm{v}_{j})=\mathcal{R}_{2}.$

As a check of the formula (19), performing an integral of
$\mathrm{p}_{V}(\bm{v}_{i},\bm{v}_{j})$ for $\bm{v}_{i}\in\mathcal{R}_{1}$ and
$\bm{v}_{j}\in\mathcal{R}_{2}$ gives 1 as expected by following similar steps
as shown in (3).

For the two different regions, the definition of the particle pair-correlation
(4) now becomes

$\mathrm{g}(\bm{x}_{1},\bm{x}_{2}):=\frac{J_{2}-1}{J_{2}}\frac{\mathrm{p}(\bm{x}_{1},\bm{x}_{2})}{\mathrm{p}(\bm{x}_{1})\mathrm{p}(\bm{x}_{2})},$
(20)

where $J_{2}=\\#\mathcal{X}^{s}$.

###  3.1 Isotropic distributions

The probability density for the two regions, and an isotropic distribution of
particles, is similar to before with

$p(\bm{x}_{1})=\frac{1}{|\mathcal{R}_{1}|}\quad\text{and}\quad
p(\bm{x}_{2})=\frac{1}{|\mathcal{R}_{2}|}.$ (21)

Using analogous steps shown in Section 2.2 to reach (18) we can reach

$\mathrm{g}(z)=\frac{1}{{\mathfrak{n}}^{2}L(z)}\frac{1}{S}\sum_{s}^{S}\sum_{r_{i}\in\mathcal{X}^{s}\cap\mathcal{R}_{1}}\sum_{r_{j}\not=r_{i},r_{j}\in\mathcal{X}^{s}}\delta(|\bm{r}_{i}-\bm{r}_{j}|-z).$
(22)

However, now we can explicitly calculate $L(z)$ given by

$L(z)=\int_{\mathcal{R}_{1}}\int_{\partial\mathcal{B}(\bm{x}_{1};z)\cap\mathcal{R}_{2}}\mathrm{d}S_{2}\mathrm{d}\bm{x}_{1}.$
(23)

Assume we want to calculate $\mathrm{g}(z)$ for $0\leq z\leq Z$, to do this we
can require that: the distance between the boundaries
$\partial\mathcal{R}_{2}$ and $\partial\mathcal{R}_{1}$ be greater than or
equal to $Z$. Then $\bm{x}_{2}\in\partial\mathcal{B}(\bm{x}_{1};z)$ implies
that $\bm{x}_{2}\in\mathcal{R}_{2}$ for $z\leq Z$. In other words, the domain
of integration over $S_{2}$ becomes
$\partial\mathcal{B}(\bm{x}_{1};z)\cap\mathcal{R}_{2}=\partial\mathcal{B}(\bm{x}_{1};z)$,
which used in (23) for three spatial dimensions leads to

$L(z)=4\pi z^{2}|\mathcal{R}_{1}|.$ (24)

As a quick check, if $\mathcal{R}_{2}$ was a sphere with radius $R_{2}$
centred at the origin then
$\int_{0}^{R_{2}}L(z)\mathrm{d}z=|\mathcal{R}_{1}||\mathcal{R}_{2}|$ as it
should. Substituting (23) into (22) leads to

$\displaystyle\mathrm{g}_{3}(z)=\frac{1}{4\pi{\mathfrak{n}}z^{2}J_{1}}\frac{1}{S}\sum_{s}^{S}\sum_{r_{i}\in\mathcal{X}^{s}\cap\mathcal{R}_{1}}\sum_{r_{j}\not=r_{i},r_{j}\in\mathcal{X}^{s}}\delta(|\bm{r}_{i}-\bm{r}_{j}|-z).\quad\text{(3D
spatial)}$ (25)

where $J_{1}=\\#[\mathcal{X}^{s}\cap\mathcal{R}_{1}]$, and we use
$\mathrm{g}_{3}$ instead of $\mathrm{g}$ to indicate that the above is for 3
spatial dimensions.

For two spatial dimensions we have $L(z)=2\pi z|\mathcal{R}_{1}|$ which
substituted into (22) leads to

$\displaystyle\mathrm{g}_{2}(z)=\frac{1}{2\pi{\mathfrak{n}}zJ_{1}}\frac{1}{S}\sum_{s}^{S}\sum_{r_{i}\in\mathcal{X}^{s}\cap\mathcal{R}_{1}}\sum_{r_{j}\not=r_{i},r_{j}\in\mathcal{X}^{s}}\delta(|\bm{r}_{i}-\bm{r}_{j}|-z).\quad\text{(2D
spatial)}$ (26)

One way to verify the formulas above is given in Section 4.2.

Next, we re-deduce the results above without using Dirac delta functions. If
this is not appealing, then skip to Section 5 where we deduce the Structure
Factor.

##  4 The discrete form for isotropic pair-correlations

In this section we redo the calculations that lead to (25) and (26) but
without the use of the Dirac delta function. That is, without taking the limit
of volume elements tending to zero as used to reach (6). Doing this serves two
purpose: 1) it can be simpler to understand, and implement as a numerical
method, these discrete formulas and 2) it helps to verify the formulas by
reaching formulas we can compare with literature and having two avenues to
deduce the same formulas. Otherwise, this section can be skipped.

###  4.1 Particles in one region

Here we consider that all particles are within one region,
$\bm{r}_{i},\bm{r}_{j}\in\mathcal{X}^{s}$, and re-deduce the results that led
to the formula (18) but without Dirac deltas.

To reach a simple formula for $\mathrm{g}(z)$ with the discrete approximation
(2), we discretise the integral in (13), substitute (2), and use the discrete
approximation for the differentials

$\mathrm{d}\bm{x}_{1}=|V(\bm{v}_{i})|,\quad\mathrm{d}z\mathrm{d}S_{2}=|V(\bm{v}_{j})|,$
(27)

to obtain

$\mathrm{g}(z)=\frac{J-1}{J}\frac{|\mathcal{R}|^{2}}{L(z)\mathrm{d}zS}\sum_{s}^{S}\sum_{i}\sum_{|\bm{v}_{i}-\bm{v}_{j}|\approx
z}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{i})\right]}{\\#\mathcal{X}^{s}}\frac{\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{j})\right]}{\\#\left[\mathcal{X}^{s}\backslash
V(\bm{v}_{i})\right]},$ (28)

where we used (8) to substitute $p(\bm{x}_{1})=p(\bm{x}_{2})=1/|\mathcal{R}|$,
and the sum over $j$ is for every $\bm{v}_{j}$ such that
$z-\mathrm{d}z/2<|\bm{v}_{i}-\bm{v}_{j}|<z+\mathrm{d}z/2$. As the minimum
value for $z=2a$, we just set $\mathrm{g}(z)=0$ for $z\leq 2a$.

At this point, due to the choice (27), the volumes $|V(\bm{v}_{i})|$ and
$|V(\bm{v}_{j})|$ do not appear explicitly in (28). This allows us to simplify
the formula by choosing $V(\bm{v}_{i})$ to be small enough so that it contains
at most one particle, and likewise for $V(\bm{v}_{j})$. For the indices $i$
where there is no particle in $V(\bm{v}_{i})$ we have
$\\#\left[\mathcal{X}^{s}\cap V(\bm{v}_{i})\right]=0$ and likewise for $j$.
This makes it convenient now to only sum over the $i$, and $j$ where

$\displaystyle\\#\left[\mathcal{X}^{s}\cap
V(\bm{v}_{i})\right]=\\#\left[\mathcal{X}^{s}\cap V(\bm{v}_{j})\right]=1,$
(29) $\displaystyle\\#\left[\mathcal{X}^{s}\backslash
V(\bm{v}_{i})\right]=\\#\mathcal{X}^{s}-1=J-1,$ (30)

which used in (28) leads to

$\mathrm{g}(z)=\frac{1}{{\mathfrak{n}}^{2}}\frac{1}{L(z)dzS}\sum_{s}^{S}\sum_{i}\sum_{|\bm{v}_{i}-\bm{v}_{j}|\approx
z}1,$ (31)

where we used the particle number density ${\mathfrak{n}}:=J/|\mathcal{R}|$.
Now, when the volume elements $V(\bm{v}_{i})$ and $V(\bm{v}_{j})$ are close to
zero, then, as there is only one particle centre in $V(\bm{v}_{i})$ and one
other particle centre in $V(\bm{v}_{j})$, we will have that
$\bm{v}_{i}\approx\bm{r}_{i}$ and $\bm{v}_{j}\approx\bm{r}_{j}$ for some
$\bm{r}_{i},\bm{r}_{j}\in\mathcal{X}^{s}$, which means we can sum over
$\bm{r}_{i}$ and $\bm{r}_{j}$ instead of $\bm{v}_{i}$ and $\bm{v}_{j}$. This
leads us to rewrite (31) in the reduced form

$\mathrm{g}(z)=\frac{1}{{\mathfrak{n}}^{2}L(z)\mathrm{d}zS}\sum_{s}^{S}\sum_{\bm{r}_{i}\in\mathcal{X}^{s}}\\#\mathcal{X}^{s}_{i}(z),\quad\text{(discrete
pair-correlation)}$ (32)

where $\mathcal{X}^{s}_{i}(z)$ are all the particles $\bm{r}_{j}$ such that
are $|\bm{r}_{i}-\bm{r}_{j}|\approx z$, or more precisely with set builder
notation:

$\mathcal{X}^{s}_{i}(z):=\\{\bm{r}_{j}\in\mathcal{X}^{s}:z-\mathrm{d}z/2\leq|\bm{r}_{i}-\bm{r}_{j}|<z+\mathrm{d}z/2\\},$
(33)

and $\\#\mathcal{X}^{s}_{i}(z)$ is the number of elements in
$\mathcal{X}^{s}_{i}(z)$.

The only change in the formula (31) when changing from three spatial
dimensions to two spatial dimensions is that $L(z)$ is given by (12) but with
$\mathcal{R}$ being two dimensional and the integral over $S_{2}$ being a line
integral.

Below we show that (32) and (18) are equivalent. To do this we need to show
that

$\lim_{\mathrm{d}z\to
0}\frac{\\#\mathcal{X}^{s}_{i}(z)}{\mathrm{d}z}=\sum_{\bm{r}_{j}\in\mathcal{X}^{s}}\delta(|\bm{r}_{i}-\bm{r}_{j}|-z).$
(34)

To show the above let us define

$f(x)=\begin{cases}1&\text{if}\quad z-\mathrm{d}z/2\leq x\leq
z+\mathrm{d}z/2,\\\ 0&\text{otherwise}.\end{cases}$

Then
$\\#\mathcal{X}^{s}_{i}(z)=\sum_{\bm{r}_{j}\in\mathcal{X}^{s}}f(|\bm{r}_{i}-\bm{r}_{j}|)$.
To complete the demonstration of (34) we note that for any $x$ we have that

$\lim_{\mathrm{d}z\to 0}\frac{f(x)}{\mathrm{d}z}=\delta(x-z).$ (35)

###  4.2 Particles in two regions

Using (21) and the steps shown in Section 3, we can reach a formula which is
analogous to (32) but for $\bm{r}_{i}\in\mathcal{R}_{1}$ and
$\bm{r}_{j}\in\mathcal{R}_{2}$ as illustrated in Figure 1 given by

$\mathrm{g}_{3}(z)=\frac{1}{S}\frac{1}{4\pi
z^{2}dz}\sum_{s}^{S}\frac{1}{{\mathfrak{n}}J_{1}}\sum_{r_{i}\in\mathcal{X}^{s}\cap\mathcal{R}_{1}}\\#\mathcal{X}^{s}_{i}(z),\quad\text{(3D
isotropic)},$ (36)

where we substituted (24) and used
$J_{1}=\\#[\mathcal{X}^{s}\cap\mathcal{R}_{1}]$. The above is the same as the
formula [3, Equation (2)], and [20, Equation (8.3.8)] when specialising their
formula for particles with the same radius, and when taking
$\mathcal{R}_{1}=\mathcal{R}_{2}$.

For two spatial dimensions the pair-correlation becomes

$\mathrm{g}_{2}(z)=\frac{1}{S}\frac{1}{2\pi
zdz}\sum_{s}^{S}\frac{1}{{\mathfrak{n}}J_{1}}\sum_{r_{i}\in\mathcal{X}^{s}\cap\mathcal{R}_{1}}\\#\mathcal{X}^{s}_{i}(z),\quad\text{(2D
isotropic)}.$ (37)

As a final check, note that if particle positions were uncorrelated then
$\\#\mathcal{X}^{s}_{i}(z)\approx 4\pi{\mathfrak{n}}z^{2}\mathrm{d}z$, for
three spatial dimensions, because the number of particles would then be
proportional to the volume times the number density. Substituting this
approximation into (36) leads to $\mathrm{g}_{3}(z)=1$, which is what is
expected from uncorrelated particles.

##  5 The structure factor

The pair-correlation $\mathrm{g}(z)$ in the form (25) shows how
$\mathrm{g}(z)$ is a discontinuous function when calculated from a finite
number of particles. It is discontinuous in the variables $z$ and the position
of the particles $\bm{r}_{i}$. As $\mathrm{g}(z)$ is not a smooth function, we
can not use techniques from local nonlinear optimisation to calculate a
particle configuration to match a specific pair-correlation. To avoid this, we
can take a transform of the pair-correlation such as the structure factor:

$\displaystyle\mathrm{S}(k)=1+{\mathfrak{n}}\int(\mathrm{g}(r)-1)\mathrm{e}^{-\mathrm{i}\bm{k}\cdot\bm{r}}\mathrm{d}\bm{r}\quad\text{(The
Structure factor)},$ (38)

where $r=|\bm{r}|$ and $k=|\bm{k}|$. The above matches the typical definition
of the structure factor [43, 46, 47, 32], noting that the notation
$h(r)=\mathrm{g}(r)-1$ is commonly used.

For some of the calculations below we will perform the integral over
$\mathrm{g}(r)$ and $-1$ separately. To do this we note that

$\delta(\bm{k})=\frac{1}{(2\pi)^{n}}\int\mathrm{e}^{-\mathrm{i}\bm{k}\cdot\bm{r}}\mathrm{d}\bm{r},$
(39)

where $n$ is the spatial dimension, which in this paper is either $n=3$ or
$n=2$.

For any given isotropic pair-correlation $\mathrm{g}(r)$ the structure factor
(38) can be simplified to a 1D integral. To do this, we assume that particles
become uncorrelated at a distance of $R$ so that $\mathrm{g}(r)=1$ for $0\leq
r\leq R$. Then from (38) we can calculate***The best way to do the above
integral is to use a spherical coordinate system for $\bm{r}$ where the
“$z$-axis” of the coordinate system is aligned with $\bm{k}$, so that
$\bm{k}\cdot\bm{r}=kr\cos\phi$. that for three spatial dimensions:

$\mathrm{S}_{3}(k)=1+\frac{4\pi}{k}{\mathfrak{n}}\int_{0}^{R}(\mathrm{g}_{3}(r)-1)\sin(kr)r\mathrm{d}r,$
(40)

and for two spatial dimensions:

$\mathrm{S}_{2}(k)=1+2\pi{\mathfrak{n}}\int_{0}^{R}(\mathrm{g}_{2}(r)-1)\mathrm{J}_{0}(kr)r\mathrm{d}r.$
(41)

We will use the above to calculate the structure-factor from a target pair-
correlation in the next section.

We can further simplify the structure-factor when calculating it from a
configuration of particles. For three spatial dimensions, we can substitute
(25) into the structure-factor (40), but with $R=\infty$, and using the
property of the Dirac delta, to obtain:

$\displaystyle\mathrm{S}_{3}(k)=1+\frac{1}{S}\sum_{s}^{S}\frac{1}{J_{1}}\sum_{\bm{r}_{i}\in\mathcal{X}^{s}\cap\mathcal{R}_{1}}\sum_{\underset{\bm{r}_{j}\not=\bm{r}_{i}}{\bm{r}_{j}\in\mathcal{X}^{s}}}\frac{\sin(k|\bm{r}_{i}-\bm{r}_{j}|)}{k|\bm{r}_{i}-\bm{r}_{j}|},$
(3D Structure factor) (42)

for $k>0$. When calculating the above a term of the form
$-(2\pi)^{3}{\mathfrak{n}}\delta(\bm{k})$ appears, according to (39), however,
as we consider only $k>0$ it has no contribution to the above.

Following analogous steps, the two dimensional structure factor calculated
from (26) becomes:

$\displaystyle\mathrm{S}_{2}(k)=1+\frac{1}{SJ_{1}}\sum_{s}^{S}\sum_{\bm{r}_{i}\in\mathcal{X}^{s}\cap\mathcal{R}_{1}}\sum_{\underset{\bm{r}_{j}\not=\bm{r}_{i}}{\bm{r}_{j}\in\mathcal{X}^{s}}}\mathrm{J}_{0}(|k||\bm{r}_{i}-\bm{r}_{j}|),$
(2D Structure factor) (43)

for $k>0$, where $\mathrm{J}_{0}$ is the Bessel function of the first kind.

##  6 Particle configurations from the structure factor

In theory, the pair-correlation is calculated by taking into account an
infinite number of different particle configurations. Yet, many exotic
material properties can be achieved by choosing specific pair-correlations,
with one example being hyperuniform disordered materials [42, 43, 52]. So
being able to calculate one configuration of particles which closely
represents any given pair-correlation would provide a route to fabricate
particulate materials which exhibit exotic properties.

To recover a specific configuration of particles from a pair-correlation, we
show a method to find a configuration of particles which is the mean particle
configuration. Suppose we are given some pair-correlation
$\mathrm{g}^{\star}(z)$, then we want one configuration of particles that when
substituted into (36) will be close to $\mathrm{g}^{\star}(z)$, when removing
the sum over the ensembles by setting $S=1$.

###  6.1 Restrictions

For any given pair-correlation $\mathrm{g}$ and $\mathrm{S}$ structure factor
there are certain restrictions [50, 10, 46] that need to satisfied. These need
to be considered when choosing a target pair-correlation $\mathrm{g}^{\star}$
or structure-factor $\mathrm{S}^{\star}$. We will only consider a few of these
restrictions, and note that there may be an infinite number of necessary,
though more complicated, conditions on the pair-correlation [10].

The simplest restrictions are that

$\mathrm{g}^{\star}(r)\geq 0\quad\text{and}\quad\mathrm{S}^{\star}(k)\geq 0,$
(44)

where the first is due to the definition (20) together with probability
functions being positive. The second is discussed in [10, 39, 46, 50].

As we focus on disordered particulates we require that particles become
uncorrelated at some distance $R$. This requirement implies that

$\mathrm{g}^{\star}(r)=1\quad\text{for}\;\;r\geq R.$ (45)

Further, from the definition (20), together with (21), we have that

$\int_{\mathcal{R}_{2}}\int_{\mathcal{R}_{1}}\mathrm{g}(\bm{x}_{1},\bm{x}_{2})\mathrm{d}\bm{x}_{1}\mathrm{d}\bm{x}_{2}=\frac{J_{2}-1}{J_{2}}|\mathcal{R}_{1}||\mathcal{R}_{2}|\int_{\mathcal{R}_{1}}\int_{\mathcal{R}_{2}}p(\bm{x}_{1},\bm{x}_{2})\mathrm{d}\bm{x}_{1}\mathrm{d}\bm{x}_{2}\\\
=\frac{J_{2}-1}{J_{2}}|\mathcal{R}_{1}||\mathcal{R}_{2}|.$ (46)

Alternatively, using isotropy so that $\mathrm{g}$ depends on only
$|\bm{x}_{1}-\bm{x}_{2}|$, the condition (45), and specialising to three
spatial dimensions, we obtain

$\int_{\mathcal{R}_{2}}\int_{\mathcal{R}_{1}}(\mathrm{g}_{3}(|\bm{x}_{1}-\bm{x}_{2}|)-1)\mathrm{d}\bm{x}_{1}\mathrm{d}\bm{x}_{2}=\int\int_{0}^{R}\int_{\mathcal{R}_{1}}(\mathrm{g}_{3}(z)-1)\mathrm{d}\bm{x}_{1}z^{2}\mathrm{d}z\mathrm{d}\Omega=\\\
4\pi|\mathcal{R}_{1}|\int_{0}^{R}(\mathrm{g}_{3}(z)-1)z^{2}\mathrm{d}z,$ (47)

where we used the change of variables $\bm{z}=\bm{x}_{2}-\bm{x}_{1}$, and the
condition (45). Combining (47) and (46) leads to the restriction:

$\int_{0}^{R}\mathrm{g}_{3}(z)z^{2}\mathrm{d}z=\frac{R^{3}}{3}-\frac{1}{4\pi{\mathfrak{n}}},$
(48)

where we used ${\mathfrak{n}}=J_{2}/|\mathcal{R}_{2}|$. For two spatial
dimensions, following analogous steps, we obtain the restriction

$\int_{0}^{R}\mathrm{g}_{2}(z)z\mathrm{d}z=\frac{R^{2}}{2}-\frac{1}{2\pi{\mathfrak{n}}}.$
(49)

We can translate the pair-correlation to satisfy the above. Let
$\mathrm{g}_{2}(r)=\mathrm{g}_{2}^{0}(r)+\alpha dp(r)$ where $dp(r)\to 0$ when
$r\to R$. Then given a number density ${\mathfrak{n}}$ we can obtain $dp$ from
the above:

$\alpha=\Big{[}\frac{R^{2}}{2}-\frac{1}{2\pi{\mathfrak{n}}}-\int_{0}^{R}\mathrm{g}_{2}(z)z\mathrm{d}z\Big{]}\Big{[}\int_{0}^{R}dp(z)z\mathrm{d}z\Big{]}^{-1}.$
(50)

A possible choice is $dp(r)=\mathrm{e}^{-6r/R}.$

###  6.2 Gradient optimisation

For some inner product

$\langle F,G\rangle_{k}=\int F(k)G(k)w(k)\mathrm{d}k,$

where $w(k)$ is some known weight. The objective is to find a particle
configuration $\mathcal{X}$ that minimises

$\min_{\bm{r}_{i}\in\mathcal{X}}f(\mathcal{X}),\quad\text{where}\quad
f(\mathcal{X}):=\langle\mathrm{S}-\mathrm{S}^{\star},\mathrm{S}-\mathrm{S}^{\star}\rangle_{k}.$
(51)

Most methods in the literature [24, 36] achieve this by using non-gradient
based methods such as genetic algorithms, Nelder-Mead, and Simulated
Annealing. However, $f$, as written above, is a smooth function of the
position of the particles in $\mathcal{X}$. Further, it is straightforward to
analytically calculate the gradient of $f$ in terms of the particle positions
$\bm{r}_{j}$. This implies the gradient based methods [37] are likely to out
perform non-gradient based methods. In specific, the objective function $f$,
and its gradient, can be computational expensive to calculate so we opt to use
the Limited memory BFGS (L-BFGS) methods [23], as it stores information on the
Hessian (the gradient of the gradient) and uses this to accelerate
convergence. This often implies that L-BFGS requires less evaluations of the
objective function and its gradient [23].

Specifically, we develop a method to minimise (51) in two steps, one global,
and one local. Separating the two steps allows us to search over a large area
of the parameter space with the global step, while still obtaining high
precision with the local step. Another clear reason, based on recovering a
configuration of particles, is that if we had the constraint of particles not
overlapping, then it can lead to particles being locked in configurations
which can be far from the global minimum. For this reason we only enforce no
particle overlapping in the local step. In more detail the two steps are:

  1. 1.

The global step. A global optimisation that completely rearranges all
particles to minimise (51). For this step, we will allow particles to overlap,
to avoid locking the particles in a configuration, and we will use a limited
range for the wavenumbers $k_{1}\leq k\leq k_{2}$ when minimizing (51) . That
is, in this step we do not want to resolve spatial details smaller than the
length scale $2a$. So the shortest wavelength $\lambda$ we consider for the
structure factor is $\lambda=2a$ which corresponds to wavenumbers $k\leq
k_{2}=\pi/a$. The smallest wavenumber $k_{1}$ is determined by the dimensions
of the material: let $D_{1}$ be the smallest dimension of $\mathcal{R}_{1}$,
then the longest wavelength we consider is $\lambda=D_{1}$ which implies that
$k\geq k_{1}=2\pi/D_{1}$.

  2. 2.

The local step. In this step, we enforce a penaliser $W$, shown below, when
the particles overlap. For this step we also want to resolve spatial details.
Suppose we want to resolve details up to $a/4$, then $k_{2}=8\pi$ and
$k_{1}\leq k\leq k_{2}$ for this step.

To use techniques from nonlinear optimisation [37] to minimise calculate (51)
we need to calculate the gradient:

$\displaystyle\frac{\partial
f(\mathcal{X})}{\partial\bm{r}_{j}}=\frac{\partial}{\partial\bm{r}_{j}}\langle\mathrm{S}-\mathrm{S}^{\star},\mathrm{S}-\mathrm{S}^{\star}\rangle_{k}=2\sum_{i}\langle\frac{\partial\mathrm{S}}{\partial
R_{ij}},\mathrm{S}-\mathrm{S}^{\star}\rangle_{k}\frac{\partial
R_{ji}}{\partial\bm{r}_{j}}$ (52)
$\displaystyle=2\langle\sum_{i}\frac{\partial\mathrm{S}}{\partial
R_{ji}}\frac{\bm{R}_{ji}}{R_{ji}},\mathrm{S}-\mathrm{S}^{\star}\rangle_{k},$
(53)

where, $\bm{R}_{ji}=\bm{r}_{j}-\bm{r}_{i}$, $R_{ji}=|\bm{r}_{j}-\bm{r}_{i}|$,
and

$\frac{\partial R_{ji}}{\partial\bm{r}_{j}}=\frac{\bm{R}_{ji}}{R_{ji}}.$

In particular, for (43) we have

$\frac{\partial\mathrm{S}_{2}}{\partial
R_{ji}}=\frac{k}{J_{1}}\begin{cases}\mathrm{J}_{0}^{\prime}(|k|R_{ji}),\quad&\text{if}\quad\bm{r}_{j}\not\in\mathcal{R}_{1},\\\
2\mathrm{J}_{0}^{\prime}(|k|R_{ji}),\quad&\text{if}\quad\bm{r}_{j}\in\mathcal{R}_{1},\end{cases}$
(54)

where there are two cases because: if $\bm{r}_{j}\in\mathcal{R}_{1}$ then
$\mathrm{J}_{0}(kR_{ji})$ gets summed twice in (43), but if
$\bm{r}_{j}\not\in\mathcal{R}_{1}$ then $\mathrm{J}_{0}(kR_{ij})$ only appears
once in the sum. Likewise for (42) we have

$\frac{\partial\mathrm{S}_{3}}{\partial
R_{ji}}=\frac{k}{J_{1}}\frac{1}{(kR_{ij})^{2}}\begin{cases}kR_{ij}\cos(kR_{ij})-\sin(kR_{ij}),\quad&\text{if}\quad\bm{r}_{j}\not\in\mathcal{R}_{1},\\\
2kR_{ij}\cos(kR_{ij})-2\sin(kR_{ij}),\quad&\text{if}\quad\bm{r}_{j}\in\mathcal{R}_{1}.\end{cases}$
(55)

For most optimisation methods, we choose to use Optim.jl [27], we need to
supply the total gradient:

$\nabla f=\left[\frac{\partial
f(\mathcal{X})}{\partial\bm{r}_{1}},\frac{\partial
f(\mathcal{X})}{\partial\bm{r}_{2}},\ldots,\frac{\partial
f(\mathcal{X})}{\partial\bm{r}_{j}}\right],$

where the block vector on the right is typically flattened to be just one
large vector.

For the local step, after the global step is complete, we add restriction that
penalises particles that are overlapping. That is instead of minimizing (51)
we minimise:

$\min_{\bm{r}_{i}\in\mathcal{X}}f(\mathcal{X})+AW(\mathcal{X}),$ (56)

where $A$ is some large positive constant that is problem dependant and

$W=\sum_{i,j\not=i}\chi_{\\{R_{ji}<2a\\}}\mathrm{e}^{-4R_{ji}^{2}/(2a)^{2}},$
(57)

The above contributes a gradient of

$\frac{\partial
W}{\partial\bm{r}_{j}}=-\frac{4}{a^{2}}\sum_{i\not=j}\chi_{\\{R_{ji}<2a\\}}\mathrm{e}^{-4R_{ji}^{2}/(2a)^{2}}\bm{R}_{ji}.$
(58)

For ease of implementation we use a discrete form:

$\langle h,g\rangle_{k}=\sum_{q}h_{q}g_{q}w_{q},$

where $w_{q}$ are some Gaussian quadrature weights. For example we can write1

$\frac{\partial
f(\mathcal{X})}{\partial\bm{r}_{j}}=2\sum_{iq}\frac{\partial\mathrm{S}_{q}}{\partial
R_{ji}}\frac{\bm{R}_{ji}}{R_{ji}}(\mathrm{S}_{q}-\mathrm{S}^{\star}_{q})w_{q}$
(59)

###  6.3 Preliminary numerical results

Here we give some preliminary results, and mention what can be achieved in the
future with our method in the conclusion.

The first step is to have a systematic way to choose candidates pair-
correlations, with one motivation being to control wave propagation [15]. Any
chosen pair-correlation needs to satisfy the restrictions given in Section
6.1. For a preliminary result we take a well studied pair-correlation Percus-
Yevick [20, 43, 1] which represents disordered particles with very short range
correlation: they are uniformly distributed, and are only correlated because
they can not overlap. An example is shown in Figure 2

Figure 2: On the left, the Percus-Yevick pair-correlation for hard discs [1],
in two spatial dimensions, where the discs occupy 15% of the volume fraction.
On the right the corresponding structure factor when using (41).

The next step is to generate an initial configuration of particles with the
correct sizes, and within a given volume fraction. The simplest way to do this
is to place particles on a grid. This also facilitates defining the regions
$\mathcal{R}_{1}$ and $\mathcal{R}_{2}$ which are needed to calculate the
pair-correlation without introducing artefacts from the boundary, as shown in
Section 3. However, placing the particles exactly on a periodic grid would
lead to problems: due to symmetry, this configuration is likely on local
maximum, and therefore not an ideal starting position for gradient based
methods. To remedy this, we just move the particles by a small distance in a
random direction.

Figure 3: The left image shows the initial position of all the particles,
while the right image shows the particle configuration that best fits the
structure factor shown in Figure 2 according to our optimisation method.

The first step, the global step, of our method minimise (51) and is able to
exactly match the structure factor given, which is shown in Figure 2. The
result of the structure factor of the optimal particle configuration (after
the local step) is shown in Section 6.3. After applying the global step the
predicted structure-factor is almost the same (at least under the eye-ball
norm), while the predicted pair-correlation is also shown in Section 6.3 on
the right. The pair-correlation for the particulate is likely rather noisy, as
there are only 600 particles that occupy 15% of the volume fraction in two-
spatial dimensions. But clearly, as a preliminary result we can see our method
is able to produce configurations of particle that match a desired pair-
correlation and structure-factor.

##  7 Discussion

Summary. In this paper, we deduce from first principals how to calculate both
the pair-correlation and structure-factor of a finite disordered particulate.
We are able to demonstrate how to do this without getting artefacts from the
boundary of the particulate, something which seems to be ignored in the
literature. This means we can calculate the pair-correlation of an infinite
particulate from a finite sample. This is generally desired, as most
theoretical methods use the pair-correlations from an infinite medium,
including the very important translation invariance which greatly simplifies
the pair-correlation.

Being able to calculate the pair-correlation from a particulate is rather
straightforward in comparison to the inverse: calculating a particulate from a
pair-correlation. We present a method to calculate one configuration of
particles that best fits a given pair-correlation or structure-factor. This
problem is called the realizability problem, and has been primarily motivated
by using Neutron and X-ray scattering to measure molecular structures [9, 45,
19], as the structure factor can be measured from these experiments. For this
reason much of the literature focused on using assumptions that are valid for
molecular structures, and use molecular dynamics [5, 48] to help identify the
particle configuration. Alternatively, we look at obtaining one configuration
of particles from the structure factor with the minimal possible assumptions,
such as particles not allowing particles to overlap.

Most methods in the literature [46, 17] that calculate particle configurations
from the pair-correlation of structure factor use non-gradient based methods
such as genetic algorithms, Nelder-Mead, and Simulated Annealing. Several
works in fact focus on trying to approximate the pair-correlation with a
particulate, which is an inherently discontinuous function as the joint
probability function (7) involves Dirac deltas. Instead, we suggest that it is
best to seek a configuration of particles to approximate the structure factor,
as it is a smooth function of the particle positions, as shown by (42) and
(43). This enables us to analytically calculate the gradient of the structure
factor, in terms of the particle positions. With an analytic gradient it is
likely that gradient based methods [37] will out perform non-gradient based
methods.

Future directions. We present a two step method to calculate the structure
factor from a configuration of particles, one step for global optimisation
that avoids particle locking, and one step for local optimisation. This method
can be further developed to easily add priors about the particle configuration
and treat this as statistical inverse problem [18]. There is a significant
amount of prior information that could be used [29, 51]. For example, when
using more than one type of particle, some particles may repel or attract each
other. Or there can be specific knowledge on chains or sub-components of
particles. This information can be added as a prior, or regulariser. This
seems to be an unexplored approach that could greatly increase the performance
of this inverse problem.

## References

  * [1] M Adda-Bedia et al.  “Solution of the Percus-Yevick equation for hard disks”  In _The Journal of chemical physics_ 128.18 American Institute of Physics, 2008, pp. 184508
  * [2] George B Arfken et al.  “Mathematical methods for physicists: a comprehensive guide”  Academic press, 2011
  * [3] JA Barker et al.  “Monte Carlo values for the radial distribution function of a system of fluid hard spheres”  In _Molecular Physics_ 21.1 Taylor & Francis, 1971, pp. 187–191
  * [4] V.V. Bringi et al.  “The effects on pair correlation function of coherent wave attenuation in discrete random media”  In _IEEE Transactions on Antennas and Propagation_ 30.4 Institute of ElectricalElectronics Engineers (IEEE), 1982, pp. 805–808  DOI: 10.1109/tap.1982.1142852
  * [5] O. L. Caballero et al.  “Molecular Dynamics Simulations for Neutrino Scattering in Heterogeneous High Dense Media”  In _AIP Conference Proceedings_ AIP, 2008  DOI: 10.1063/1.2905137
  * [6] Rémi Carminati et al.  “Principles of Scattering and Transport of Light”  Cambridge University Press, 2021  DOI: 10.1017/9781316544693
  * [7] F. Carsughi et al.  “Small-angle neutron scattering from silica particles in solution with different concentrations”  In _Physica B: Condensed Matter_ 234-236 Elsevier BV, 1997, pp. 343–346  DOI: 10.1016/s0921-4526(97)00980-0
  * [8] F. S. Carvalho et al.  “Radial distribution function for liquid gallium from experimental structure factor: a Hopfield neural network approach”  In _Journal of Molecular Modeling_ 26.8 Springer ScienceBusiness Media LLC, 2020  DOI: 10.1007/s00894-020-04436-y
  * [9] L Cormier et al.  “A reverse Monte Carlo study of a titanosilicate glass”  In _Journal of Physics: Condensed Matter_ 9.46 IOP Publishing, 1997, pp. 10129–10136  DOI: 10.1088/0953-8984/9/46/011
  * [10] O. Costin et al.  “On the Construction of Particle Distributions with Specified Single and Pair Densities”  In _The Journal of Physical Chemistry B_ 108.51 American Chemical Society (ACS), 2004, pp. 19614–19618  DOI: 10.1021/jp047793m
  * [11] Jenness Crawford et al.  “Aspects of correlation function realizability”  In _The Journal of Chemical Physics_ 119.14 AIP Publishing, 2003, pp. 7065–7074  DOI: 10.1063/1.1606678
  * [12] M. Fábián et al.  “Network structure of multi-component sodium borosilicate glasses by neutron diffraction”  In _Journal of Non-Crystalline Solids_ 353.18-21 Elsevier BV, 2007, pp. 2084–2089  DOI: 10.1016/j.jnoncrysol.2007.02.030
  * [13] Hiroyuki Fujii et al.  “Structural properties of liquid Au–Si and Au–Ge alloys with deep eutectic region”  In _Journal of Non-Crystalline Solids_ 353.18-21 Elsevier BV, 2007, pp. 2094–2098  DOI: 10.1016/j.jnoncrysol.2007.02.031
  * [14] Artur L. Gower et al.  “Multiple waves propagate in random particulate materials”  In _SIAM J. Appl. Math._ 79.6 SIAM, 2019, pp. 2569–2592
  * [15] Artur Lewis Gower et al.  “Effective waves for random three-dimensional particulate materials”  In _New Journal of Physics_ IOP Publishing, 2021
  * [16] C. J. Horowitz et al.  “Neutrino pasta scattering: The opacity of nonuniform neutron-rich matter”  In _Physical Review C_ 69.4 American Physical Society (APS), 2004  DOI: 10.1103/physrevc.69.045804
  * [17] C. J. Horowitz et al.  “Nonuniform neutron-rich matter and coherent neutrino scattering”  In _Physical Review C_ 70.6 American Physical Society (APS), 2004  DOI: 10.1103/physrevc.70.065806
  * [18] Jari Kaipio et al.  “Statistical and computational inverse problems”  Springer Science & Business Media, 2006
  * [19] D. A. Keen et al.  “Structural modelling of glasses using reverse Monte Carlo simulation”  In _Nature_ 344.6265 Springer ScienceBusiness Media LLC, 1990, pp. 423–425  DOI: 10.1038/344423a0
  * [20] Jin Au Kong et al.  “Scattering of electromagnetic waves: numerical simulations”  John Wiley & Sons, 2004
  * [21] T. Kuna et al.  “Realizability of Point Processes”  In _Journal of Statistical Physics_ 129.3 Springer ScienceBusiness Media LLC, 2007, pp. 417–439  DOI: 10.1007/s10955-007-9393-y
  * [22] Tobias Kuna et al.  “Necessary and sufficient conditions for realizability of point processes”  In _The Annals of Applied Probability_ 21.4 Institute of Mathematical Statistics, 2011  DOI: 10.1214/10-aap703
  * [23] Dong C Liu et al.  “On the limited memory BFGS method for large scale optimization”  In _Mathematical programming_ 45.1-3 Springer, 1989, pp. 503–528
  * [24] C. Manwart et al.  “Reconstruction of random media using Monte Carlo methods”  In _Physical Review E_ 59.5 American Physical Society (APS), 1999, pp. 5596–5599  DOI: 10.1103/physreve.59.5596
  * [25] R L McGreevy  “Reverse Monte Carlo modelling”  In _Journal of Physics: Condensed Matter_ 13.46 IOP Publishing, 2001, pp. R877–R913  DOI: 10.1088/0953-8984/13/46/201
  * [26] R. L. McGreevy et al.  “Reverse Monte Carlo Simulation: A New Technique for the Determination of Disordered Structures”  In _Molecular Simulation_ 1.6 Informa UK Limited, 1988, pp. 359–367  DOI: 10.1080/08927028808080958
  * [27] Patrick K. Mogensen et al.  “Optim: A mathematical optimization package for Julia”  In _Journal of Open Source Software_ 3.24 The Open Journal, 2018, pp. 615  DOI: 10.21105/joss.00615
  * [28] Rana Nandi et al.  “Transport Properties of the Nuclear Pasta Phase with Quantum Molecular Dynamics”  In _The Astrophysical Journal_ 852.2 American Astronomical Society, 2018, pp. 135  DOI: 10.3847/1538-4357/aa9f12
  * [29] Edoardo Patelli et al.  “On optimization techniques to reconstruct microstructures of random heterogeneous media”  In _Computational Materials Science_ 45.2 Elsevier BV, 2009, pp. 536–549  DOI: 10.1016/j.commatsci.2008.11.019
  * [30] R.K. Pathria et al.  “Statistical Mechanics of Interacting Systems: The Method of Cluster Expansions”  In _Statistical Mechanics_ Elsevier, 2011, pp. 299–343  DOI: 10.1016/b978-0-12-382188-1.00010-4
  * [31] Oliver H. E. Philcox et al.  “Disordered Heterogeneous Universe: Galaxy Distribution and Clustering across Length Scales”  In _Physical Review X_ 13.1 American Physical Society (APS), 2023  DOI: 10.1103/physrevx.13.011038
  * [32] J. Piekarewicz et al.  “Proton fraction in the inner neutron-star crust”  In _Physical Review C_ 85.1 American Physical Society (APS), 2012  DOI: 10.1103/physrevc.85.015807
  * [33] Szilvia Pothoczki et al.  “Partial radial distribution functions of methylene halide molecular liquids”  In _Journal of Molecular Liquids_ 153.2-3 Elsevier BV, 2010, pp. 112–116  DOI: 10.1016/j.molliq.2010.01.011
  * [34] L. Pusztai et al.  “MCGR: An inverse method for deriving the pair correlation function from the structure factor”  In _Physica B: Condensed Matter_ 234-236 Elsevier BV, 1997, pp. 357–358  DOI: 10.1016/s0921-4526(96)00986-6
  * [35] L. Pusztai et al.  “Reverse Monte Carlo model calculations on a-C:H two-component systems”  In _Zeitschrift für Physik B Condensed Matter_ 101.4 Springer ScienceBusiness Media LLC, 1996, pp. 631–636  DOI: 10.1007/s002570050256
  * [36] Jacques A. Quiblier  “A new three-dimensional modeling technique for studying porous media”  In _Journal of Colloid and Interface Science_ 98.1 Elsevier BV, 1984, pp. 84–102  DOI: 10.1016/0021-9797(84)90481-8
  * [37] Andrzej Ruszczynski  “Nonlinear optimization”  Princeton university press, 2011
  * [38] Bo Sjöberg et al.  “Interparticle interactions and structure in nonideal solutions of human serum albumin studied by small-angle neutron scattering and Monte Carlo simulation”  In _Biophysical Chemistry_ 52.2 Elsevier BV, 1994, pp. 131–138  DOI: 10.1016/0301-4622(94)00089-1
  * [39] Frank H. Stillinger et al.  “Pair Correlation Function Realizability: Lattice Model Implications”  In _The Journal of Physical Chemistry B_ 108.51 American Chemical Society (ACS), 2004, pp. 19589–19594  DOI: 10.1021/jp0478155
  * [40] S. Torquato et al.  “Controlling the Short-Range Order and Packing Densities of Many-Particle Systems”  In _The Journal of Physical Chemistry B_ 106.33 American Chemical Society (ACS), 2002, pp. 8354–8359  DOI: 10.1021/jp0208687
  * [41] S. Torquato et al.  “Exactly solvable disordered sphere-packing model in arbitrary-dimensional Euclidean spaces”  In _Physical Review E_ 73.3 American Physical Society (APS), 2006  DOI: 10.1103/physreve.73.031106
  * [42] Salvatore Torquato  “Hyperuniform states of matter”  In _Physics Reports_ 745 Elsevier BV, 2018, pp. 1–95  DOI: 10.1016/j.physrep.2018.03.001
  * [43] Salvatore Torquato et al.  “Random heterogeneous materials: microstructure and macroscopic properties”  In _Appl. Mech. Rev._ 55.4, 2002, pp. B62–B63
  * [44] G. Torrie et al.  “Monte Carlo calculation of y(r) for the hard-sphere fluid”  In _Molecular Physics_ 34.6 Informa UK Limited, 1977, pp. 1623–1628  DOI: 10.1080/00268977700102821
  * [45] Matthew G Tucker et al.  “RMCProfile: reverse Monte Carlo for polycrystalline materials”  In _Journal of Physics: Condensed Matter_ 19.33 IOP Publishing, 2007, pp. 335218  DOI: 10.1088/0953-8984/19/33/335218
  * [46] O.U. Uche et al.  “On the realizability of pair correlation functions”  In _Physica A: Statistical Mechanics and its Applications_ 360.1 Elsevier BV, 2006, pp. 21–36  DOI: 10.1016/j.physa.2005.03.058
  * [47] Kevin Vynck et al.  “Light in correlated disordered media”  In _arXiv preprint arXiv:2106.13892_ , 2021
  * [48] G. A. Wijs et al.  “First-principles molecular-dynamics simulation of liquid Li12Si7”  In _Physical Review B_ 48.18 American Physical Society (APS), 1993, pp. 13459–13468  URL: https://doi.org/10.1103/physrevb.48.13459
  * [49] Wikipedia contributors  “Spherical cap — Wikipedia, The Free Encyclopedia” [Online; accessed 12-January-2022], https://en.wikipedia.org/w/index.php?title=Spherical_cap&oldid=1062812955, 2021
  * [50] Masami Yamada  “Geometrical Study of the Pair Distribution Function in the Many-Body Problem”  In _Progress of Theoretical Physics_ 25.4 Oxford University Press (OUP), 1961, pp. 579–594  DOI: 10.1143/ptp.25.579
  * [51] C. L. Y. Yeong et al.  “Reconstructing random media”  In _Physical Review E_ 57.1 American Physical Society (APS), 1998, pp. 495–506  DOI: 10.1103/physreve.57.495
  * [52] Ge Zhang et al.  “Realizable hyperuniform and nonhyperuniform particle configurations with targeted spectral functions via effective pair interactions”  In _Physical Review E_ 101.3 American Physical Society (APS), 2020  DOI: 10.1103/physreve.101.032124


# SkillRec: A Data-Driven Approach to Job Skill Recommendation for Career
Insights

Xiang Qian Ong  Information Systems Technology and Design Pillar
Singapore University of Technology and Design
Singapore
<EMAIL_ADDRESS>Kwan Hui Lim  Information Systems Technology
and Design Pillar
Singapore University of Technology and Design
Singapore
<EMAIL_ADDRESS>
###### Abstract

Understanding the skill sets and knowledge required for any career is of
utmost importance, but it is increasingly challenging in today’s dynamic world
with rapid changes in terms of the tools and techniques used. Thus, it is
especially important to be able to accurately identify the required skill sets
for any job for better career insights and development. In this paper, we
propose and develop the Skill Recommendation (SkillRec) system for
recommending the relevant job skills required for a given job based on the job
title. SkillRec collects and identify the skill set required for a job based
on the job descriptions published by companies hiring for these roles. In
addition to the data collection and pre-processing capabilities, SkillRec also
utilises word/sentence embedding techniques for job title representation,
alongside a feed-forward neural network for job skill recommendation based on
the job title representation. Based on our preliminary experiments on a
dataset of 6,000 job titles and descriptions, SkillRec shows a promising
performance in terms of accuracy and F1-score.

###### Index Terms:

Recommendation Systems, Skill Recommendation, Occupational Analysis, Natural
Language Processing, Neural Networks

##  I Introduction

Technological advancements and industrial changes are constantly evolving at a
rapid pace [1]. With this rapid change, new job scopes are constantly being
generated at companies along with a dynamic list of skill sets to meet these
ever-growing changes. As such, individuals that are looking to change their
career or seeking a new job often face the problem of skill mismatch due to
those constant changes. This problem is further highlighted by the recent
COVID-19 pandemic which has caused drastic changes in many industries, such as
aviation and tourism. During this period, there was a spike in the number of
retrenched workers seeking to find employment in other fields [2] and
potentially facing skill mismatch.

On top of that, jobs are tightly coupled with skill sets, which are essential
for fulfilling the job requirements. Individuals who are transiting into a new
career and industry would need to conduct extensive research or seek advice on
the latest job requirements and the relevant skills. This process is commonly
performed to improve the chances of applicants securing a new job, especially
in a new industry. This is necessary due to the huge barrier of entry for
individuals transiting over to new jobs or new industries, given their
existing skill sets. Different industry fields and even roles within the same
industry require different sets of skills to carry out tasks associated with
the job. Likewise, hiring managers also look out for those skill sets when
assessing suitable candidates for an open role. As such, identifying the skill
sets required for a given job is especially important for potential job
seekers, even for those staying in the same industry.

###  I-A Related Work

In terms of academic research, this problem and research area have also
garnered keen interest in recent years due to these various developments [3,
1, 4, 5, 6, 7]. For example, there are various streams of research that
investigate diverse but related topics such as: (i) recommending job
candidates to companies based on their suitability [8, 9, 10]; (ii)
recommending jobs to candidates based on the candidates skill sets and career
history [11, 12]; (iii) modelling the career trajectories of individuals [13,
14]; (iv) predicting the future career transitions and movements of workers
[15, 16].

Among these works, the most closely related works to our research are those on
job skill recommendations [17, 18, 19] and those using various embedding
techniques [4, 20, 21, 22] for occupational analysis and mining. For example,
Dave et al. proposed a representation learning model using job transitions,
job-skills and skill co-occurence networks for recommending both jobs and
skills [17]. Others have used reinforcement learning in the form of a Deep-Q
Network for skill recommendation and better interpretation [19]. Liu et al.
developed a bi-directional language model for contextualised job tile
embeddings that was demonstrated on a downstream occupational named entity
recognition task [4]. Skills2Graph [21] used FastText embeddings to better
maintain semantic similarity between taxonomic components of job postings.
Also closely related to our research are various works that have developed
systems and implemented approaches for related occupational analysis and
mining tasks [23, 20, 24].

###  I-B Motivation and Problem Statement

The fundamental problem is that identifying the relevant skill sets required
for a particular job title is both time-consuming and a complex process.
However, this process is essential to job seekers looking to enter new jobs
and yet very challenging as they are new to both the field and industry. The
problem can be broken down into two sub-problems: the lack of knowledge to
learn the relevant skills and the lack of proficiency in the relevant skills,
which will be elaborated next.

####  I-B1 Lack of knowledge to learn the relevant skills

To learn the relevant skills required for the particular job, there is a need
to understand the current job workflow, the tools that are used on the job and
knowledge in that field of interest. The root of the problem is the lack of
domain expertise in the new industry to truly understand the skill sets
required for the job [25].

There are various different ways in which the problem can be mitigated. The
first is to seek advice from experts who have worked in the field for a long
time. They understand the job requirements and the necessary skill sets
needed, from their extensive years of experience. Alternatively, governments’
skills recommendations serve as a good indicator of possible general relevant
skill sets that are required in the industry for the near future [26, 27].
Lastly, individuals could perform online research on the industry, given that
there is an abundance of articles and forums that share insights and
recommendations.

This step is crucial as it serves as the basis for the next step. Only after
knowing the skills that need to be learned, then the individual can embark on
the next step.

####  I-B2 Lack of proficiency in the relevant skills

After identifying the required skill sets to transit to a new industry or
role, it is crucial for individuals to have access to resources to learn those
skills. Proficiency in these relevant skills is an essential requirement that
leads to successful employment and subsequent progression. The resources could
come in many different forms, both physical and digital, which we discuss
next.

A classroom setting is a common structure whereby individuals will attend
classes and undergo a structured curriculum to gain knowledge and experience
from the instructors. At the end of such courses, there are formal
certificates that act as proof of completion of the course. Classroom formats
could also be virtual, ranging from online courses such as Udemy and Coursera.
Similarly, these classes provide both the content and practical lessons, with
the added bonus of online certificates upon completion. Openly available
resources are unstructured curricula that could comprise different resources
from both online and offline. This gives the individual the flexibility of
time and learning opportunity but typically do not result in a formal
certification upon completion.

###  I-C Problem Definition

In this paper, our main focus is on the first sub-problem, which is the lack
of knowledge to learn the relevant skills required for a job role. This
problem is significantly more crucial than the second sub-problem, as the
former’s outcome will determine the skill sets that an individual is required
to learn. Failing at this step will render useless any future efforts at both
a industry transition or career progression. In addition, there are multiple
open-source and affordable resources online to address the problem of a lack
of proficiency, thus this second sub-problem is thus not as crucial as
compared to the first sub-problem that we focus on. More formally, we
investigate the problem of identifying and recommending the required job skill
sets required for a role given the job title.

###  I-D Main Contributions

The current ways to address the problems are mainly human-driven insights,
such as seeking the advice of current practitioners in the new industry or
through self-reading of openly available resources. However, we could tap into
data-driven insights that perform the analysis and recommend the right skill
sets to pick up for the industry or role we are interested in. This could
reduce the time taken or the cost incurred in the process.

In this paper, we make the following contributions:

  1. 1.

We investigate the problem of job skill recommendation, which is to recommend
a list of required skill sets given role in the form of a job title.

  2. 2.

To address the above problem, we propose and develop the Skill Recommendation
(SkillRec) system. SkillRec comprises the capability to collect openly
available job listing information, process these data to extract job title and
required skill sets, and subsequently recommend skill sets given a job title.

  3. 3.

As part of SkillRec, we also develop two approaches to skill recommendation,
utilising word/sentence embedding techniques (using BERT and FastText) for job
title representation, before inputing this job title representation to a
simple but effective feed-forward neural network for job skill recommendation.

  4. 4.

Using a dataset of 6,000 job listings, we perform preliminary experiments to
evaluate the performance of our proposed SkillRec system. Our preliminary
results show promising performance of SkillRec in terms of accuracy and
F1-score.

Figure 1: System Architecture of Our Proposed Skill Recommendation (SkillRec)
System

##  II Proposed Skill Recommendation System

Our proposed Skill Recommendation (SkillRec) system comprises four main
components, including a web collection component, data pre-processing
component, job title representation component and job skill recommendation
component. Figure 1 shows an overview of our proposed system. Next, we provide
a brief summary of each component, followed by a more detail description in
the later sections.

  1. 1.

Web Collection Component. This component retrieves job listings from various
web sources, such as online job posting platforms.

  2. 2.

Data Pre-processing Component. This component performs the necessary data pre-
processing and cleaning to extract relevant information from the job
descriptions, such as the job title, job description, required skills, company
name, etc.

  3. 3.

Job Title Representation Component. This component utilises word/sentence
embedding techniques, such as BERT and FastText, to model job titles in term
of its vector representations.

  4. 4.

Job Skill Recommendation Component. The job title representation from the
previous step is then fed into this component, which comprises a feed-forward
neural network that outputs a set of recommended skills given this job title.

###  II-A Web Collection Component

The web collection component retrieves various types of web-based information
such as job descriptions from popular online job posting platforms, as well as
skill sets from massive open online courses platforms.

The job description data was retrieved over a period of six months at a weekly
interval to ensure that the dataset is kept up-to-date. A total of 6,000 job
descriptions were retrieved from three popular online job posting platform.

Massive open online courses platforms provide a good source of information on
popular educational courses and the skill sets they aim to impart. Using two
popular massive open online courses platforms, we collected and identified 589
unique skill sets that are relevant for this work.

###  II-B Data Pre-processing Component

The data pre-processing component works on the job title and job description
data collected by the web collection component. This component first parses
the collected webpage to extract the relevant information, i.e., the job title
and job description. Thereafter, the standard text pre-processing steps, such
as conversion to lowercase, removal of non-textual symbols, tokenisation, are
performed on these two fields. The skill sets for a job listing is derived
from its corresponding job description based on an exact match of a mentioned
skill with our earlier list of 589 unique skill sets.

###  II-C Job Title Representation Component

As our main task is to determine the required skill sets for a given job
title, we proceed to model the job titles using various embedding techniques.
In this study, we experiment with two popular word/sentence representation
techniques:

  * •

BERT Embedding [28]. The Bidirectional Encoder Representation of Transformer
(BERT) is a transformer-based language model comprising multiple layers of
bidirectional Transformer encoders, along with self-attention mechanism. BERT
is then pre-trained with the Masked Language Model task of predicting hidden
tokens using other non-hidden tokens, and/or Next Sentence Prediction task of
predicting the next sentence given an input sentence.

  * •

FastText Embedding [29, 30]. FastText builds upon the skip-gram approach that
is popularly used by the Word2Vec model [31, 32]. Using FastText, each word is
represented by a bag-of-character n-grams, and in turn each character n-grams
is represented by a vector. The vector representation of each word is then
derived from the sum of these character n-grams vector representations.

BERT was chosen due to its strong performance in various language modelling
tasks, such as text classification, question answering, among others. In
addition, BERT has been frequently used in a variety of other recommendation
and prediction tasks, such as location prediction [33, 34] and tourism
recommendation [35, 36].

After extracting the job title in the previous step, either BERT or FastText
is then used to obtain the vector representation of the job title. This job
title representation is subsequently used as input to our neural network for
job skill recommendation, which we describe next.

###  II-D Job Skill Recommendation Component

Using the vector representation of the job title from the previous step, this
job title representation is then fed into a feed-forward neural network for
job skill recommendation based on the job title. We utilise a feed-forward
neural network with two fully-connected hidden layer, with the first hidden
layer comprising 1280 nodes, the second hidden layer comprising 640 nodes and
a final output layer of 589 nodes corresponding to the 589 unique skill sets.
We employ Sigmoid as our activation function, AdamW as our optimizer and a
dropout of 0.5 to reduce overfitting.

##  III Experiments and Results

Using our collected dataset, we use 80% for training our model as described in
Section II and use the remaining 20% for evaluation. For evaluation purposes,
we use the standard metrics of accuracy and F1-score, based on the recommended
job skill against the ground truth job skill in the original job description.

We experimented with the following algorithms and baselines:

  * •

Bag-of-Words. A classical bag-of-words model that is traditionally used for
classification tasks involving text-based data. In our application, we utilise
a frequency-based approach to model the frequency of which a particular skill
is mentioned in relation to a specific job. This model is used as a simple
baseline to compare our two proposed approaches.

  * •

FastText+NN. As described earlier in our SkillRec system, this model uses the
FastText embedding for job title representation (Section II-C), followed by a
neural network (NN) for job skill recommendation (Section II-D).

  * •

BERT+NN. Similar to the earlier FastText+NN model, this model uses the BERT
embedding for job title representation (Section II-C), followed by the same
neural network (NN) for job skill recommendation (Section II-D).

TABLE I: Experimental Results Algorithm | Accuracy | F1-score
---|---|---
Bag-of-Words | 0.4768 | 0.2887
FastText+NN | 0.9728 | 0.4931
BERT+NN | 0.9870 | 0.4973

Table I shows the average results of our experiments for the three models in
terms of accuracy and F1-score. The results show that our proposed SkillRec
system using BERT+NN is the best performing model, with FastText+NN being a
extremely close second. BERT has been shown to perform well for a range of
natural language processing tasks, from text classification to question
answering. Similarly for this task of occupational skill recommendation, BERT
as well as FastText perform well with the add-on of a straightforward neural
network. Unsurprisingly, Bag-of-Words performs the worst by a large margin due
to the simple job title representation approach that is unable to capture
similar job titles obtained.

TABLE II: Example Job Title and Relevant Skills Job Title | Related Jobs and Skills
---|---
developer | java, python, fullstack, android, software
software | stack, java, python, fullstack, android
chef | pastry, kitchen, hair, chemist, accountant
senior director | director, partner, pacific, planner, strategy
software engineer | software, engineer, python, java, embedded

In addition to the experiments measuring accuracy and F1-score of the job
skill recommendations, we also performed a qualitative check on the job title
representations learnt and the similar job titles and skill sets based on
these representations. Table II shows a sample of some job titles and the
related job titles and skill sets. This simple analysis shows the
effectiveness of our approach as occupations such as developer and software
engineer results in relevant skills like Java, Python, fullstack, embedded,
etc. Similarly, non-technical roles such as senior director also resulted in
similar roles like partner and planner and relevant skills like strategy. For
certain roles like chef, the results are more mixed with relevant terms like
pastry and kitchen being mentioned but also non-relevant ones like chemist and
accountant. This issue highlights a limitation due to our preliminary
experiments and limited dataset, which contains more technical roles than non-
technical ones.

##  IV Conclusion and Future Work

In this paper, we proposed the SkillRec system to address the problem of job
skill recommendation given a job title. SkillRec comprises four main
components for web collection, data pre-processing, job title representation
and job skill recommendation. In particular, our preliminary experiments show
that an approach that combines job title representation using BERT, along with
job skill recommendation using a straightforward feed-forward neural network
works effectively in terms of accuracy and F1-score.

Future enhancements to SkillRec can involve a more dynamic recommendation of
job skills to account for new skill sets that may emerge in future or existing
skill sets that evolve over time. Similarly, we can also include more context
in addition to the job titles, such as the industry, country, job applicant
details, etc, for making more contextualised and personalised job skill
recommendations.

##  V Acknowledgment

This research is funded in part by the Singapore University of Technology and
Design under grant RS-MEFAI-00005-R0201.

## References

  * [1] I. Khaouja, I. Kassou, and M. Ghogho, “A survey on skill identification from online job ads,” _IEEE Access_ , vol. 9, pp. 118 134–118 153, 2021.
  * [2] W. E. Forum, “The future of jobs report 2020,” _World Economic Forum_ , 2020\.
  * [3] A. Gugnani and H. Misra, “Implicit skills extraction using document embedding and its use in job recommendation,” in _Proceedings of the AAAI Conference on Artificial Intelligence_ , vol. 34, no. 08, 2020, pp. 13 286–13 293.
  * [4] J. Liu, Y. C. Ng, Z. Gui, T. Singhal, L. Blessing, K. L. Wood, and K. H. Lim, “Title2vec: a contextual job title embedding for occupational named entity recognition and other applications,” _Journal of Big Data_ , vol. 9, no. 1, pp. 1–16, 2022.
  * [5] Z. Ao, G. Horváth, C. Sheng, Y. Song, and Y. Sun, “Skill requirements in job advertisements: A comparison of skill-categorization methods based on wage regressions,” _Information Processing & Management_, vol. 60, no. 2, p. 103185, 2023.
  * [6] A. Barducci, S. Iannaccone, V. La Gatta, V. Moscato, G. Sperlì, and S. Zavota, “An end-to-end framework for information extraction from italian resumes,” _Expert Systems with Applications_ , vol. 210, p. 118487, 2022\.
  * [7] J. Liu, Y. C. Ng, K. L. Wood, and K. H. Lim, “Ipod: a large-scale industrial and professional occupation dataset,” in _Conference Companion Publication of the 2020 on Computer Supported Cooperative Work and Social Computing_ , 2020, pp. 323–328.
  * [8] C. Zhu, H. Zhu, H. Xiong, C. Ma, F. Xie, P. Ding, and P. Li, “Person-job fit: Adapting the right talent for the right job with joint representation learning,” _ACM Transactions on Management Information Systems (TMIS)_ , vol. 9, no. 3, p. 12, 2018.
  * [9] D. Shen, H. Zhu, C. Zhu, T. Xu, C. Ma, and H. Xiong, “A joint learning approach to intelligent job interview assessment.” in _IJCAI_ , 2018, pp. 3542–3548.
  * [10] C. Qin, H. Zhu, T. Xu, C. Zhu, L. Jiang, E. Chen, and H. Xiong, “Enhancing person-job fit for talent recruitment: An ability-aware neural network approach,” in _The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval_, 2018, pp. 25–34.
  * [11] J. Malinowski, T. Keim, O. Wendt, and T. Weitzel, “Matching people and jobs: A bilateral recommendation approach,” in _Proceedings of the 39th Annual Hawaii International Conference on System Sciences (HICSS’06)_ , vol. 6, 2006, pp. 137c–137c.
  * [12] Y. Zhang, C. Yang, and Z. Niu, “A research of job recommendation system based on collaborative filtering,” in _2014 Seventh International Symposium on Computational Intelligence and Design_ , vol. 1, 2014, pp. 533–538.
  * [13] Y. Liu, L. Zhang, L. Nie, Y. Yan, and D. S. Rosenblum, “Fortune teller: predicting your career path,” in _Thirtieth AAAI conference on artificial intelligence_ , 2016.
  * [14] D. Mimno and A. McCallum, “Modeling career path trajectories,” 2008.
  * [15] C. James, L. Pappalardo, A. Sîrbu, and F. Simini, “Prediction of next career moves from scientific profiles,” _arXiv preprint arXiv:1802.04830_ , 2018.
  * [16] Y. Yang, D.-C. Zhan, and Y. Jiang, “Which one will be next? an analysis of talent demission,” _Proceedings of the 1st International Workshop on Organizational Behavior and Talent Analytics_ , 2018.
  * [17] V. S. Dave, B. Zhang, M. Al Hasan, K. AlJadda, and M. Korayem, “A combined representation learning approach for better job and skill recommendation,” in _Proceedings of the 27th ACM International Conference on Information and Knowledge Management_ , 2018, pp. 1997–2005.
  * [18] A. Maurya and R. Telang, “Bayesian multi-view models for member-job matching and personalized skill recommendations,” in _2017 IEEE International Conference on Big Data (Big Data)_. IEEE, 2017, pp. 1193–1202.
  * [19] Y. Sun, F. Zhuang, H. Zhu, Q. He, and H. Xiong, “Cost-effective and interpretable job skill recommendation with deep reinforcement learning,” in _Proceedings of the Web Conference 2021_ , 2021, pp. 3827–3838.
  * [20] A. Giabelli, L. Malandri, F. Mercorio, M. Mezzanzanica, and A. Seveso, “Skills2graph: Processing million job ads to face the job skill mismatch problem.” in _Proceedings of the Thirtieth International Joint Conference on Artificial Intelligence_ , 2021.
  * [21] ——, “Skills2job: A recommender system that encodes job offer embeddings on graph databases,” _Applied Soft Computing_ , vol. 101, p. 107049, 2021.
  * [22] C. M. Jaramillo, P. Squires, H. G. Kaufman, A. M. Da Silva, and J. Togelius, “Word embedding for job market spatial representation: tracking changes and predicting skills demand,” in _2020 IEEE International Conference on Big Data (Big Data)_. IEEE, 2020, pp. 5713–5715.
  * [23] X. J. S. Ashok, E.-P. Lim, and P. K. Prasetyo, “Jobsense: A data-driven career knowledge exploration framework and system,” in _2018 IEEE International Conference on Data Mining Workshops (ICDMW)_. IEEE, 2018, pp. 1411–1416.
  * [24] C. Qin, H. Zhu, C. Zhu, T. Xu, F. Zhuang, C. Ma, J. Zhang, and H. Xiong, “Duerquiz: A personalized question recommender system for intelligent job interview,” in _Proceedings of the 25th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining_, 2019, pp. 2165–2173.
  * [25] P. Restrepo, “Skill mismatch and structural unemployment,” _MIT Working Paper_ , 2015.
  * [26] A. Felstead, D. Gallie, F. Green, and Y. Zhou, _Skills at Work in Britain, 1986 to 2006_. ESRC Centre on Skills, Knowledge and Organisational Performance, 2007.
  * [27] A. Felstead, D. Gallie, F. Green, and H. Inanc, “Skills at work in britain: First findings from the skills and employment survey 2012,” 2013.
  * [28] J. Devlin, M.-W. Chang, K. Lee, and K. Toutanova, “Bert: Pre-training of deep bidirectional transformers for language understanding,” in _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , 2019, pp. 4171–4186.
  * [29] P. Bojanowski, E. Grave, A. Joulin, and T. Mikolov, “Enriching word vectors with subword information,” _Transactions of the association for computational linguistics_ , vol. 5, pp. 135–146, 2017.
  * [30] A. Joulin, É. Grave, P. Bojanowski, and T. Mikolov, “Bag of tricks for efficient text classification,” in _Proceedings of the 15th Conference of the European Chapter of the Association for Computational Linguistics: Volume 2, Short Papers_ , 2017, pp. 427–431.
  * [31] T. Mikolov, K. Chen, G. Corrado, and J. Dean, “Efficient estimation of word representations in vector space,” _arXiv preprint arXiv:1301.3781_ , 2013\.
  * [32] T. Mikolov, I. Sutskever, K. Chen, G. S. Corrado, and J. Dean, “Distributed representations of words and phrases and their compositionality,” _Advances in neural information processing systems_ , vol. 26, 2013.
  * [33] M. Li, K. H. Lim, T. Guo, and J. Liu, “A transformer-based framework for poi-level social post geolocation,” _Proceedings of the 45th European Conference on Information Retrieval_ , 2023.
  * [34] L. F. Simanjuntak, R. Mahendra, and E. Yulianti, “We know you are living in bali: Location prediction of twitter users using bert language model,” _Big Data and Cognitive Computing_ , vol. 6, no. 3, p. 77, 2022.
  * [35] J. Arreola, L. Garcia, J. Ramos-Zavaleta, and A. Rodriguez, “An embeddings based recommendation system for mexican tourism,” in _IberLEF@SEPLN_ , 2021, pp. 110–117.
  * [36] N. L. Ho and K. H. Lim, “Poibert: A transformer-based model for the tour recommendation problem,” _Proceedings of the 2022 IEEE International Conference on Big Data_ , 2022.