# Fast dynamics and high effective dimensionality of liquid fluidity

C. Cockrell1,∗, O. Dicks2, I. T. Todorov3, A. M. Elena3 and K. Trachenko4 1
Department of Materials, Imperial College London, Exhibition Road, London, SW7
2AZ, UK  2 Department of Physics and Astronomy, University of British
Columbia, Agricultural Road, Vancouver, V6T 1Z1, Canada  3 Daresbury
Laboratory, Keckwick Ln, Daresbury, Warrington, WA4 4AD, UK  4 School of
Physical and Chemical Sciences, Queen Mary University of London, Mile End
Road, London, E1 4NS, UK

###### Abstract

Fluidity, the ability of liquids to flow, is the key property distinguishing
liquids from solids. This fluidity is set by the mobile transit atoms moving
from one quasi-equilibrium point to the next. The nature of this transit
motion is unknown. Here, we show that flow-enabling transits form a
dynamically distinct sub-ensemble where atoms move on average faster than the
overall system, with a manifestly non-Maxwellian velocity distribution. This
is in contrast to solids and gases where no distinction of different ensembles
can be made and where the distribution is always Maxwellian. The non-
Maxwellian distribution is described by an exponent $\alpha$ corresponding to
high dimensionality of space. This is generally similar to extra synthetic
dimensions in topological quantum matter, albeit higher dimensionality in
liquids is not integer but is fractional. The dimensionality is close to 4 at
melting and exceeds 4 at high temperature. $\alpha$ has a maximum as a
function of temperature and pressure in liquid and supercritical states,
returning to its Maxwell value in the solid and gas states.

###### pacs:

65.20.De 65.20.JK 61.20Gy 61.20Ja

##  I Introduction

Similarly to gases and differently from solids, liquids flow. However,
important liquid properties such as density, bulk modulus and heat capacity,
especially at conditions not far from melting, are very close to those of
solids but are very different from those of gases. This dichotomy has
contributed to the problem of theoretical description of liquids and the long-
held belief that a general liquid theory is unworkable landaustat  and not
tractable due to the absence of a small parameter Pitaevskii ; akhiezer ; ropp
. This is markedly different from solids and gases where a general theory,
universally describing their properties, was developed long ago.

The dynamics of individual liquid particles also shows a dichotomy related to
solid and gas dynamics first introduced by Frenkel frenkel : the motion of
liquid particles combines small solidlike oscillations around quasi-
equilibrium positions inside a cage formed by particle neighbours and larger-
distance jumps to new positions. Observed in molecular dynamics simulations
and referred to as “transits” wallacemd , these jumps endow liquids with their
ability to flow dyre . The distance travelled by transits is close to the
inter-particle separation as envisaged by Frenkel frenkel ; wallacemd .

In order to understand liquids at the level similar to that existing in solids
and gases, we need to understand microscopic dynamics of liquid particles.
This dynamics includes individual particle velocities and sets the system
trajectory in phase space and all of system properties landaustat . This
dynamics also underpins properties averaged over different particles and time
such as correlation functions often discussed in the area of liquids hansen1 ;
balucani . Useful in some cases, averaging in liquids, whose dynamics combine
oscillation and transits, inevitably lose detailed information about these
distinct components of motion contained in the primary phase trajectory data.

Conceptually, an oscillatory motion is well-understood and is at the
cornerstone of physics. While it took time to ascertain this motion in
liquids, accumulating experiments show that the oscillatory motion in liquids
is very similar to that in solids. Strong evidence for this includes close
similarity of phonon dispersion curves measured in liquids and solids copley ;
pilgrim ; burkel ; pilgrim2 ; water-tran ; hoso ; hoso3 ; monaco1 ; monaco2 ;
sn ; ropp . This similarity extends to the zone boundary as in solids,
corresponding to particle vibrating inside its cage as envisaged by Frenkel.

Differently from the oscillatory component of motion, the microscopic nature
and properties of diffusive transits in liquids remain largely unknown and are
unclear even at the conceptual level. This includes the transit energy,
distribution of this energy across different transits and so on. The problem
is that we can’t use results for flowing diffusive transits from the dense gas
theory chapman : differently from gases, liquids are condensed state of
matter. This implies that transits in liquids constantly move in a strong
fluctuating field created by other particles, with the inter-particle
interaction being as strong as in solids. For this reason, useful concepts
such as sudden collisions and related mean free path used in the theory of
gases do not apply to dense strongly-interacting liquids chapman .

Here, we find that transits in liquids have several unusual and unexpected
properties. First, they move on average faster than the the whole system.
Second, the distribution of transit velocities is non-Maxwellian. Third, this
distribution is characterised by an exponent $\alpha$ corresponding to
effective space dimensionality larger than $3$. This is generally similar to
extra synthetic dimensions in topological quantum matter naturereview ;
phystoday , albeit extra dimensionality in liquids is not integer but is
fractional. Interestingly, $\alpha$ has a maximum as a function of temperature
and pressure, reducing to the normal Maxwell value in solid and gas states.
This effect is unrelated to dynamic heterogeneity discussed in the slow
supercooled regime of liquids approaching glass transition.

##  II Methods

We perform molecular dynamics (MD) simulations using the DL$\\_$POLY package
dlpoly . The interactions in our system are set by the widely-used Lennard-
Jones potential, with parameters $\epsilon=0.003177$ eV and $\sigma=2.775$ Å
corresponding to neon. We equilibrate systems with periodic boundaries along
isobars and isotherms in the NPT ensemble, and use the average density of
these runs to generate configurations for use in production runs in the NVE
ensemble. We select the NVE ensemble because we are interested in the
individual dynamics of atoms, and therefore cannot allow thermostats or
barostats interfering with the natural Hamiltonian evolution of the system. In
all simulations a timestep of 1 fs was employed. Systems with 500 to 32000
atoms were used in this study. Results discussed below are size-independent.

The main property we will be concerned with is the distribution of particle
velocities of transit atoms, $\rho(v)$, where $v=|{\bf v}|$ and ${\bf v}$ is
velocity. We start with describing an algorithm to detect transits. As
mentioned earlier, the main property of a transit is its large displacement
(compared to a smaller oscillatory motion inside the cage): the transit moves
from one quasi-equilibrium position to another, endowing the liquid with the
ability to flow. This implies that transits are to be detected based on the
distance travelled, $d$: $d>d_{c}$, where $d_{c}$ is the distance cutoff.

With this definition, we calculate $\rho(v)$ from the molecular dynamics
trajectory, which we print every 10 fs. At the moment an atom’s position
satisfies the transit criterion, we log its instantaneous velocity. Because we
work with fluids which lack a crystal structure, the reference structure for a
transit at the start is simply the initial condition of the simulation. The 10
fs printing frequency is a compromise between accuracy and computational load.
Frequent printing allows the capture of a jump as it happens, rather than
after its conclusion. Test runs with more frequent printing - every timestep
(1 fs) - did not yield any difference in the results which follow. We then
construct a histogram from this list of transit velocities with a bin width of
1 m/s, which is normalised to create a probability distribution $\rho(v)$.
Once an atom has been flagged as in-transit, its starting position could be
reset within the reference structure in order to collect its velocity again
for detection of subsequent transits. In practice, we improve statistics by
limiting simulation time to 30 ps and running between 200 and 4000 different
initial conditions (using different velocity seeds), such that the final
distributions are composed of, typically, 500000 to 2000000 velocities,
depending on the desired degree of noise. This gives the same distribution as
a single long simulation where atomic reference positions are reset after a
transit.

We note that there is no time averaging involved in detecting transits: we do
not add a condition of the particle moving a certain distance within a certain
timeframe. Instead, we log the instantaneous velocity at the moment the atom
is designated as a transit. This is essential in order to study the dynamics
of mobile transits.

We show the radial distribution function (RDF) at 100 bar and two different
temperature in Fig. 1a. The low temperature $T=25$ K marks several points,
including the first maxima and minima. We use these distances as useful
reference points for cutoffs $d_{c}$. $\rho(v)$ of transits calculated for
$d_{c}$ in the range 3-10 Å are shown in Fig. 1b (using larger $d_{c}$ is
limited by the simulation time needed to acquire long diffusion distances,
especially at low temperature). We observe that the distributions are very
close to each other. This was tested for other conditions explored in this
study. For sufficiently large $d_{c}$, the algorithm catches an atom
undergoing its 1st, 2nd, 3rd and other subsequent transit motions. These
consecutive transits are physically equivalent to each other. Hence it doesn’t
matter if we record an atom’s velocity during its first or any subsequent
transits. This explains the insensitivity of $\rho(v)$ to $d_{c}$. We note
that the algorithm could catch atoms in the middle of oscillatory motion
between the transits, however a large-amplitude transit is much more likely to
cross a large distance as compared to a small-amplitude oscillating atom. The
convergence of $\rho(v)$, starting from sufficiently large $d_{c}$, into a
single function indicates that these distributions represent a physically
meaningful sub-system of atoms.

Differently from large $d_{c}$, $\rho(v)$ deviate at small $d_{c}$, as
displayed in Fig. 1c, and start to visibly converge towards the Maxwell
distribution of the whole system below about 1.0 Å and for $d_{c}$ closer to
amplitude of the the oscillatory motion. This is easy to understand: short
$d_{c}$ include oscillatory motion, and therefore the motion of nearly any
atom will satisfy the transit criterion. This cutoff at which the distribution
of transit atoms becomes cutoff-independent, 1.0 Å, corresponds to the
approximate maximum amplitude of oscillatory motion. An atom which meets this
displacement is therefore undergoing a non-oscillatory transit. This
interpretation is ratified by the plot of the mean square displacement (MSD)
as a function of time in Fig. 1d. The nonlinear short-time behaviour is well-
known to correspond to the free motion of a particle before it interacts with
and is accelerated by its neighbours. This nonlinear regime therefore roughly
corresponds to the amplitude of oscillation, and we fit a quadratic function
to it. The long-time behaviour is the combination of oscillation and transits
regime, and to it we fit a linear function. The departure from the quadratic
ballistic regime visibly begins at around 0.6 Å, and concludes shortly after
1.0 Å, which coincides with the convergence of $\rho(v)$ to its high-$d_{c}$
value at around 0.8-1.0 Å. The distribution $\rho(v)$ therefore becomes
constant when the MSD switches from its ballistic to its diffusive regime.

Figure 1: (a) RDF at 100 bar in the liquidlike state at 25 K and the gaslike
state at 120 K; (b) Distribution of transiting atoms with cutoffs
corresponding to the maxima and minima of the RDF as well as a large cutoff of
10 Å; (c) Distribution of transiting atoms with small cutoffs, showing
convergence onto the Maxwell distribution of the whole system.

This implies that $\rho(v)$ are insensitive to the definitions of transits for
sufficiently large $d_{c}$ and can be consistently used in our analysis. We
note a related way to characterise transit motion used earlier. These
characterisations involved the calculation of the liquid relaxation time
$\tau$, the time it takes the transit to move from one quasi-equilibrium point
to the next frenkel . Tuning $d_{c}$ gives $\tau$ similar to those based on
other methods such as decay of correlation functions sastry . Previous work
egami  used a topological definition of $\tau$ as time between an atom gaining
or losing a nearest neighbour. We calculated velocity distributions of atoms
which satisfied this criterion, again using a variety of $d_{c}$. We found
that such a distribution is unfailingly identical to the Maxwell distribution
because the histogram $\rho(v)$ now includes oscillating atoms whose motion
contributes to an atom gaining or losing a neighbour. We combined this
topological definition with our geometric definition: an atom was defined as
in-transit if it underwent a neighbour change and had been displaced beyond a
certain $d_{c}$. The displacement cutoff and the neighbour cutoff were defined
to be equal to one another. This definition yielded the same distributions as
those resulting from the displacement criterion alone, at all pressure and
temperature conditions.

Recall that our selection criterion of transits is not based on particle’s
velocity and therefore is not related to selecting faster atoms directly. The
resulting $\rho(v)$ converge into a single cutoff-independent function, which
is not the case for the criteria based on velocity. In Fig. 2a, we plot the
velocity distribution formed by logging the velocities of the fastest subsets
of atoms in the system at each timestep. These distributions depend heavily on
the subset size: as this size becomes large, the distributions become the tail
of the Maxwell distribution, as they should. This is demonstrated in Fig. 2b
where $\rho(v)$ from Fig. 2a are re-scaled to the tail of the Maxwell
distribution. Fig. 2b also includes a distribution formed from atoms moving
faster than 300 m/s. We therefore find that velocity-based criteria do not
converge on a single distribution, whereas the physically justified criterion
based on large distance covered by transits do.

Figure 2: (a) Probability density of the fastest 0.02%, 4% and 20% atoms in
the system; (b) Distributions of fastest subsets rescaled such that their
convergence onto the tail of the Maxwell distribution is evident, as well as
the distribution formed by explicitly selecting atoms faster than 300 m/s.

##  III Results and Discussion

We now analyse the velocity distributions of transits. The central observation
from Fig. 1b is that $\rho(v)$ of transits is shifted to the right relative to
the Maxwell distribution, implying that transits move on average faster than
atoms in the overall system. We use nonlinear fitting in Mathematica with a
generalised Maxwell distribution as a fitting function:

$\rho(v)\propto v^{\alpha}\exp(-\frac{mv^{2}}{2k_{\rm B}T})$ (1)

When $\alpha=2$, Eq. (1) is the standard Maxwell distribution in three
dimensions landaustat . If $\rho(v)$ is calculated over all atoms in the
system, such as seen in Fig. 1 with a cutoff of 0 Å, this gives the fitted
temperature very close to the simulation temperature and the Maxwell
distribution of velocities with $\alpha=2$.

As seen in Fig. 3, $\rho(v)$ of transits can not be described by Eq. (1) with
$\alpha=2$: the fit is very poor. The fit markedly improves if we let $\alpha$
vary. This includes both cases where we fix $T$ to the simulation temperature
and where we let $T$ vary. In the latter case, the fit quality is slightly
better (this is not visible in Fig. 3) and returns temperature within 4%
agreement with the simulation temperature. At conditions in Fig. 3, $\alpha$
is about 3.25.

Figure 3: Generalised Maxwell distributions from Eqs. (1) and (3) fit to the
velocity distribution of transiting atoms at 100 bar and 25 K.

According to Fig. 1b, transits move on average faster than atoms in the
overall system. The maximum of distribution (1) corresponds to velocity
$v_{m}$

$v_{m}^{2}=\alpha\frac{k_{\rm B}T}{m}$ (2)

Larger $v_{m}$ correspond to larger $\alpha$, consistent with our earlier
funding of $\alpha=$3.25.

Other modified Maxwell distributions were also employed in order to fit to the
transit velocity distributions. To potentially describe the shift of
velocities of faster transits to the right of the distribution, we fixed
$\alpha$ and $T$, and offset the distribution by a fitting parameter:

$\rho(v)\propto\Theta(v-\mu)\
(v-\mu)^{2}\exp\left(-\frac{m(v-\mu)^{2}}{2k_{\mathrm{B}}T}\right)$ (3)

where $\mu$ is the offset parameter and the Heaviside step function $\Theta$
deletes the unphysical behaviour of the function outside the domain $v<\mu$.
This resulted in a better fit, however the low-velocity tail of the
distribution still does not match that of the transit velocity distribution.

The above analysis indicates two important points. First, $\rho(v)$ of
transits are not the same shape as that of a standard Maxwell distribution.
This may not be unexpected in itself, in view that the detection algorithm
imposes a condition on atoms to qualify as transits. What is notable, however,
is that $\rho(v)$ are remarkably close to a generalised Maxwell distribution
with $\alpha>2$.

Second, we find that transits are on average faster than average motion of the
system and therefore represent a distinct dynamical sub-ensemble of the
overall system. This is seen in Fig. 1b and is consistent with $\alpha>2$. In
this picture, a liquid is a mixture of two distinct dynamical ensembles with
different velocity distributions, albeit these ensembles are not fixed but
dynamically change: an atom performs an oscillatory motion for time $\tau$
($\tau$ is liquid relaxation time frenkel ; dyre ), then becomes a transit
atom, moving to a new position where it oscillates again and so on. These two
distinct ensembles correspond to faster transit atoms and slower oscillating
atoms, respectively. We note that this distinction between two sub-ensembles
is lost if an average property over the entire system is calculated as is
usually the case hansen1 ; balucani .

We characterise the nature of a transit by examining the speed of atoms in a
time window around the moment of transit. This is done by collecting velocity
data from atoms in times preceding and following the transit as per our above
definition. At long times before and after the transit, the mean speed is
expected to return to the mean speed defined by the Maxwell distribution, as
transiting is expected to happen to all atoms with equal probability, such
that individual atoms are not more prone to transiting than other atoms. We
plot the speed profile averaged over 200,000 jumping atoms at 100 bar and 20
K, just above the melting line, in Fig. 4. This profile extends 2.4 ps before
and after the moment of transit, defined by the above condition, which defines
$\Delta t=0$. This is done for several transit cutoffs, demonstrating the same
cutoff-independence seen in Fig. 1. We observe that the mean speed is peaked
just before the moment of transit, implying that atoms accelerate out of their
local environment into a new environment, facilitated by the opening up of a
fast diffusion pathway. The subsequent minimum and secondary maximum after
this peak represents the deceleration as the atom settles into its new
oscillatory quasi-equilibrium position. At long times, the mean speed becomes
identical to the mean speed of the Maxwell distribution, seen in the flat line
of the 0 Å cutoff. Small but nonzero cutoffs deviate from this velocity and
converge on the invariant speed profile at larger cutoffs, as was seen in
distributions in Fig. 1. We note that this speed profile is not time-
symmetric: the minimum and secondary maximum are not seen before the transit
time, and the peak is not centred on $\Delta t=0$. We also note that the peak
in Fig. 4 is seen for small as well as large cutofffs. Recall that the transit
is identified by its displacements from the starting position (in the
transiting atom’s past, which is the source of the asymmetry). Hence both fast
and slower atoms are counted as transits for short cutoffs, resulting in
smaller average peak velocity as compared to larger peak related to long
cutoffs.

Figure 4: Mean speed profile of transiting atoms for a window before and after
the moment of transit as a function of transit cutoff definition.

These separate ensembles are a contrast to the motion in solids and gases. In
solids, all motion is oscillatory and no separate regime can be defined. As a
result, the velocity distribution is Maxwellian. In gases where collisions are
rare and effectively instantaneous, atoms spend most of their time moving
freely with velocities similarly distributed according to the Maxwell
distribution.

We add three remarks regarding $\rho(v)$ of transits. First, there is no
apriori clear reason to assume that transits in liquids are on average faster.
One general consideration might be as follows. No jumps can take place in any
reasonable time in a liquid at fixed cage volume because this would involve
large energy barriers and long waiting times. Instead, Frenkel argued frenkel
, a volume increase of a surrounding fluctuating cage provides a fast
diffusion pathway for an atom to escape and become a transit. Traversal of
these diffusion pathways is associated with a higher kinetic energy due to
activation energy barriers involved. If we assume the speed of an oscillating
atom takes the form $v=v_{0}\sin(\omega t)$, then the average speed is
$\frac{2}{\pi}v_{0}$ over the oscillation period. If the transit speed, on the
other hand, is comparable to the amplitude, $v_{0}$ transits are expected to
be on average faster than oscillating atoms by a factor of $\frac{\pi}{2}$.
Whereas this is clearly an approximate argument, the direct quantitative
evidence for faster transits comes from our current MD simulations.

Second, the non-Maxwellian $\rho(v)$ of transits does not contradict the
Maxwell distribution in the entire classical system originating from the
factorisation of distributions dependent on momenta and coordinates,
respectively landaustat . As mentioned earlier, $\rho(v)$ of transits
corresponds to a dynamically changing sub-ensemble of the whole system which
has the Maxwell distribution overall.

Third, the ratio of transits identified by the algorithm to the overall number
of atoms is small: for example, at 40 K and 1 kbar, just above the melting
line, about 1% of atoms within a given 0.1 ps window, approximately
corresponding to a typical Debye vibration period, are in-transit, meaning
that over the course of an oscillation an atom has a 1% probability to be a
transit in the detection algorithm. In the gaslike state at 200 K and 1 kbar,
this increases to 8%. For this reason, the distribution of atoms tagged as
non-transits is very close to the velocity distribution in the overall system
which is Maxwellian.

We now address the variation of $\alpha$ with temperature and pressure. In
Fig. 5, we plot $\alpha$ as a function of temperature varied in a wide range
including both the subcritical liquid and the supercritical fluid state. We
observe an interesting non-monotonic behavior with a maximum. At low
temperature where the system is a solid, our algorithm does not detect
transits, and therefore the effective dimensionality is undefined. The
calculated distribution of velocities over all atoms in the solid phase is
Maxwellian with $\alpha=2$ as expected, which is what is shown below the
melting line in Fig. 5 for reference.

Figure 5: Parameter $\alpha$ in the generalised Maxwell distribution as a
function of: (a) temperature; (b) pressure on the 1 kbar isobar and 240 K
isotherm respectively.

At high temperature in the gaslike state, $\alpha$ tends to $2$ as well. This
takes place on both isobars and isotherms. At these extreme conditions above
the Frenkel line flreview , an atom lacks any quasi-equilibrium environment,
corresponding to the complete loss of all minima in the RDF, as seen in Fig.
1, with only a very broad and low first peak remaining. In this high-
temperature/low-pressure gaslike regime, particle dynamics is governed by the
gaslike ballistic motion where collisions set the particle mean free path
chapman . Velocity distribution is Maxwellian in this case, causing the return
of $\alpha$ to 2.

The value of $\alpha$ in Eq. (1) is related to space dimensionality: the
Maxwell distribution of velocity components $v_{x}$, $v_{y}$ and $v_{z}$ in
three-dimensional space corresponds to $\rho(v)\propto v^{2}$ in spherical
coordinates landaustat . A Maxwell distribution in two spatial dimensions has
the exponent $\alpha=1$, and in four dimensions $\alpha=3$. Therefore,
$\alpha>2$ corresponds to the motion in effectively higher dimensionality.

This higher effective dimensionality is related to the additional condition
imposed upon a transit. This condition is associated with the higher energy
involved in the particle escape from the surrounding cage: the velocity
distribution of transits is deformed towards faster $v$ as is seen in Fig. 1b.
We have seen that this deformation can not be represented by a higher
temperature, but only by a higher space dimensionality. $\alpha=2$ in the
Maxwell distribution in Eq. (1) corresponds to a geometric factor $v^{2}$,
arising from the two-dimensional surface corresponding to a single speed in a
uniform three-dimensional velocity space. Deformed towards faster $v$, the
speed distribution of transits can also be represented in a uniform velocity
space, albeit with higher dimensionality.

As is seen in Fig. 5, $\alpha$ is generally not an integer but is fractional.
The largest $\alpha$ is close to $4$ ($\alpha=3.7$) and corresponds to an
effective space dimensionality close to five. In view that $\alpha$ is not
generally an integer, it is interesting to note that close and above the
melting line along the isobar and isotherm in Fig. 5, $\alpha$ is close to 3,
such that the subsystem of transits follows a Maxwell distribution in a four-
dimensional space. $\alpha=3$ at melting could be a potentially useful insight
for a theory of the solid-liquid melting transition, in view that mobile
transiting defects aid the formation of the liquid nucleus and lower the
barrier for melting tuckerman .

In a general sense, the existence of two dynamical sub-ensembles might
resemble dynamic heterogeneity (DH) noted in slow glassy dynamics of viscous
supercooled liquids approaching glass transition ediger . This DH can
originate from the expected correlation of mobility and structure: higher
mobility is caused by fewer neighbours and wider surrounding cage due to the
nature of the model used, which was heterogeneous with two different atomic
species hetero1 ; hetero2 . This is a model-specific and hence a non-generic
effect. The effect reported here is unrelated to DH in glassy dynamics and
operates in a markedly different part of the phase diagram. For example, the
maximum of $\alpha$, corresponding to fast dynamics where the deviation from
the Maxwell distribution is the largest, takes place well above the critical
point, namely at $6.8T_{c}$ in Fig. 5a and $4.5P_{c}$ in Fig. 5b, where
$T_{c}$ and $P_{c}$ are critical temperature and pressure, 45 K and 27 bar
respectively in this system. This effective dimensionality is therefore a
universal effect, operating both in subcritical liquids and deep in the
supercritical state, rather than solely in supercooled liquids. Differently
from effective higher dimensionality reported here, DH in supercooled liquids
due to structural restraints can be anisotropic and lower-dimensional
harrowell2022 . The velocity profile of transiting atoms, presented in Fig. 4,
is calculated over transiting atoms, with no reference to their local
environments, and shows the essence of these transits: a sharp increase in
speed during the transit with rapid decorrelation to the average thermal speed
on either side. Indeed, transiting atoms are not spatially correlated with
each other. Every atom in the fluid will undergo a transit, given enough time,
and the average speed of this atom during the transit will obey the profile
displayed in Fig. 4.

We note that topological quantum matter in $D$ real spatial dimensions can be
described by effective high dimensionality $D+d$, where $d$ are synthetic
dimensions related to additional distinct degrees of freedom due to, for
example, coupling, and arise through reinterpretation of underlying physics
naturereview ; phystoday . This is generally similar to our current finding:
higher effective dimensionality in liquids is related to a distinct (faster)
dynamical sub-ensemble in the system. Differently to higher synthetic
dimensions discussed earlier naturereview , higher dimensionality in liquids
does not need to be integer and is fractional.

##  IV Summary

In summary, we have shown that flow-enabling transit atoms in liquids form a
dynamically distinct sub-ensemble where atoms move on average faster than the
overall system and which has a manifestly non-Maxwellian velocity
distribution. This is in contrast to solids and gases where no such
distinction related to different ensembles can be made. The non-Maxwellian
distribution of transits can be described by an exponent $\alpha$
corresponding to high space dimensionality. This dimenstionality is close to 4
at melting and exceeds 4 at high temperature. $\alpha$ in the liquid state
interestingly has a maximum as a function of temperature and pressure,
returning to its Maxwell value in the solid and gas states.

##  V Data Availability

The authors declare that the data supporting the findings of this study are
available within the paper.

##  VI Acknowledgements

We thank Peter Harrowell for critical reading of the manuscript and for useful
comments.

##  VII Author Contributions

All authors conceptualised the study. C.C. performed simulations and
subsequent data analysis. All authors designed and wrote the manuscript.

##  VIII Competing Interests

The authors declare no competing interests.

## References

  * (1) L. D. Landau and E. M. Lifshitz.  Statistical Physics, part 1. Pergamon Press, 1970.
  * (2) L. P. Pitaevskii.  Statistical physics. liquids.  Physical Encyclopedia, ed. by A. M. Prokhorov, 4:670, 1994.
  * (3) A. I. Akhiezer.  Memories about jacob ilyich frenkel.  Low Temp. Phys., 20:194, 1994.
  * (4) K. Trachenko and V. V. Brazhkin.  Collective modes and thermodynamics of the liquid state.  Reports on Progress in Physics, 79:016502, 2016.
  * (5) J. Frenkel.  Kinetic Theory of Liquids.  Oxford University Press, 1947.
  * (6) D. C. Wallace, E. D. Chisolm, and B. E. Clements.  Observation of single transits in supercooled monatomic liquids.  Phys. Rev. E, 64:011205, 2001.
  * (7) J. C. Dyre.  Colloquium: The glass transition and elastic models of glass-forming liquids.  Rev. Mod. Phys., 78:953, 2006.
  * (8) J. P. Hansen and I. R. McDonald.  Theory of Simple Liquids.  Elsevier, 2013.
  * (9) U. Balucani and M. Zoppi.  Dynamics of the Liquid State.  Oxford University Press, 2003.
  * (10) J. R. D. Copley and J. M. Rowe.  Short-wavelength collective excitations in liquid rubidium observed by coherent neutron scattering.  Phys. Rev. Lett., 32:49, 1974.
  * (11) W. C. Pilgrim, S. Hosokawa, H. Saggau, H. Sinn, and E. Burkel.  Temperature dependence of collective modes in liquid sodium.  J. Non-Cryst. Sol., 250-252:96, 1999.
  * (12) E. Burkel.  Phonon spectroscopy by inelastic x-ray scattering.  Rep. Prog. Phys., 63:171, 2000.
  * (13) W. C. Pilgrim and C. Morkel.  State dependent particle dynamics in liquid alkali metals.  J. Phys.: Condens. Matt., 18:R585, 2006.
  * (14) A. Cunsolo, C. N. Kodituwakku, F. Bencivenga, M. Frontzek, B. M. Leu, and A. H. Said.  Transverse dynamics of water across the melting point: a parallel neutron and x-ray inelastic scattering study.  Phys. Rev. B, 85:174305, 2012.
  * (15) S. Hosokawa et al.  Transverse acoustic excitations in liquid ga.  Phys. Rev. Lett., 102:105502, 2009.
  * (16) S. Hosokawa, M. Inui, Y. Kajihara, S. Tsutsui, and A. Q. R. Baron.  Transverse excitations in liquid sn.  J. Phys.: Condens. Matt., 27:194104, 2015.
  * (17) V. M. Giordano and G. Monaco.  Fingerprints of order and disorder on the high-frequency dynamics of liquids.  PNAS, 107:21985, 2010.
  * (18) V. M. Giordano and G. Monaco.  Inelastic x-ray scattering study of liquid ga: Implications for the short-range order.  Phys. Rev. B, 84:052201, 2011.
  * (19) S. Hosokawa et al.  Transverse excitations in liquid sn.  J. Phys.: Condens. Matt., 25:112101, 2013.
  * (20) Sydney Chapman and T. G. Cowling.  The mathematical theory of non-uniform gases.  Cambridge University Press, 1990.
  * (21) T. Ozawa and H. M. Price.  Topological quantum matter in synthetic dimensions.  Nature Rev. Phys., 1:349, 2019.
  * (22) K. R. A. Hazzard and B. Gadway.  Synthetic dimensions.  Physics Today, 76(4):62, 2023.
  * (23) I. T. Todorov, B. Smith, M. T. Dove, and K. Trachenko.  J. Mater. Chem., 16:1911, 2006.
  * (24) S. Sengupta, F. Vasconcelos, F. Affouard, and S. Sastry.  Dependence of the fragility of a glass former on the softness of interparticle interactions.  The Journal of Chemical Physics, 135:194503, 2011.
  * (25) T. Iwashita, D. M. Nicholson, and T. Egami.  Elementary excitations and crossover phenomenon in liquids.  Physical Review Letters, 110:205504, 2013.
  * (26) C. Cockrell, V. V. Brazhkin, and K. Trachenko.  Transition in the supercritical state of matter: review of experimental evidence.  Physics Reports, 941:1, 2021.
  * (27) A. Samanta, M. E. Tuckerman, T. Q. Yu, and E. Weinan.  Microscopic mechanisms of equilibrium melting of a solid.  Science, 346:729, 2014.
  * (28) M. D. Ediger.  Spatially heteregogeneous dynamics in supercooled liquids.  Annu. Rev. Phys. Chem., 51:99, 2000.
  * (29) K. Vollmayr-Lee, W. Kob, K. Binder, and A. Zippelius.  Dynamical heterogeneities below the glass transition.  J. Chem. Phys., 116:5158, 2002.
  * (30) M. S. G. Razul, G. S. Matharoo, and P. H. Poole.  Spatial correlation of the dynamic propensity of a glass-forming liquid.  J. Phys.: Condens. Matt., 23:235103, 2011.
  * (31) G. Sun and P. Harrowell.  A general structural order parameter for the amorphous solidification of a supercooled liquid.  J. Chem. Phys., 157:024501, 2022.


#  A Recommender System for Equitable Public Art Curation and Installation

Anna Haensch  Tufts University, Data Intensive Studies Center
<EMAIL_ADDRESS>,  Abiy Tasissa  Tufts University, Department of
Mathematics<EMAIL_ADDRESS>and  Dina Deitsch  Tufts University Art
Galleries<EMAIL_ADDRESS>
(Date: 19 July 2022)

###### Abstract.

The placement of art in public spaces can have a significant impact on who
feels a sense of belonging. In cities, public art communicates whose interests
and culture are being favored. In this paper, we propose a graph matching
approach with local constraints to build a curatorial tool for selecting
public art in a way that supports inclusive spaces. We develop a cost matrix
by drawing on Schelling’s model of segregation. Using the cost matrix as an
input, the optimization problem is solved via projected gradient descent to
obtain a soft assignment matrix. We discuss regularization terms to set
curatorial constraints. Our optimization program allocates artwork to public
spaces and walls in a way that de-prioritizes “in-group” preferences, by
satisfying minimum representation and exposure criteria. We draw on existing
literature to develop a fairness metric for our algorithmic output. Using
Tufts University as a testbed, we assess the effectiveness of our approach and
discuss its potential pitfalls from both a curatorial and equity standpoint.

##  1\. Introduction

The placement of art in public spaces can have a significant impact on the
individuals who view it. In the urban landscape, public art plays an important
role in establishing the inclusionary/exclusionary role of urban regeneration,
signaling whose interests and culture are being favored [6]. An academic
campus acts as a kind of city with its own physical configuration and
established hierarchies, and the design of these physical spaces can deeply
inform one’s sense of belonging [1]. Historically, most of the campuses across
the United States were designed to support a majority white, male academic
community. More recently, there has been a movement towards creating more
inclusive campuses, which call for a broad re-imaginging of the spaces of
campus and the people they serve [7]. As part of this, many college campuses
have been turning to the visual arts as a means to address representation,
using campus and public art collections to assist in dismantling the
exclusionary spaces on campus.

Figure 1. Five Mauve Heads, Two Yellow Heads, Mary Tillman Smith, 1989

Much of the historic art on the Tufts University campus highlights the
predominantly white, male history of the university. The university’s
permanent collection, catalogued in the Tufts University Art Galleries Artist,
Subject, $\&$ Donor Database [9], consists of 2,146 works, the great majority
of which are created by white men (see figure 3). This disproportionate
representation, especially when restricted to the historic university
portraiture collection is a reflection of the institution’s past. The
collection is not enyclopedic, but rather it reflects the legacy of the
institution. Works such the portrait of Dr. Morton Henry Prince, a former
chairman of the departments of psychology and neurology, seen in Figure 2,
comprise a significant portion of the portrait collection. This work is
important to the identity of the university both because of the institutional
history it captures, but also because of the prominence of its creator, the
renown American portraitist, John Singer Sargent (1856-1925). However,
preserving and displaying this legacy falls out of step with the increasingly
loud demands from the academic community to create a more equitable and
inclusive space. A work such as Figure 1 by the prominent Mississippi born
Black American artist Mary Tillman Smith (1904-1995), while not a portrait,
depicts human figures from a viewpoint that has potential to connect with
viewers and subtly change the perception of both belonging on campus and the
value of different artistic forms.

Comparing the demographics in the campus art collection to the campus body, we
immediately observe from Figure 3, that the gender representation through
artists shown on campus is not consistent with the gender representation
within the campus student body [10]. This naturally raises a question of the
impact of the art hanging on campus and how optimizing the choices of campus
art in public view might shape narratives about who belongs at the university
[8].

Figure 2. Dr. Morton Henry Prince, John Singer Sargeant, 1895.

In this paper, we explore one particular notion of optimal art hanging,
motivated by the work of economist Thomas Schelling. Schelling’s model of
segregation is an agent based model which demonstrates how individual
preference for neighbors of the same race, even when mild, leads to
segregation [5]. Schelling’s work shows that in order to desegregate
neighborhoods it is necessary to not only be tolerant of people different from
oneself, but actively seek out diverse neighbors. Applying these principles to
representation in public art, we develop our model of hanging art with a
strong preference against “in group” representation.

Formally we approach this as a graph matching problem with local constraints.
First we consider the buildings on campus as nodes on a graph, and consider
the works of art in the collection as nodes on another graph. Classically, the
graph matching problem would endeavor to match the art graph with the building
graph in a way that preserves closeness within the graphs respectively. In our
approach, we also care about the local characteristics of the nodes, and
therefore create a graph matching that preserves with inter- and intra- graph
notions of closeness. This will all be formalized in section 2.1.

While the questions of value in art can be framed in multiple ways, through
this work we hope to expand upon notions of value in art as a reflection of
community identity, and as a mechanism for rethinking community belonging. We
will use the graph matching algorithm to create an art recommendation system
which we call OpArt: The Optimal Art Curation Tool. In contrast to
conventional assignment algorithms that match students to art pieces that most
resemble their attributes and preferences, our optimization problem directly
integrates diversity in the model and attempts to bring algorithmic fairness
to the assignment problem. In section 2 we describe our data sources and our
methods for formalizing the problem. In section 2.1 we define our model
parameters and describe the OpArt optimization algorithm. In section 3 we
present validation metrics and describe the results from our optimized model.
Finally, in section 4 we discuss the implications of the model output in the
context of curatorial practice and public art. All of the tools and and data
for this project are publicly available on Github at
https://github.com/annahaensch/OpArt_Curation_Tool.

##  2\. Methods

We frame our inquiry as a problem of resource allocation, if we consider the
objects in the collection as resources that are being allocated to points on
campus for display. Of the 2,146 works in the Tufts permanent collection,
there are 392 works on display across buildings on Tufts’ Medford and Boston
campuses. There are 23 buildings on campus that regularly display works on
loan from the permanent collection. This includes the main campus library
which has 58 pieces on loan, the presidential residence which has 24 pieces on
loan, as well as numerous other academic, residential, and administrative
buildings.

Figure 3. Left: Each dot represents one work of art in the Artist, Subject, &
Donor Database where position and color indicated the inferred race and gender
of the creator. All identity attributes are inferred rather than self-
reported, consequently the data does not capture the full diversity of artists
in the collection. Right: Each dot represents one student registered in the
Tufts Fall 2021 Enrollment Calculator [10]. Values are self-reported, although
limitations in university reporting and classification often suppress certain
non-binary gender identities. To acknowledge the limitations in data
collection, a label has been created for transgender and non-binary students
and artists, although none currently appear among the inferred/reported gender
identities in either database.

For each of these 23 buildings with art on loan, we can approximate who passes
through the building each day by determining which subsets of the student body
are served by this building and how. This is achieved by running a script
which queries the Tufts University mobile mapping site
https://m.tufts.edu/tufts_mobile/map_all/ to gather building names, GPS
coordinates, and Office & Departments. From this we create a one-hot table
which classifies buildings by school and student population served. Buildings
can contain academic space (e.g. faculty offices, classrooms, labs) for one of
the 11 specific schools at the university, public use space (e.g. libraries,
auditoriums, dining), administrative space (e.g. registrar, parking office) or
residential space.

There are currently 13,293 students enrolled across the 11 schools of the
university. Using the Tufts Fall 2021 Enrollment Calculator [10] we construct
a data frame where each row corresponds to one enrolled student with student
school affiliation, race, and gender as the features of interest. With this,
we fill the space according to the following rules:

  1. 1.

Each student is assigned to all of the buildings that correspond to their
primary school of enrollment

  2. 2.

1% of students at random are assigned to buildings that house administrative
space.

  3. 3.

2% of students at random are assigned to buildings that house public space

  4. 4.

Students are split among residence halls at random according to the known
number of available beds.

We note that students can be assigned to multiple buildings at once (since
people often pass between buildings during a given day). Once the students are
assigned to buildings, the buildings are filled. To deal with the randomness
this assignment, we carry out this filling procedure multiple times to obtain
averages.

Instructions for recreating the results here are included in the README of the
OpArt Github repository. In the following section we will generalize the
methodology described above to the more flexible framework of objects (e.g.
pieces of art) to locations (e.g. buildings).

###  2.1. Mathematical model

In our model, we consider the assignment of $M$ objects to $N$ locations,
where the set of all objects is given by

$\mathcal{A}=\\{\mathbf{a}_{1},...,\mathbf{a}_{M}\\}$

and the set of all locations by

$\mathcal{B}=\\{\mathbf{b}_{1},...,\mathbf{b}_{N}\\}.$

We will suppose that there are more objects than locations, that is, $M>N$. We
will also suppose that each location, $\mathbf{b}_{n}$, has a unique capacity
for objects given by $h(\mathbf{b}_{n})$. In this way,

$\sum_{n=1}^{N}h(\mathbf{b}_{n})=\text{total number of objects needed},$

and we defined $\mathbf{h}:=(h(\mathbf{b}_{1}),...,h(\mathbf{b}_{N}))$.
Similarly, we will suppose that each object $\mathbf{a}_{m}$ has a unique
capacity for assignment given by $k(\mathbf{a}_{m})$, whereby

$\sum_{m=1}^{M}k(\mathbf{a}_{m})=\text{total number of objects available},$

and we define $\mathbf{k}:=(k(\mathbf{a}_{1}),...,k(\mathbf{a}_{M}))$. We will
also consider users (e.g. the students) as,

$\mathcal{S}=\\{\mathbf{s}_{1},...,\mathbf{s}_{T}\\}$

where $T$ is the total size of the user base. Moreover, for each location we
will define an associate subset of users,
$\mathcal{S}_{n}\subseteq\mathcal{S}$ which denote the set of all users
passing through location $\mathbf{b}_{n}$ in a given day. In the specific case
of students on campus, the precise methodology used to compose the sets
$\mathcal{S}_{n}$ will be described in Section 2 above, but we will note here
that the $\mathcal{S}_{n}$ are not necessarily disjoint.

To further characterize the locations, we will consider the users as vectors
in the $D$-dimensional space

$\mathcal{X}=\mathcal{X}_{1}\times...\times\mathcal{X}_{D}$

where each $\mathcal{X}_{i}$ is a categorical space with
$\mid\mathcal{X}_{i}\mid=T_{i}$. In the case of students in buildings, these
could be categorical assignments for e.g. race and gender or in the case of
artwork these could be categorical assignments for e.g. artist race and
gender. In this way, user $s_{i}\in\mathcal{S}$ and object
$a_{i}\in\mathcal{A}$, can be written as the vectors,

$\mathbf{s}_{i}=(s_{i1},...,s_{iD})\,\,\,\text{ and
}\,\,\,\mathbf{a}_{i}=(a_{i1},...,a_{iD}).$

For user $\mathbf{s}_{i}\in\mathcal{S}$ we will also define the quantized
vector

$Q_{n}(\mathbf{s}_{i}):=(q_{n}(s_{i1}),...,q_{n}(s_{iD}))\in[0,1]^{D},$

where $q_{n}(s_{id})$ is the proportional share of users in $\mathcal{S}_{n}$
with attribute $s_{id}$ in dimension $d$, in other words,

$q_{n}(s_{id})=\frac{\mid\\{\mathbf{s}_{j}\in\mathcal{S}_{n}:s_{jd}=s_{id}\\}\mid}{\mid\mathcal{S}_{n}\mid}$

Similarly, for an object $\mathbf{a}_{i}\in\mathcal{A}$ we define the
quantized vector

$Q_{\mathcal{A}}(\mathbf{a}_{i}):=(q_{\mathcal{A}}(a_{i1}),...,q_{\mathcal{A}}(a_{iD}))\in[0,1]^{D},$

where $q_{\mathcal{A}}(a_{id})$ is the proportional share of objects in
$\mathcal{A}$ with attribute $a_{id}$ in dimension $D$, in other words,

$q_{\mathcal{A}}(a_{id})=\frac{\mid\\{\mathbf{a}_{j}\in\mathcal{A}:a_{jd}=a_{id}\\}\mid}{\mid\mathcal{A}\mid}.$

Since it will be of use later on, we define component-wise product of this
quantized vector as

(1) $\displaystyle\rho$ $\displaystyle:=$
$\displaystyle\prod_{d=1}^{D}q_{\mathcal{A}}(a_{id})$

Of interest is to find the “optimal” assignment of the objects in
$\mathcal{A}$ to the locations $\mathcal{B}$ that maximizes the reach of
artwork with respect to the relevant student body. This optimality will be
made precise in the discussion to follow.

We define a metric $\|\cdot\|_{\mathcal{S}}$ which will be based on a weighted
norm that takes into account the diversity attributes in $\mathcal{B}$ and
$\mathcal{S}$. For a vectors $\bf x,\bf y\in\mathcal{X}$ we define our metric
as

$\displaystyle\|{\bf x}-{\bf
y}\|_{n}:=\|\delta_{x_{1},y_{1}}\cdot(q_{n}(x_{1})-q_{n}(y_{1})),...,\delta_{x_{D},y_{D}}\cdot(q_{n}(x_{D})-q_{n}(y_{D}))\|$

where $\|\cdot\|$ denotes the usual Euclidean norm, and $\delta$ denotes the
Kronecker delta function. In this way, vectors ${\bf x}$ and ${\bf y}$ which
share all attributes, and therefore correspond to the same point in
$\mathcal{X}$, will have distance 0 from one another. In cases where some
attributes are not shared, the contribution to the distance will be small for
any attribute, $\mathcal{X}_{d}$, where $x_{d}$ and $y_{d}$ come from classes
that hold a similar proportional share of that attribute. Here on, closeness
is based on this weighted metric.

###  2.2. OpArt optimization program

We now design an optimization program based on the following objective: we
want to assign the art pieces to the buildings in such a way that we maximize
dissimilarity. However, care is needed to further put restrictions on this
objective. Otherwise, an optimal solution of the optimization program would be
to assign the most dissimilar art piece to a given building. Given this, our
goal is to increase exposure to new items for the users of the building and
not necessarily preclude from placing an art piece consistent with the
demographics. In particular, our dissimilarity measure aims to capture this
using weights which inform the extent to which we want to balance exposure to
diverse items with exposure to a familiar art piece. More specifically, the
similarity of a building to a given art piece depends on two factors. The
first factor simply considers the similarity of the attributes of the
particular art piece and attributes of the occupants of the building. The
second factor weights this similarity by considering the extent to which a
given occupant is a representative/outlier of that building. To detect outlier
occupants in a particular building, we employ the deviation of their
attributes from the mode of the building where the latter is computed by
considering all occupants, and will be denoted $\text{mode}(\mathbf{b}_{n})$.
While we employ mode based outlier detection algorithms [4], our goal here is
not to remove the outliers but rather identify them and define appropriate
weighted metric.

We propose the following cost function,

(2) $\displaystyle C(\mathbf{b}_{n},\mathbf{a}_{m})$ $\displaystyle=$
$\displaystyle\frac{\exp\left\\{\sum_{\mathbf{s}\in\mathcal{S}_{n}}\frac{-\alpha\cdot\|\mathbf{s}-\text{mode}(\mathbf{b}_{n})\|_{n}}{\beta\cdot\rho\cdot\|\mathbf{s}-\mathbf{a}_{m}\|}\right\\}}{\sum_{i=1}^{M}\exp\left\\{\sum_{\mathbf{s}\in\mathcal{S}_{n}}\frac{-\alpha\cdot\|\mathbf{s}-\text{mode}(\mathbf{b}_{n})\|_{n}}{\beta\cdot\rho\cdot\|\mathbf{s}-\mathbf{a}_{i}\|}\right\\}}$

where $\alpha$ and $\beta$ are positive parameters to be tuned and $\rho$ is
as in equation (1). This cost matrix is derived from the standard multinomal
logistic regression model, where $C(\mathbf{b}_{n},\mathbf{a}_{m})$ is taken
as the probability of object $\mathbf{a}_{m}$ given location $\mathbf{b}_{n}$.
This choice of cost function is motivated by the observation that the distance
of a particular user to an object is exponentially weighted by considering the
extent to which the user is representative of the location. The cost matrix is
then defined as $\mathbf{C}$ where the $nm^{\text{th}}$ entry is
$C(\mathbf{b}_{n},\mathbf{a}_{m})$. The proposed optimization program is

(3) $\underset{\mathbf{P}\in\mathcal{R}^{N\times
M},\mathbf{P}\succeq\mathbf{0},\mathbf{P}\mathbf{1}=\mathbf{h}}{\min}\quad\textrm{trace}(\mathbf{C}^{\intercal}\mathbf{P}),$

where $\mathbf{P}$ is an assignment matrix and $\mathbf{P}\mathbf{1}$ denotes
the $m$-fold sum taken over the columns of $\mathbf{P}$. We note that
$\sum_{m}\mathbf{P}_{n,m}=h(\mathbf{b}_{n})$ for all $1\leq n\leq N$. To avoid
the case where a given object is overly allocated, we further limit the number
of times an object can be utilized. We do this via a penalty term that
discourages the $m^{\text{th}}$ object from being assigned to more than
$\mathbf{k}_{m}$ locations where $\mathbf{k}$ is as defined above. With this,
the full optimization program is given by

(4)
$\underset{\mathbb{S}}{\min}\quad\textrm{trace}(\mathbf{C}^{\intercal}\mathbf{P})+\frac{\lambda}{2}\|\mathbf{P}^{\intercal}\mathbf{1}-\mathbf{k}\|_{F}^{2},$

where $\mathbb{S}$ is the set defined as
$\mathbb{S}=\\{\mathbf{P}\in\mathcal{R}^{N\times
M}:\mathbf{P}\succeq\mathbf{0},\mathbf{P}\mathbf{1}=\mathbf{h}\\}$ and
$\lambda$ is a regularization parameter to be tuned. The optimization program
in (4) is convex. We use a projected gradient descent scheme to solve the
algorithm in (4). The algorithm alternates between taking a gradient descent
step and projecting onto the set $\mathbb{S}$. We use the algorithm in [2] to
exactly project onto $\mathbb{S}$.

Algorithm 1 A projected gradient descent algorithm to solve (4)

1:Initialization: Set $\mathbf{P}^{(0)}$, $\lambda$, $\mathbf{k}$,
$\mathbf{h}$, $\epsilon$ (step size), maxiterations (maximum number of
iterations)

2:for k = 1:maxiterations do

3: Gradient descent:
$\bar{\mathbf{P}}^{(k)}=\mathbf{P}^{(k-1)}-\epsilon(\mathbf{C})-\epsilon\lambda\mathbf{1}(\mathbf{1}^{\intercal}\mathbf{P}^{(k-1)}-\mathbf{k}^{\intercal}).$

4: Project $\bar{\mathbf{P}}^{(k)}$ onto the set $\mathbb{S}$ and obtain
$\mathbf{P}^{(k)}$ as follows:

5: for i= 1:N do

6: Sort the entries of row $i$ of $\bar{\mathbf{P}}$ in decreasing order:
$u_{1}\geq u_{2}\geq...\geq u_{M}$.

7: Define $L$: $L=\underset{1\leq k\leq
M}{\max}\,\\{k\mid\frac{\sum_{r=1}^{k}u_{k}-h_{n}}{k}<u_{k}\\}$

8: Define $\tau$: $\tau=\frac{\sum_{\ell=1}^{L}u_{\ell}-\mathbf{h}_{n}}{L}$

9: Update row $i$ of $\bar{\mathbf{P}}^{(k)}$:
$(\mathbf{P}^{(k)})(i,j)=\max\left((\bar{\mathbf{P}}^{(k})(i,j)-\tau,0\right)$
for all $j=1:M$

10: end for

11:end for

Below, we summarize important aspects of the algorithms: initialization,
convergence, step size and stopping criterion.

Initialization: Since (4) is a convex optimization program, all the minimizers
$\mathbf{P}^{*}$, i.e. the optimal assignment matrices, have the same
objective value. However, these global minimizers themselves could be
different from each other. Given that, it is of interest to assess the role of
initialization on what kind of global minimizer is attained. Towards
quantifying this, we consider three initializations for all experiments
discussed in the paper. The first is initializing the permutation based on the
ground truth assignment. In our setting, the ground truth refers to the
existing assignment of artworks to buildings. Another initialization we
consider is uniform assignment where each art piece is equally likely to be
assigned to each building. Finally, we initialize by sampling uniformly
randomly from the set $\mathbb{S}$ using the algorithm in [12].

Convergence and step size: To choose the step size, we first compute the
Lipschitz constant of the gradient of the objective function in (4). This
quantity is useful in establishing the convergence of gradient descent. Let
$F(\mathbf{P})$ denote the objective function in (4). The gradient of $F$ is
Lipschitz continuous with parameter $L>0$ if $\|\nabla
F(\mathbf{P}_{1})-\nabla
F(\mathbf{P}_{2})\|_{F}\leq\|\mathbf{P}_{1}-\mathbf{P}_{2}\|_{F}$ for all
assignment matrices $\mathbf{P}_{1},\mathbf{P}_{2}$. It can be verified that
the Lipschitz constant of the gradient of $F$ is $L=\lambda N$. We use the
following standard result on the convergence of projected gradient descent
[13].

######  Theorem 1 ([13] (Theorem 3.7)).

The projected gradient descent algorithm with $\epsilon=\frac{1}{L}$ satisfies

$f(\mathbf{P}^{(k)})-f(\mathbf{P}^{*})\leq\frac{3\epsilon\|\mathbf{P}^{(0)}-\mathbf{P}^{*}\|_{F}^{2}+f(\mathbf{P}^{0})-f(\mathbf{P}^{*})}{k},$

where $\mathbf{P}^{(0)}$ denotes the initial assignment matrix,
$\mathbf{P}^{*}$ denotes the optimal solution and $\mathbf{P}^{(k)}$ denotes
the estimate of the assignment matrix at the k-th iteration.

Given the above convergence result, for all our experiments, we set
$\epsilon=\frac{1}{2N}$.

Stopping criterion: The stopping criterion is the maximum number of iterations
set to $1000$ which for all experiments is consistent with convergence in
objective value.

###  2.3. Leveraging assignment priors

From a practical point of view, if the optimization algorithm obtains an
assignment matrix that differs significantly from the current assignment,
significant costs might be incurred to make the necessary changes. With that,
it might be useful in certain cases to make a gradual change. Rather than
seeking globally optimal solutions that might be hard to realize in a timely
manner, the practitioner might be interested in improving the current
assignment gradually. To reflect this in our model, we make use of the
existing assignment matrix, denoted by $\mathbf{P}^{\text{current}}$ and
modify (4) as follows

(5)
$\underset{\mathbb{S}}{\min}\quad\textrm{trace}(\mathbf{C}^{\intercal}\mathbf{P})+\frac{\lambda}{2}\|\mathbf{P}^{\intercal}\mathbf{1}-\mathbf{k}\|_{F}^{2}+\frac{\tau}{2}\|\mathbf{P}-\mathbf{P}^{\text{current}}\|_{F}^{2}.$

Above the last term is a regularization term that encourages the optimized
assignment matrix to not be far away from the current assignment. From an
optimization perspective, (7) is strongly convex for which there is a unique
globally optimal solution. By varying the parameter $\tau$, the practitioner
will obtain assignment matrices that are close to the current assignment while
also optimizing the assignment problem of matching art pieces to buildings.
Complete details of the optimization program are given in Appendix References.

##  3\. Results

All the algorithms in this paper are implemented in Python. For the
experiments, it is necessary to suitably set the hyper-parameters $\alpha$,
$\beta$, $\lambda$ and $\tau$ that determine the optimal assignment matrix to
be obtained from solving (7). For all our numerical experiments, $\alpha=-1$.
For the experiments, we set $\beta$ to values in the range $10^{-1}$ to
$10^{15}$. For a fixed value of $\beta$, before setting $\lambda$ and $\tau$,
it is crucial to find the scaling of the terms in the objective (7). In
particular, let
$f_{1}(\mathbf{P})=\textrm{trace}(\mathbf{C}^{\intercal}\mathbf{P})$,
$f_{2}(\mathbf{P})=\|\mathbf{P}^{\intercal}\mathbf{1}-\mathbf{k}\|_{F}^{2}$
and $f_{3}(\mathbf{P})=\|\mathbf{P}-\mathbf{P}^{\text{current}}\|_{F}^{2}$
denote the first three terms of the objective in (7) respectively. For a fixed
value of $\alpha$ and $\beta$, $f_{1}$ is determined. To determine appropriate
scaling for $\lambda$ and $\tau$, we sample $r$ assignment matrices uniformly
from the set $\mathbb{S}$ and compute $f_{2}$ and $f_{3}$ for each sampled
matrix. We now define the scaling factors $\lambda_{s}$ and $\tau_{s}$ as
follows: $\lambda_{s}=\frac{1}{r}\sum_{i=1}^{r}f_{2}(i)/f_{1}$ and
$\tau_{s}=\frac{1}{r}\sum_{i=1}^{r}f_{3}(i)/f_{1}$. In our experiments, we
have used $r=50$. The scaling factors ensure that the three terms are
comparable. With that, to set $\lambda$ and $\tau$, we proceed as
$\lambda=\lambda_{s}\cdot\bar{\lambda}$ and $\tau=\tau_{s}\cdot\bar{\tau}$
where $\bar{\lambda}$ and $\bar{\tau}$ are multiplicative factors that are
chosen from the range $[1,10000]$.

The goal of this work is to allocate objects to locations with the objective
or prioritizing a self-representative allocation for individuals who are both
historically underrepresented in terms of the object and the location. To
evaluate the success of this objective we will use an unfairness metric
similar to what is introduced in [11]. We define $\mathcal{G}$ to be the
disadvantaged group, and $\neg\mathcal{G}$ to be the advantaged group, so
$\mathcal{G}\cup\neg\mathcal{G}=\mathcal{S}$.

We define the fairness of our assignment as

$\displaystyle U=E_{\mathcal{G}}[y]-E_{\neg\mathcal{G}}[y]$

where $E_{\mathcal{G}}[y]$ as average number of self-representative objects
seen by a person in group $\mathcal{G}$,

(6) $\displaystyle E_{\mathcal{G}}[y]$ $\displaystyle=$
$\displaystyle\frac{1}{\mid{\mathcal{G}}\mid}\sum_{\mathbf{s}\in{\mathcal{G}}}r(\mathbf{s})$

where $r(\mathbf{s})$ is the expected number of representative objects seen by
person $\mathbf{s}$ which can be computing as an expected value conditioned on
the buildings. These results are shown for variable optimization methods in
Table 1. The primary difference between our fairness metric and that proposed
in [11] is that we don’t take the absolute value since our present work is
specifically focused on elevating underrepresented groups.

Table 1. Mean $\pm$ standard deviation for $E_{\mathcal{G}}[y]$ for baseline (i.e. current) and optimized assignment where we consider $\mathcal{G}$ as the disadvantaged group by gender (i.e. “Non-Man”) and race (i.e. “Non-White”). The first row shows the baseline, and each subsequent row shows an optimization with a different initialization procedure setting $P^{0}$ equal to a uniform matrix, the current assignment matrix, or a random projection. Statistics in this table are gathered over 50 re-samplings of the student body. | Man | Non-Man | White | Non-White
---|---|---|---|---
Baseline | $12.280\pm 0.134$ | $2.456\pm 0.021$ | $13.309\pm 0.129$ | $1.591\pm 0.016$
Uniform Initialization | $10.203\pm 0.111$ | $4.018\pm 0.032$ | $6.594\pm 0.065$ | $7.614\pm 0.068$
Current Initialization | $10.869\pm 0.119$ | $3.375\pm 0.044$ | $7.780\pm 0.107$ | $6.552\pm 0.059$
Random Projection | $10.357\pm 0.273$ | $3.871\pm 0.259$ | $6.860\pm 0.523$ | $7.362\pm 0.474$

From Table 1 we can clearly see that our proposed methods of optimization not
only increase the expected average self-representative objects seen by
individuals in the disadvantaged group, but also increase fairness, $U$, over
the baseline. Moreover, the overall fairness of the optimization is not
particularly sensitive to our choice of initialization and is relatively
stable over multiple rounds of sampling, where each round of sampling can be
considered as a different “day” on campus. The first row of the table gives
the baseline stats for the current assignment, and subsequent rows apply the
different initialization procedures described in section 2.2. Before carrying
out simulations, we scale $\lambda$ (resp. $\tau$) by the order of magnitude
of quotient of the first and second (resp. first and third) terms of 7 to
ensure that the scales are comparable between terms. We replace ${\lambda}$
and ${\tau}$ with the scaled coefficients, $\overline{\lambda}$ and
$\overline{\tau}$.

###  3.1. Group structure of optimal assignment matrices

As remarked in the previous section, the optimization program in (4) is convex
with all the minimizers sharing the same objective value. For a varying set of
the hyper-parameters, it is of interest to assess how “different” the optimal
permutation matrices are. To realize this goal, we use the spectral embedding
algorithm in [14]. Spectral embedding takes as an input a similarity matrix
and outputs an embedding, in a certain dimension, in terms of point
coordinates. In our setting, the similarity matrix would inform us how two
optimal assignment matrices are similar to each other. First, under different
choices of the hyper-parameters, we generate $r$ optimal permutation matrices.
Let $\mathbf{W}\in\mathcal{R}^{r\times r}$ denote the similarity matrix with
$\mathbf{W}_{i,j}$ denoting the similarity between the $i$-th optimal
permutation matrix and the $j$-th optimal permutation matrix. For our
experiments, $\mathbf{W}$ is defined as follows:
$\mathbf{W}_{i,j}=\frac{1}{\|\mathbf{P}_{i}-\mathbf{P}_{j}\|_{1}}$ if $i\neq
j$ and $0$ if $i=j$ where $\|\cdot\|_{1}$ the entry-wise $\ell_{1}$ norm. A
larger entry would indicate that the two assignment matrices are more similar.
Given this, we run the spectral embedding algorithm and set the embedding
dimension to $2$ for visualization purposes. In Figure 4, we show the spectral
embeddings results. The first case concerns when $\bar{\lambda}$ is fixed and
$\bar{\tau}$ and $\beta$ are varied. As $\bar{\lambda}$ is relatively larger,
the art capacity constrains are met. We can clearly see that the different
permutation matrices cluster based on the size of $\bar{\tau}$. In the right
most plot, $\bar{\lambda}$ and $\bar{\tau}$ are fixed while $\bar{\beta}$ is
varied. For this setting, we observe that there are different clusters
corresponding to the different range of $\beta$.

Figure 4. Spectral embedding of optimal assignment matrices by varying the
hyper-parameters: Left, Middle: $\beta$ and $\overline{\tau}$ for fixed
$\overline{\lambda}$, Right: $\beta$ while $\overline{\tau}$ and
$\overline{\lambda}$ is fixed.

##  4\. Discussion

In Figure 5 we can see how the relative importance of model parameters are
born out in the optimization algorithm. In particular, when $\lambda$ is
small, the algorithm ignores artwork capacity constraints, which means that
assignments might perform quite well in terms of fairness (see for example row
3, column 1 of Figure 5 ), but are in-practice impossible, since they assign
artworks beyond what the collection contains. For example, such an assignment
might place a high priority on artwork by Native Hawaiian Women across several
buildings on campus, despite the fact that the Tufts Permanent collection
holds no such works. Perhaps this is a tool that could be used to inform new
artwork acquisition.

Figure 5. Each plot in the grid above shows the $U$ score by demographic
category. The blue line indicates the score with $\mathcal{G}$ as the ”non-
male” group and the orange line indicates the score with $\mathcal{G}$ as the
”non-white” group. The dashed gray line at 0 is there to highlight when the
assignment exhibits preferential treatment to the underrepresented group.

In Figure 5, we see that as $\tau$ increases the artwork hanging has a strong
preference towards the status-quo. In this case, disadvantaged groups continue
to see very little art that reflects their own identity, until the $\lambda$
terms becomes large enough to overpower this preference for the status-quo
(see for example row 1, column 4). From a curatorial standpoint, however, this
slow movement towards more equitable art might be preferential, since
acquiring and installing new exhibits takes both time and money.

On the one hand this optimization tool is highly effective in highlighting the
glaring blind spots of artists collected in Tufts’ permanent art collection
and can indeed function as a guide for future acquisitions. The history of
collecting at Tufts is a short and passive one in that it has historically
relied on gifts from alumni and would therefore tend to reflect those donors,
who, following the historical demographics of the school, were and are
majority white male.

On the other hand, curatorial work is multifaceted and the selection of
artwork for campus spaces has to take into account a number of factors—from
the basic logistics of security and direct sunlight (which will fade an
artwork over time) to the content and subject matter of an artwork and its
relevancy to the selected location. Similar to the methods used to leverage
assignment priors in Section 3.2, we could add a penalty term to leverage
curatorial preferences by building and artwork. For instance, it might be poor
judgment to hang a small, intimate photograph of a nude figure in a large
dining hall where the space was large and loud and 1000s of students passed
each day. The work would be subject to decay and fading from the light,
possible damage from the proximity to the food, and the subject might not land
well with a group of rowdy students after a long day, eating dinner.
Installing artwork in public spaces – be they a large lawn or library wall –
requires a sensitivity to context that needs to run from the pragmatic to the
conceptual – while yes, taking into account the maker of the artwork and its
capacity to reflect our campus diversity as effectively as humanly possible.

## Acknowledgments

The authors wish to thank the Data Intensive Studies Center at Tufts
University for generous seed funding. The authors also wish to thank The Tufts
University Art Galleries for inviting the authors to explore their collections
and for granting access to their metadata.

## References

  * [1] L. Carvalho, C.G. Freeman, A. Kearney, M. Mentis, R. Martinez-Maldonado.  Spaces of inclusion and belonging: The learning imaginaries of doctoral students in a multi-campus and distance university. Australasian Journal of Educational Technology 34.6 (2018).
  * [2] L. Condat.  Fast projection onto the simplex and the $l_{\mathbf{1}}$ ball. Mathematical Programming 158.1 (2016): 575-585.
  * [3] D. Greenwald.  What Can Data Teach Us about Museum Collections? American Alliance of Museums, 28 Oct. 2021.  https://www.aam-us.org/2020/04/27/what-can-data-teach-us-about-museum-collections/?gclid=Cj0KCQiA3rKQBhCNARIsACUEW_bL5VwnsDav3L8ZYyObnkVVWTSo9Vsl_IRhLp3z4yHnJ9GYgSJMVvgaAkZVEALw_wcB,  last accessed: 02/16/2022.
  * [4] C. Leys, C. Ley, O. Klein, P. Bernard, L. Licata.  Detecting outliers: Do not use standard deviation around the mean, use absolute deviation around the median. Journal of experimental social psychology 49.4 (2013): 764-766.
  * [5] T.C. Schelling.  Micromotives and macrobehavior.  WW Norton & Company, 2006.
  * [6] J. Sharp, V. Pollock, R. Paddison.  Just art for a just city: Public art and social inclusion in urban regeneration. Urban Studies 42.5-6 (2005): 1001-1023.
  * [7] R. Slee.  The irregular school: Exclusion, schooling and inclusive education.  Routledge, 2011.
  * [8] Tufts University.  Diversity and Inclusion. 10 Feb. 2021, https://www.tufts.edu/strategic-themes/diversity-and-inclusion,  last accessed: 07/26/2022.
  * [9] Tufts University Art Gallery.  Tufts Collection: Artist, Subject, and Donor Data. 12 Feb. 2022, https://github.com/annahaensch/Tufts_Art_Gallery_Optimization/blob/main/data/TUAG_Artist_Subject_Donor_Data.xlsx,  last accessed: 07/26/2022.
  * [10] Tufts University Office of Institutional Research.  Fall 2021 Enrollment Calculator. 31 Oct. 2021, https://provost.tufts.edu/institutionalresearch/enrollment/,  last accessed: 07/26/2022.
  * [11] S. Yao, B. Huang.  Beyond parity: Fairness objectives for collaborative filtering.  Advances in neural information processing systems, 30 (2017).
  * [12] Noah Smith, Roy Tromble.  Sampling uniformly from the unit simplex.  Johns Hopkins University, Tech. Rep, 29 (2004).
  * [13] Sébastien Bubeck.  Convex optimization: Algorithms and complexity.  Foundations and Trends in Machine Learning, 8.3-4 (2015): 231–235.
  * [14] Andrew Ng, Michael Jordan, Yair Weiss.  On spectral clustering: Analysis and an algorithm.  Advances in neural information processing systems, 14 (2001).

##  Appendix A Details of optimization program

We revisit the main optimization program.

(7)
$\underset{\mathbb{S}}{\min}\quad\textrm{trace}(\mathbf{C}^{\intercal}\mathbf{P})+\frac{\lambda}{2}\|\mathbf{P}^{\intercal}\mathbf{1}-\mathbf{k}\|_{2}^{2}+\frac{\tau}{2}\|\mathbf{P}-\mathbf{P}^{\text{current}}\|_{F}^{2}.$

We note that $\mathbb{S}=\\{\mathbf{P}\in\mathcal{R}^{N\times
M}:\mathbf{P}\succeq\mathbf{0},\mathbf{P}\mathbf{1}=\mathbf{h}\\}$ is a convex
set. Let $f_{1}(\mathbf{P})=\textrm{trace}(\mathbf{C}^{\intercal}\mathbf{P})$,
$f_{2}(\mathbf{P})=\|\mathbf{P}^{\intercal}\mathbf{1}-\mathbf{k}\|_{F}^{2}$
and $f_{3}(\mathbf{P})=\|\mathbf{P}-\mathbf{P}^{\text{current}}\|_{F}^{2}$
denote the first three terms of the objective in (7) respectively. The first
term is linear in $\mathbf{P}$ and hence convex. The second and third terms
are composition of a convex function with affine functions yielding convex
functions. It follows that the objective in (7) is convex. All in all, we have
a convex program. As discussed in the main part of the paper, we employ
projected gradient descent to optimize (7). From our numerical experiments,
convergence in objective is attained before the maximum iterations set to
$1000$. A typical instance of the convergence plot is shown in Figure 6.

Figure 6. Convergence of the projected gradient descent algorithm with
$\beta=100$, $\overline{\tau}=1$ and $\overline{\lambda}=1$

To derive the step size for projected gradient descent, we rely on the
Lipschitz constant of the gradient. Below, we detail the computation.
Hereafter, the objective in (7) is denoted by $f(\mathbf{P})$.
Computation of $\nabla f_{1}(\mathbf{P})$: We compute $\nabla
f_{1}(\mathbf{P})$ as follows.

$\displaystyle\frac{\partial}{\partial
P_{i,j}}(\textrm{trace}(\mathbf{C}^{\intercal}\mathbf{P}))$
$\displaystyle=\frac{\partial}{\partial
P_{i,j}}\sum_{mn}C_{mn}P_{m,n}=\sum_{mn}C_{m,n}P_{m,n}\delta_{m,i}\delta_{n,j}=C_{i,j}$

Therefore, $\nabla f_{1}(\mathbf{P})=\mathbf{C}$.

Computation of $\nabla f_{2}(\mathbf{P})$: We compute $\nabla
f_{2}(\mathbf{P})$ as follows.

$\displaystyle\frac{\partial}{\partial
P_{i,j}}\left(\|\mathbf{P}^{\intercal}\mathbf{1}-\mathbf{k}\|_{2}^{2}\right)$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
P_{i,j}}\left(\sum_{m=1}^{M}((\mathbf{P}^{\intercal}\mathbf{1})_{m}-k_{m})^{2}\right)$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
P_{i,j}}\left(\sum_{m=1}^{M}\left[\sum_{r=1}^{N}P_{r,m}-k_{m}\right]^{2}\right)$
$\displaystyle=$ $\displaystyle 2\left(\sum_{r=1}^{N}P_{r,j}-k_{j}\right)$
$\displaystyle=$ $\displaystyle 2(\mathbf{P}^{T}\mathbf{1})_{j}-2k_{j}$

We note that the last term is independent of $i$. It follows that $\nabla
f_{2}(\mathbf{P})=2\mathbf{1}\mathbf{1}^{T}\mathbf{P}-\mathbf{1}\mathbf{k}^{T}=2\mathbf{1}(\mathbf{1}^{T}\mathbf{P}-\mathbf{l}^{T})$.

Computation of $\nabla f_{3}(\mathbf{P})$: We compute $\nabla
f_{3}(\mathbf{P})$ as follows.

$\displaystyle\frac{\partial}{\partial
P_{i,j}}\left(|\mathbf{P}-\mathbf{P}^{\text{current}}\|_{F}^{2}\right)$
$\displaystyle=$ $\displaystyle\frac{\partial}{\partial
P_{i,j}}\sum_{m,n}\left(P_{m,n}-\mathbf{P}^{\text{current}}_{m,n}\right)^{2}$
$\displaystyle=$ $\displaystyle
2\sum_{m,n}\left(P_{m,n}-\mathbf{P}^{\text{current}}_{m,n}\right)\cdot\delta_{m,i}\delta_{n,j}$
$\displaystyle=$ $\displaystyle
2\left(P_{i,j}-P^{\text{current}}_{i,j}\right)$

Therefore, $\nabla
f_{3}(\mathbf{P})=2(\mathbf{P}-\mathbf{P}^{\text{current}})$.
We now proceed to compute the Lipschitz constant of the gradient as follows.
For any $\mathbf{P}_{1},\mathbf{P}_{2}$,

$\displaystyle\|\nabla f(P_{1})-\nabla
f(\mathbf{P}_{2})\|_{F}=\|(\lambda\mathbf{1}\mathbf{1}^{T}-\tau\mathbf{I})(\mathbf{P}_{1}-\mathbf{P}_{2})\|_{F}\leq\|\lambda\mathbf{1}\mathbf{1}^{T}-\tau\mathbf{I}\|_{2}\,\|\mathbf{P}_{1}-\mathbf{P}_{2}\|_{F}$

The eigenvalues of $\lambda\mathbf{1}{1}^{T}$ are $0$ and $\lambda N$. It
follows that $\|\nabla f(P_{1})-\nabla
f(\mathbf{P}_{2})\|_{F}\leq(N-\tau)\,\|\mathbf{P}_{1}-\mathbf{P}_{2}\|_{F}$.
Hence, the Lipschitz constant of the gradient is $\lambda N-\tau$.