# Emotion Conditioned Creative Dialog Generation

Khalid Alnajjar
University of Helsinki
Finland
<EMAIL_ADDRESS>
Mika Hämäläinen
University of Helsinki
Finland

<EMAIL_ADDRESS>

###### Abstract

We present a DialGPT based model for generating creative dialog responses that
are conditioned based on one of the following emotions: anger, disgust, fear,
happiness, pain, sadness and surprise. Our model is capable of producing a
contextually apt response given an input sentence and a desired emotion label.
Our model is capable of expressing the desired emotion with an accuracy of
0.6. The best performing emotions are neutral, fear and disgust. When
measuring the strength of the expressed emotion, we find that anger, fear and
disgust are expressed in the most strong fashion by the model.

##  1 Introduction

Dialog systems and different kinds of natural language interfaces are all
around us. When we seek information or contact customer support, we are
increasingly more often first greeted by a bot rather than a person. Bots are
typically stiff and not lifelike making communication with them an awkward
experience. This is because their goal oriented nature typically sets some
constraints in terms of how creative a system can be.

We can perceive a gap between dialog systems that are designed to convey a
certain message in a goal oriented fashion and dialog systems that generate
chit-chat. Chit-chat can be generated rather freely because any topical
response is valid, where as a certain degree of factual correctness is to be
expected from a goal oriented system.

In this paper, we seek to bring the two lines of dialog generation research
closer together. We implement a system that can generate topical responses (in
the sense of chit-chat) with a fixed emotional content. Thus the goal of the
system is to convey a desired emotion in its response, no matter what the
actual textual content ended up being. This means that part of the semantics
of the out is fixed, while a part as to how to contextually adapt the
emotional content. is still up to computational creativity.

We base our experiments on a recently published dialog dataset that contains
sentiment annotations111https://zenodo.org/record/6990638. The dataset is
based on a video game called Fallout New Vegas and it has dialog where each
line of the dialog is annotated as containing one of the following emotions:
anger, disgust, fear, happiness, neutral, pain, sadness or surprise. This
dataset makes for an optimal training data for the task we seek to solve.

| Anger | Disgust | Fear | Happy | Neutral | Pained | Sad | Surprised
---|---|---|---|---|---|---|---|---
Sentences | 3 335 | 932 | 1 620 | 4 029 | 8 802 | 994 | 1 055 | 1 649
Table 1: Number of sentences per emotion Prompt | Response | Emotion
---|---|---
I hear you’ve been causing trouble. | Oh yeah? Fuck off, asshole. | Anger
I was hoping you’d be that stupid. | What? Hey, guys! Help me here! | Surprise
I’ve dealt with those newcomers. | I take care of those who help with that. Here, you earned it | Happiness
Table 2: Examples of the training data

##  2 Related Work

In terms of computational creativity and natural language generation, there
are several papers out there that present work conducted on a variety of
different creative language generation tasks such as poem generation Hegade et
al. (2021); Hämäläinen et al. (2022), humor generation Weller et al. (2020);
Alnajjar and Hämäläinen (2021), news generation Shu et al. (2021); Koppatz et
al. (2022) and story generation Vicente et al. (2018); Concepción et al.
(2019). In this section, we will take a closer look at the work conducted on
dialog generation.

Xie and Pu (2021) present work on generating empathetic dialog. Their model
deals with the following categories of empathetic intent: questioning,
agreeing, acknowledging, sympathizing, encouraging, consoling, suggesting and
wishing. They base their model on the transformer architecture and they use
RoBERTa for input encoding. In addition, they train a classifier that predicts
the salient empathetic intent.

Dialog generation has also been tackled in a context-controlled and topic-
aware manner Ling et al. (2021). Their model consists of four parts: a
hierarchical context encoder, a contex dependent topic representation module,
a context guided topic transition module and a joint attention based response
decoder.

Chen et al. (2021) present their work on multi-turn dialog generation. They
use a cross-hierarchical encoder that encodes a sentence for an answer
selector model, after this a response generator model is used to generate the
final output. The initial encoding is done by a transformer based model while
the final genration is done by an LSTM model.

Dialog adaptation has been studied before in the context of video games
Hämäläinen and Alnajjar (2019). The authors use an LSTM model to paraphrase
the syntax of existing dialog to introduce diversity and a word2vec model to
adapt the meaning of the sentence towards a desired player attribute.

##  3 Creativity and Emotion

There are several takes on creativity in a computational setting. In this
section, we cover some of these theoretical ways of understanding
computational creativity. Theoretical foundation has been, for a long time, at
the very core of computational creativity to combat systems that do mere
generation. That is generation for the sake of outputting something by any
means necessary.

A computationally creative system should exhibit skill, imagination and
appreciation according to Colton (2008). He argues that all of these three
components are a strict requirement for creativity to exist in a system. Skill
refers to the system’s capability of producing a creative artifact whereas
appreciation means that the system should also know why its creation is good.
Imagination requires the system to be capable of garnering a lot of diverse
output for one input.

Creativity can also be modeled through the FACE theory Colton et al. (2011).
This theory states creativity comes from the interplay between framing,
aesthetics, concept and expression. Expressions are the creative output
produced by the system. The system itself is called concept. Aesthetics is
similar to appreciation in the previous theory; it means that the system
should be able to appreciate the creative value of its output. Framing
highlights the fact that creativity does not take place in a vacuum but is
presented in a context. In our case, framing would be the entire dialog
between a human user and the machine.

Boden (1998) identifies three types of creativity; exploratory creativity,
transformational creativity and combinatory creativity. In combinatory
creativity, a system forms new artifacts by combining old ones in novel ways.
In exploratory creativity, a system is conducing a search in a conceptual
space discovering new creative artifacts. A system that can achieve
transformational creativity can change its search space.

A system is considered to be autonomously creative if it can change its own
standards without being explicitly told to do so Jennings (2010). The change
cannot occur at random either because a simple random change at random
intervals would otherwise be enough to satisfy the criterion.

Emotion is considered as a higher level cognitive phenomenon than a feeling
(see Shouse 2005). Feelings are seen to be universally felt in a similar
fashion as a response to some external or internal stimulus. Emotions, on the
other hand, are culturally and socially represented and their existence in all
cultures in a similar way is not a given thing (see Lim 2016).

One way of seeing emotions is that they rely on affect (see Russell 2003),
which is a state our mind is continuously in. An affect can move across two
axes: positive-negative and arousal-relaxed. The affect we feel is
contextually resolved to a higher level emotion based on the context we find
ourselves in. For instance, a high arousal and negativity could be interpreted
either as anger or disgust among others.

Ekman (1992) has identified six basic emotions: anger, disgust, fear,
happiness, sadness and surprise. These emotions are also present in our
dataset, which makes this theory optimal to build upon. The basic emotions are
considered to be universal across cultures based on studies conducted on
facial expressions.

Prompt | Emotion | Output
---|---|---
Doesn’t look like the talkative type. | Anger | What the fuck are you looking at?
Wait a minute! I ain’t done nothin’ to you! | Fear | I don’t want any trouble, okay?
Asked you a question. How many? | Sad |  | I’m sorry, but I don’t have the time
---
to talk right now.
| Need coin, fast? Those blasted nonhumans at the
---
Vivaldi Bank say your credit’s no good? Come see me!
Surprise | I’ve never heard of that before.
Table 3: Examples of the results

##  4 Data

Our dataset consists of dialog in English where each line is annotated with an
emotion label. Table 1 shows the data size and how many sentences were there
in the dataset for each emotion. The dialog is from a video game called
Fallout New Vegas222https://fallout.bethesda.net/en/games/fallout-new-vegas by
Obsidian Entertainment. The game is set in a post-apocalyptic world inhabitet
by creatures that were born as a result of nuclear radiation such as mutants,
ghouls and feral ghouls in addition to humans.

The game dialog consists mostly of player prompts that result in a response by
an NPC (non-player character). In Table 2 we can see some examples of the
dialog. In our case, we train our model by using the prompt and the emotion
label as an input to predict the response.

Fallout New Vegas is a relatively large video game because it is an open world
RPG (role-playing game) where the player can roam freely from one place to
another. As a result of this, the game has a variety of different scripted
characters which means that there is no bias in terms of having the same
characters speaking with each other all the time. There is, however, a bias in
the topic of conversations given that they take place in a fictional world.
Many of the dialogs deal with fictional places, characters, items and so on.

##  5 Dialog Generation

In this section, we outline our approach to emotionally conditioned dialog
generation. We conduct our experiment using Transformers Wolf et al. (2020)
and Datasets Lhoest et al. (2021) Python libraries. We base our model on a
pretrained model called DialoGPT
medium333https://huggingface.co/microsoft/DialoGPT-medium.

DialoGPT Zhang et al. (2020) is based on the GPT-2 Radford et al. (2019)
architecture, which in turn is based on the generic transformer language model
Vaswani et al. (2017). The transformer model leverages a stack of masked
multi-head self-attention layers to train on large datasets. DialoGPT employs
a maximum mutual information scoring function Li et al. (2015); Zhang et al.
(2018) during the training phase optimizing the reward with a policy gradient
Williams (1992) with a sample averaged baseline Zhang et al. (2018).

We changed the format of the data to be as “EMOTION1: SENTENCE1. EMOTION2:
SENTENCE2 [EOS]“, where EMOTIONn indicates the emotion of the $n$th sentence
in the conversation. This way, the model is exposed to emotional knowledge
regarding the preceding locution of the response, as well. We use 90% of the
data for training and 10% for validation. We train the model for 5 epochs.

##  6 Results and Evaluation

In order to evaluate our model, we sample dialog from an unrelated open-world
RPG called The Witcher 3: Wild Hunt444https://www.thewitcher.com/en/witcher3
by CD Projekt. The dialog is extracted using a w3strings decoder
tool555https://www.nexusmods.com/witcher3/mods/1055/. By using dialog from a
different video game, we can see whether our model works in a different
domain. The Witcher 3 is set to a medieval fantasy world inhabited by magical
creatures. Some of the generated output can be seen in Table 3.

We sample randomly 15 sentences for each emotion, and produce emotion
conditioned replies. In the input, “EMOTION1:“ was not present as the dialogs
in The Witcher did not include emotion labels. These are then evaluated by
crowd-workers on a crowd-sourcing platform called Appen666https://appen.com/.
For every output we evaluate, we show the input sentence and the output to the
user and ask them if the output expresses the desired emotion as a simple
yes/no question. In case, the user selected ”yes”, we also asked how strongly
the line expressed the desired emotion on a 5-point scale (the stronger, the
higher), the scale is 0-based from 0 to 4. Each generated output was evaluated
by 5 different judges.

Figure 1: Percentage of ”yes” answers per emotion

All in all, the judges agreement with the emotion label was 0.6. In Figure 1,
we can see the judges’ agreement per emotion label. The worst performing
emotion is pained and the best performing ones are neutral, fear and disgust.
Sad, surprise and anger are all rated expressing the emotion over half of the
time. These results are based on the aggregated results provided by Appen.

Figure 2: Strength of each emotion

Figure 2 shows how strong the judges rated each emotion to be. 0 indicates not
strong at all and 4 very strong. As we can see, all emotions achieve above
mid-level performance. Anger, fear and disgust being the most strongly
expressed emotions. All in all, based on the results, it seems that the model
is capable of producing emotional responses in a suitably strong manner.

##  7 Conclusions

Emotion conditioned dialog generation remains a challenge. In this paper, we
have shown the viability of using a DialoGPT model for this task. The results
are promising especially given that the evaluation was run on out-of-the-
domain data.

The video game that was used in this paper also comes with audio files for
each line of dialog. This is an interesting direction for future research
because we could generate emotion conditioned dialog with audio as well. The
same sentence can be said with multiple different tones and intonations to
express different emotions. Building a system that can express emotion in
generated speech as well would have great application potentials in speech
oriented dialog systems.

## Acknowledgments

This work was partially financed by the Society of Swedish Literature in
Finland with funding from Enhancing Conversational AI with Computational
Creativity, and by the Ella and Georg Ehrnrooth Foundation for Modelling
Conversational Artificial Intelligence with Intent and Creativity.

## References

  * Alnajjar and Hämäläinen (2021) Khalid Alnajjar and Mika Hämäläinen. 2021.  When a computer cracks a joke: Automated generation of humorous headlines.  In _Proceedings of the 12th International Conference on Computational Creativity (ICCC 2021)_. Association for Computational Creativity.
  * Boden (1998) Margaret A Boden. 1998.  Creativity and artificial intelligence.  _Artificial intelligence_ , 103(1-2):347–356.
  * Chen et al. (2021) Xiuying Chen, Zhi Cui, Jiayi Zhang, Chen Wei, Jianwei Cui, Bin Wang, Dongyan Zhao, and Rui Yan. 2021.  Reasoning in dialog: Improving response generation by context reading comprehension.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pages 12683–12691.
  * Colton (2008) Simon Colton. 2008.  Creativity versus the perception of creativity in computational systems.  In _AAAI spring symposium: creative intelligent systems_ , volume 8, page 7. Palo Alto, CA.
  * Colton et al. (2011) Simon Colton, John William Charnley, and Alison Pease. 2011.  Computational creativity theory: The face and idea descriptive models.  In _ICCC_ , pages 90–95. Mexico City.
  * Concepción et al. (2019) Eugenio Concepción, Pablo Gervás, and Gonzalo Méndez. 2019.  Evolving the ines story generation system: From single to multiple plot lines.  In _ICCC_ , pages 220–227.
  * Ekman (1992) P Ekman. 1992.  Are there basic emotions?  _Psychological review_ , 99(3):550–553.
  * Hämäläinen and Alnajjar (2019) Mika Hämäläinen and Khalid Alnajjar. 2019.  Creative contextual dialog adaptation in an open world rpg.  In _Proceedings of the 14th International Conference on the Foundations of Digital Games_ , pages 1–7.
  * Hämäläinen et al. (2022) Mika Hämäläinen, Khalid Alnajjar, and Thierry Poibeau. 2022.  Modern french poetry generation with roberta and gpt-2.  In _13th International Conference on Computational Creativity (ICCC) 2022_.
  * Hegade et al. (2021) Prakash Hegade, Rajaram M Joshi, Vibha G Hegde, Tejaswini Kale, and Srushti Basavaraddi. 2021.  Po-miner: A web mining poem generator and its security model.  _SN Computer Science_ , 2(5):1–14.
  * Jennings (2010) Kyle E Jennings. 2010.  Developing creativity: Artificial barriers in artificial intelligence.  _Minds and Machines_ , 20(4):489–501.
  * Koppatz et al. (2022) Maximilian Koppatz, Khalid Alnajjar, Mika Hämäläinen, and Thierry Poibeau. 2022.  Automatic generation of factual news headlines in finnish.  In _15th International Natural Language Generation Conference (INLG)_.
  * Lhoest et al. (2021) Quentin Lhoest, Albert Villanova del Moral, Yacine Jernite, Abhishek Thakur, Patrick von Platen, Suraj Patil, Julien Chaumond, Mariama Drame, Julien Plu, Lewis Tunstall, Joe Davison, Mario Šaško, Gunjan Chhablani, Bhavitvya Malik, Simon Brandeis, Teven Le Scao, Victor Sanh, Canwen Xu, Nicolas Patry, Angelina McMillan-Major, Philipp Schmid, Sylvain Gugger, Clément Delangue, Théo Matussière, Lysandre Debut, Stas Bekman, Pierric Cistac, Thibault Goehringer, Victor Mustar, François Lagunas, Alexander Rush, and Thomas Wolf. 2021.  Datasets: A community library for natural language processing.  In _Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 175–184, Online and Punta Cana, Dominican Republic. Association for Computational Linguistics.
  * Li et al. (2015) Jiwei Li, Michel Galley, Chris Brockett, Jianfeng Gao, and Bill Dolan. 2015.  A diversity-promoting objective function for neural conversation models.  _arXiv preprint arXiv:1510.03055_.
  * Lim (2016) Nangyeon Lim. 2016.  Cultural differences in emotion: differences in emotional arousal level between the east and the west.  _Integrative medicine research_ , 5(2):105–109.
  * Ling et al. (2021) Yanxiang Ling, Fei Cai, Xuejun Hu, Jun Liu, Wanyu Chen, and Honghui Chen. 2021.  Context-controlled topic-aware neural response generation for open-domain dialog systems.  _Information Processing & Management_, 58(1):102392.
  * Radford et al. (2019) Alec Radford, Jeffrey Wu, Rewon Child, David Luan, Dario Amodei, Ilya Sutskever, et al. 2019.  Language models are unsupervised multitask learners.  _OpenAI blog_ , 1(8):9.
  * Russell (2003) James A Russell. 2003.  Core affect and the psychological construction of emotion.  _Psychological review_ , 110(1):145.
  * Shouse (2005) Eric Shouse. 2005.  Feeling, emotion, affect.  _M/c journal_ , 8(6).
  * Shu et al. (2021) Kai Shu, Yichuan Li, Kaize Ding, and Huan Liu. 2021.  Fact-enhanced synthetic news generation.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 35, pages 13825–13833.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  _Advances in neural information processing systems_ , 30.
  * Vicente et al. (2018) Marta Vicente, Cristina Barros, and Elena Lloret. 2018.  Statistical language modelling for automatic story generation.  _Journal of Intelligent & Fuzzy Systems_, 34(5):3069–3079.
  * Weller et al. (2020) Orion Weller, Nancy Fulda, and Kevin Seppi. 2020.  Can humor prediction datasets be used for humor generation? humorous headline generation via style transfer.  In _Proceedings of the Second Workshop on Figurative Language Processing_ , pages 186–191.
  * Williams (1992) Ronald J Williams. 1992.  Simple statistical gradient-following algorithms for connectionist reinforcement learning.  _Machine learning_ , 8(3):229–256.
  * Wolf et al. (2020) Thomas Wolf, Lysandre Debut, Victor Sanh, Julien Chaumond, Clement Delangue, Anthony Moi, Pierric Cistac, Tim Rault, Remi Louf, Morgan Funtowicz, Joe Davison, Sam Shleifer, Patrick von Platen, Clara Ma, Yacine Jernite, Julien Plu, Canwen Xu, Teven Le Scao, Sylvain Gugger, Mariama Drame, Quentin Lhoest, and Alexander Rush. 2020.  Transformers: State-of-the-art natural language processing.  In _Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing: System Demonstrations_ , pages 38–45, Online. Association for Computational Linguistics.
  * Xie and Pu (2021) Yubo Xie and Pearl Pu. 2021.  Empathetic dialog generation with fine-grained intents.  In _Proceedings of the 25th Conference on Computational Natural Language Learning_ , pages 133–147.
  * Zhang et al. (2018) Yizhe Zhang, Michel Galley, Jianfeng Gao, Zhe Gan, Xiujun Li, Chris Brockett, and Bill Dolan. 2018.  Generating informative and diverse conversational responses via adversarial information maximization.  _Advances in Neural Information Processing Systems_ , 31.
  * Zhang et al. (2020) Yizhe Zhang, Siqi Sun, Michel Galley, Yen-Chun Chen, Chris Brockett, Xiang Gao, Jianfeng Gao, Jingjing Liu, and William B Dolan. 2020.  Dialogpt: Large-scale generative pre-training for conversational response generation.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics: System Demonstrations_ , pages 270–278.


# Sonic Gradient Index Lens for Aqueous Applications

Theodore P. Martin  Acoustics Division, Naval Research Laboratory, Washington,
DC 20375, USA  Michael Nicholas  Acoustics Division, Naval Research
Laboratory, Washington, DC 20375, USA  Gregory J. Orris  Acoustics Division,
Naval Research Laboratory, Washington, DC 20375, USA  Liang-Wu Cai  Department
of Mechanical and Nuclear Engineering, Kansas State University, Manhattan,
Kansas 66506, USA  Daniel Torrent  Wave Phenomena Group, Department of
Electronic Engineering, Universidad Politecnica de Valencia, C/ Camino de Vera
s7n, E-46022 Valencia, Spain  José Sánchez-Dehesa  Wave Phenomena Group,
Department of Electronic Engineering, Universidad Politecnica de Valencia, C/
Camino de Vera s7n, E-46022 Valencia, Spain

###### Abstract

We study the acoustic scattering properties of a phononic crystal designed to
behave as a gradient index lens in water, both experimentally and
theoretically. The gradient index lens is designed using a square lattice of
stainless-steel cylinders based on a multiple scattering approach in the
homogenization limit. We experimentally demonstrate that the lens follows the
graded index equations derived for optics by mapping the pressure intensity
generated from a spherical source at 20 kHz. We find good agreement between
the experimental result and theoretical modeling based on multiple scattering
theory.

###### pacs:

43.20.Fn, 43.58.Ls,43.20.Dk

††preprint: NRL-TMartin-v1.4

Composed of ordered arrays of scatterers similar to atoms in a conventional
solid, phononic crystals (PnC) are a class of metamaterial designed to control
acoustic wave propagation in a medium. PnCs have been proposed for a broad
range of applications in wave acoustics, with acoustic lensing Cervera  _et
al._ (2002); Hu and Chan (2005); Yang  _et al._ (2004); Torrent and Sánchez-
Dehesa (2007); Cai  _et al._ (2007); Li  _et al._ (2009); Qiu  _et al._
(2005); Deng  _et al._ (2009); Zhang  _et al._ (2009); Lin  _et al._ (2009);
Ke  _et al._ (2007) featuring prominently in the literature due in part to the
ease with which focusing can be achieved by altering a crystal’s gemoetric
shape Cervera  _et al._ (2002); Hu and Chan (2005) or compositional structure.
Cai  _et al._ (2007); Torrent and Sánchez-Dehesa (2007); Li  _et al._ (2009)
Although negative index lenses have received much attention due to their
potential for near-field imaging, Qiu  _et al._ (2005); Deng  _et al._ (2009);
Zhang  _et al._ (2009) some positive index solutions such as the acoustic
analogue of the optical graded index lens Torrent and Sánchez-Dehesa (2007);
Lin  _et al._ (2009) have not yet been explored experimentally. In addition,
the majority of PnC experiments have been performed in air, Cervera  _et al._
(2002); Li  _et al._ (2009); Torrent  _et al._ (2006); Ke  _et al._ (2007)
where the large density contrast with respect to the constituent scatterers in
the crystal (typically metals) allows the scatterers to be treated as rigid.
We demonstrate below that despite the physical limitation in impedance
contrasts between an aqueous medium and the scattering elements, it is
possible to design a PnC that behaves as an ideal graded index lens (GIL) in
water based on a fully elastic multiple scattering theory (MST). Torrent  _et
al._ (2006); Torrent and Sánchez-Dehesa (2007); Krokhin  _et al._ (2003)

Figure 1(a) shows a plan view schematic of the GIL design. The axes of Fig.
1(a) and throughout the paper are oriented with the lens center at position
$(x,y)=(0,0)$. The GIL is made up of 75 stainless steel cylinders (T-316) that
are 75 cm in length and arranged in a square lattice with spacing $a=1.8$ cm
and dimensions $5a\times 15a$. Figure 1(b) plots the cylinder radii $R(y)$,
which are stepped toward zero at each successive layer above and below the
central axis ($y=0$) of the GIL. In the homogenization limit (propagation
wavelength $\lambda\gtrsim 4a$) each stratified layer in Fig. 1(a) can be
treated as an effective medium. MST Torrent and Sánchez-Dehesa (2007) is used
to calculate each layer’s effective sound speed $c_{eff}$, which is inversely
proportional to the filling fraction of the cylinders. The layers will have an
effective refractive index $n_{eff}=c_{b}/c_{eff}$ ($c_{b}=1470$ m/s is the
sound speed in water) that is maximal at the center of the GIL and decreases
to that of water at the edges. Our choice of $R(y)$ in Fig. 1(b) produces a
graded $n_{eff}$ that obeys the same relation as an optical GIL, Smith (2000)
$n_{eff}=n_{0}(1-\alpha^{2}y^{2})^{1/2}$, where $n_{0}$ is the refractive
index at the central layer. Our design results in $n_{0}=1.2$ and
$\alpha=0.04$ cm-1.

Figure 1(c) shows an image of the GIL submerged in a $6\times 6\times 4$ m3
isolation tank. The cylinders are mounted between reinforced Plexiglas plates
to provide stability. Acoustic waves are produced by a 10 cm-diameter
spherical source at 20 kHz ($\lambda\simeq 4a$). The wave propagation is
measured in the time-domain using monitoring hydrophones at a sampling rate of
1 MHz. Hydrophones are mounted to the source and onto a translational 3-axis
Velmex VXM® positioning system. The transmission intensity is measured by
averaging over a 10-cycle pulse from the source; this pulse is long enough to
approximate a continuous wave measurement, while being short enough to prevent
contamination from reflections off the surfaces of the tank. We have
experimentally verified that the intensity $P_{0}$ produced by the source in
the absence of the GIL drops radially in proportion to $1/r^{2}$.

Figure 2(a) shows the normalized pressure amplitude $P/P_{0}$ measured after
transmission through the GIL on the side opposite the source ($x>0$). The
source is located at $(x,y)=(-196,0)$ cm, and both the source and the
translational hydrophone are positioned in the plane bisecting the axial
center of the cylinders. The GIL is shown schematically to scale and at its
proper location in each figure throughout the paper. The data in Fig. 2(a) is
measured 2.144 ms after the initial cycle began to leave the source. This time
gives a snapshot when the pulse is centered on an enhancement in signal
amplitude observed in the vicinity of $x\simeq 80$ cm.

Figure 2(b) shows the normalized intensity averaged over the 10-cycle pulse
and obtained from the same data set shown in Fig. 2(a). As in Fig. 2(a), a
clear focusing peak is observed centered close to $x\simeq 80$ cm. In Fig.
2(c) we show a two-dimensional MST calculation Torrent  _et al._ (2006);
Torrent and Sánchez-Dehesa (2007) of the total pressure intensity (incident
plus scattered) derived by placing a continuous-wave cylindrical source at
$(x,y)=(-196,0)$ cm. The calculation assumes the cylinders to be a penetrable
elastic. Torrent and Sánchez-Dehesa (2007); Krokhin  _et al._ (2003) As with
the experimental data, the simulated pressure intensity is normalized to that
of the source in the absence of the GIL. The source amplitude is a Hankel
function $P_{0}=H^{(1)}_{0}(kr)$ with wavevector $k=\pi/2a$. The MST
simulation also shows a clear focusing peak, but with two important
differences: (1) the measured intensity is $\sim 2$ times larger than the
simulation, and (2) the simulated focusing peak is slightly farther from the
GIL and decays more slowly.

To quantify whether our GIL design behaves as an ideal lens, in Fig. 3(a) we
present measurements of the focusing peak along the central axis of the lens
($y=0$) for different source positions $d_{s}$. For each $d_{s}$, a large-
amplitude peak is observed above $x\gtrsim 60$ cm, while smaller peaks are
also observed closer to the GIL. As the source is moved closer to the GIL, the
large-amplitude peak moves away in qualitative agreement with the expected
behavior of a lens. In Fig. 3(b) we show MST calculations along $y=0$ for
source positions similar to those in the experiment. On initial inspection it
appears that the theory shows slowly decaying focusing peaks that change very
little with $d_{s}$. However, expansion of the region around the focusing
peaks [Fig. 3(b) inset] reveals that the peak positions move away in a manner
similar to the experiment.

We now analyze the experimental data in Fig. 3(a) above $x>62$ cm to determine
whether the focusing positions in this region follow the ideal lens equations.
For an ideal lens, the focusing peak positions $d_{p}$ should scale with
$d_{s}$ as $1/d_{p}=1/f-1/d_{s}$. The focal length $f$ of a GIL can be
approximated as, Smith (2000)

$\displaystyle f\approx\frac{1}{n_{0}\alpha\sin\alpha t}$ (1)

where $t=5a$ is the thickness of the GIL. Equation (1) gives an estimate of
$f=58.9$ cm using values of $n_{0}$ and $\alpha$ calculated in the effective
medium approximation. 111Although we are technically measuring the _back focal
length_ of the GIL, $bfl=f\cos\alpha t$, the difference between our estimate
of $f$ and the sum $bfl+t/2$ is less than $1$ cm, which is less than both our
experimental error and our data resolution ($=1$ cm). Thus we ignore this
subtlety and compare to $f$ directly. See Ref. 14 for details.

A close inspection of Fig. 3(a) reveals that the data above $x>62$ cm is
actually composed of two superimposed peaks that both move to larger $x$ as a
function of $d_{s}$. Figure 3(c) shows two examples of a double-gaussian fit
to the data in this region using a standard unconstrained, nonlinear
optimization routine. The gaussians resulting from the fit are shown
individually (blue and red) in addition to the combined fit (black). Although
the fit equation
$\gamma_{1}e^{-\beta_{1}(x-d_{p1})^{2}}+\gamma_{2}e^{-\beta_{2}(x-d_{p2})^{2}}$
contains six free parameters, the purpose of the fit is to obtain an estimate
of the peak positions $d_{p1,2}$ and their relative amplitudes $\gamma_{1,2}$.
Fig. 3(c) demonstrates that the data is well described by the double-gaussian,
with a low-amplitude peak (Peak 1) closer to the GIL and a larger-amplitude
peak (Peak 2) farther away. In both cases the amplitude of Peak 2 is about
three times larger than Peak 1, suggesting that Peak 2 is the main focusing
peak of the GIL. The relative amplitudes of Peaks 1 and 2 are observed to
follow the same qualitative behavior for all the source positions in Fig.
3(a).

Figure 3(d) plots the inverse positions $1/d_{p1,2}$ extracted from the
gaussian fits as a function of $1/d_{s}$. An ideal lens will produce a linear
trend with a slope of $-1$ and an intercept of $1/f$. Although the trends for
both peaks are linear and have intercepts that yield similar focal lengths,
the slope of Peak 1 is much less than that of an ideal lens. However, the
trend for Peak 2 results in a slope of $-1$ and its focal length $f=59.3\pm
1.5$ cm agrees with the estimate of $f$ calculated using Eqn. (1). The dashed
line in Fig. 3(d) plots the peak locations obtained from the MST calculations
in Fig. 3(b). The theory produces a slope of $-1$ and focal length $f=61.1$ cm
that closely match both the measured data and the ideal lens equations.

We propose that Peak 1 and the other low-amplitude peaks in Fig. 3(a) are the
result of constructive interference between waves scattered from the support
structure of the lens. Low-amplitude, circular interference fringes can be
observed in Fig. 2(a) emanating from above and below the plot area centered at
$x\simeq 25$ cm. These fringes are the result of scattering off of stabilizing
pillars at the corners of the GIL support structure. While averaging over a
few initial pulse cycles will reduce the interference, a small number of
cycles gives a poor approximation to a continuous wave measurement and limits
the number of multiple scattering events that contribute to the focusing peak.
Therefore we have chosen to average over many cycles and rely on the gaussian
fitting routine to remove the spurious interference.

Figure 4 demonstrates that our GIL design acts as a lens with the source off
the central axis. Figure 4(a) shows $(P/P_{0})^{2}$ measured with the
spherical source located at a $14.7\,^{\circ}$ angle with respect to the
origin. Figure 4(b) shows the MST calculation for the same source location.
Thin white lines in Figs. 4(a,b) indicate a $14.7\,^{\circ}$ angle with
respect to the $x$-axis and extend to the expected focusing positions of an
ideal lens with $f\approx 60$ cm. Both the measured data and the MST
calculation demonstrate a strong focusing peak at the expected location.
Interference fringes from the support pillar can also be observed superimposed
on the focusing peak in Fig. 4(a).

In summary, we have designed and constructed a gradient index lens that
operates in water at sonic frequencies. Our transmission measurements
demonstrate that our GIL design focuses as an ideal lens based on the optical
GIL equations. Our measurements are also consistent with the focusing
positions obtained from two-dimensional models using multiple scattering
theory. We emphasize that our GIL behaves as an ideal lens at the limit of
homogenization ($\lambda\simeq 4a$) and with a thickness on the order of a
wavelength ($t=5\lambda/4$). Such performance at the limit of homogenization
theory demonstrates the versatility of phononic crystals designed using
multiple scattering theory.

This work was supported by the U.S. Office of Naval Research.

## References

  * Cervera  _et al._ (2002) F. Cervera, L. Sanchis, J. V. Sánchez-Pérez, R. Martínez-Sala, C. Rubio, F. Meseguer, C. López, D. Caballero, and J. Sánchez-Dehesa, Phys. Rev. Lett., 88, 023902 (2002).
  * Hu and Chan (2005) X. Hu and C. T. Chan, Phys. Rev. Lett., 95, 154501 (2005).
  * Yang  _et al._ (2004) S. Yang, J. H. Page, Z. Liu, M. L. Cowan, C. T. Chan, and P. Sheng, Phys. Rev. Lett., 93, 024301 (2004).
  * Torrent and Sánchez-Dehesa (2007) D. Torrent and J. Sánchez-Dehesa, New J. Phys., 9, 323 (2007).
  * Cai  _et al._ (2007) F. Cai, F. Liu, Z. He, and Z. Liu, Appl. Phys. Lett., 91, 203515 (2007).
  * Li  _et al._ (2009) J. Li, L. Fok, X. Yin, G. Bartal, and X. Zhang, Nature Mater., 8, 931 (2009).
  * Qiu  _et al._ (2005) C. Qiu, X. Zhang, and Z. Liu, Phys. Rev. B, 71, 054302 (2005).
  * Deng  _et al._ (2009) K. Deng, Y. Ding, Z. He, H. Zhao, J. Shi, and Z. Liu, J. Phys. D: Appl. Phys., 42, 185505 (2009).
  * Zhang  _et al._ (2009) S. Zhang, L. Yin, and N. Fang, Phys. Rev. Lett., 102, 194301 (2009).
  * Lin  _et al._ (2009) S.-C. S. Lin, T. J. Huang, J.-H. Sun, and T.-T. Wu, Phys. Rev. B, 79, 094302 (2009).
  * Ke  _et al._ (2007) M. Ke, Z. Liu, P. Pang, C. Qiu, D. Zhao, S. Peng, and J. Shi, Appl. Phys. Lett., 90, 083509 (2007).
  * Torrent  _et al._ (2006) D. Torrent, A. Hakansson, F. Cervera, and J. Sánchez-Dehesa, Phys. Rev. Lett., 96, 204302 (2006).
  * Krokhin  _et al._ (2003) A. A. Krokhin, J. Arriaga, and L. N. Gumen, Phys. Rev. Lett., 91, 264302 (2003).
  * Smith (2000) W. J. Smith,  _Modern Optical Engineering_ (McGraw-Hill, Inc., 2000).
  * Note (1) Although we are technically measuring the _back focal length_ of the GIL, $bfl=f\mathop{cos}\nolimits\alpha t$, the difference between our estimate of $f$ and the sum $bfl+t/2$ is less than $1$ cm, which is less than both our experimental error and our data resolution ($=1$ cm). Thus we ignore this subtlety and compare to $f$ directly. See Ref. 14 for details.

Figure 1: (a) Plan schematic of the gradient index lens. (b) Cylinder radius
$R$ plotted vs position along the $y$-axis. (c) Digital photograph of the GIL
in the isolation tank.

Figure 2:  (a) Normalized pressure amplitude $P/P_{0}$ plotted vs $x$ and $y$,
measured 2.144 ms after the initial pulse leaves the source. (b) Measured,
normalized pressure intensity $(P/P_{0})^{2}$ plotted vs $x$ and $y$ after
averaging over a 10-cycle pulse. (c) Normalized pressure intensity
$(P/P_{0})^{2}$ calculated using MST. The focusing peak maximum is marked by a
$+$.

Figure 3:  (a) Measured, normalized pressure intensity vs $x$ for source
positions $d_{s}=196$, $185.8$, $175.7$, $165.5$, $155.4$, $145.2$, $135$,
$124.9$, and $114.7$ cm. (b) Normalized pressure intensity calculated using
MST for source positions $d_{s}=109a$, $103a$, $98a$, $92a$, $86a$, $81a$, and
$75a$. Inset: expanded region of the $y$-axis showing the focusing peak
positions. (c) Two focusing peaks from panel (a) are replotted as circles
(upper region offset for clarity). Black lines indicate a double-gaussian fit,
with blue (Peak 1) and red (Peak 2) lines showing the component gaussians
individually. (d) Inverse peak positions $1/d_{p1,2}$ plotted vs inverse
source positions $1/d_{s}$. Blue and red lines are fits to the trends of Peaks
1 and 2 respectively. The dashed line plots the peak positions of the
simulated data in panel (b).

Figure 4:  (a) Measured, normalized pressure intensity $(P/P_{0})^{2}$ plotted
vs $x$ and $y$ with the source at a $14.7\,^{\circ}$ angle with respect to the
$x$-axis. (b) Normalized pressure intensity $(P/P_{0})^{2}$ calculated using
multiple scattering theory with the source at a $15\,^{\circ}$ angle. White
lines indicate the off-axis angle; positions of the expected focusing peaks
are marked with a $+$.


[1]Richard Tattersall

1]Joint Quantum Centre (JQC) Durham-Newcastle, School of Mathematics,
Statistics and Physics, Newcastle University, Newcastle upon Tyne, NE1 7RU,
United Kingdom

# Non-equilibrium Dynamics of Vortices in Two-Dimensional Quantum Gases:
Determining the Dynamical Scaling Region Using the Mahalanobis Distance

<EMAIL_ADDRESS>Andrew Baggaley  Thomas Billam  [

###### Abstract

When a two-dimensional system undergoes a rapid quench from a disordered to an
ordered phase, it does not order instantly but instead relaxes towards
equilibrium over time. During this relaxation, the dynamical scaling
hypothesis predicts that the length scale of ordered regions increases, with
later patterns statistically similar to earlier ones except for a change in
global scale. Quantum gases are one of many systems in which such out of
equilibrium behaviour is predicted to occur. Here we present a method for
systematically testing when dynamical scaling is taking place, by quantifying
the similarity of the rescaled two-point correlation function over time using
the Mahalanobis distance. Data on the velocity field of a quantum fluid,
generated from point vortex simulations, are used to illustrate the
application of this method.

###### keywords:

universal, dynamical scaling, vortex dynamics, non-equilibrium, quantum fluid

##  1 Introduction

People search for evidence of dynamical scaling in systems such as the Ising
model of magnetic spins, binary alloys [1, 2] and Bose-Einstein Condensates
(BECs) [3]. A key prediction of the dynamical scaling hypothesis is that the
correlation function, $G(r,t)$, between two points separated by a distance
$r$, changes in time, $t$, only by a rescaling in length scale by $L_{c}(t)$,
where $L_{c}(t)$ is a correlation length that depends on time only.

$G(r,t)=G(r/L_{c}(t))\,.$ (1)

To test if this prediction is met, the correlation functions at different
times are rescaled by $L_{c}(t)$ and plotted on the same axes. If they
collapse onto a single curve, this is evidence of dynamical scaling.
Typically, this collapse is checked visually; however, a more rigorous,
algorithmic test is preferable in order to determine the times at which
dynamical scaling applies. In this paper we outline an approach that makes use
of the Mahalanobis distance [4] to quantify the deviation of a curve from the
mean, eliminate outliers and identify the scaling region.

When dynamical scaling holds, the dependence of correlation length on time is
predicted to be a power law of the form:

$L_{c}(t)\sim t^{1/z}\,,$ (2)

where $z$ is the dynamical critical exponent and is expected to have a value
of two for a two-dimensional system with a non-conserved order parameter [1,
2]. A recent numerical study of coarsening dynamics, following the rapid
quench of a Bose gas into the ordered phase, found clear evidence of dynamical
scaling [3]. The values of $z$ found were close to two but dependent upon
dissipation, parameterized by a dimensionless constant, $\gamma$. This study
used the projected Gross-Pitaevskii equation (PGPE) and the stochastic
projected Gross-Pitaevskii equation (SPGPE) to model the evolution of the
wavefunction $\psi(\mathbf{r},t)$, both in the absence and the presence of
dissipation. The authors observed growing regions of constant phase following
the quench, as quantum vortices and antivortices annihilated, as shown in Fig.
2. They found the correlation functions at different times using the
wavefunction and observed them to collapse, within a given scaling window,
onto a single curve once rescaled by $L_{c}(t)$. Their method for identifying
the scaling window was to judge by eye initially and later check that the
maximum discrepancy between rescaled correlation functions was sufficiently
small. An earlier study [5], simulating a similar system using a truncated
Wigner approach, also found evidence of universal scaling behaviour in the
time evolution of the single-particle momentum spectrum. They selected a
scaling window by inspecting the rescaled spectra and found a value of
$z=1.8(3)$ (avoiding a strongly anomalous non-thermal fixed point associated
with vortex clustering).

Figure 1: Evolution of the phase of the field $\psi$ of a BEC in the case of
zero dissipation, reproduced from [3] under CC-BY-4.0 international license
222https://creativecommons.org/licenses/by/4.0/. Distances are quoted in terms
of the healing length, $\xi=\hbar/\sqrt{m\mu}$, and times are given in terms
of $\tau=\hbar/\mu$. These are the characteristic length and time scales used
when solving the Gross-Pitaevskii equation, where $m$ is the atomic mass,
$\hbar$ is the reduced Planck constant and $\mu$ is the chemical potential.
The BEC was modelled using the projected Gross-Pitaevskii equation (PGPE) in a
doubly-periodic square box of width $\approx 363\xi$. The BEC was created in a
far from equilibrium configuration with a high density of quantised vortices
and antivortices. It then “relaxes”, in this case through annihilation of
vortex-antivortex pairs. (a), (b) and (c) correspond to times
$t\approx\\{200,2000,20000\\}\tau$ respectively. The open white squares show
the positions of vortices (positive circulation), the filled black squares
show the positions of antivortices (negative circulation) and the background
colour indicates the phase.

To illustrate our method to test for the collapse of rescaled correlation
functions, we also consider a two-dimensional quantum fluid. However, rather
than simulate the system using the Gross-Pitaevskii equation (GPE) [6, 7], we
use the point vortex model (PVM) instead [8, 9]. The PVM is often used to
simulate dynamics in two-dimensional quantum fluids [10, 11, 12] such as those
formed in ultracold atomic gases [13, 14]. Although vortices in atomic gases
are not point-like and the gas itself is compressible, the PVM can give
reasonably accurate predictions for vortex behaviour. It also benefits over
the more directly applicable GPE because only the coordinates of the vortices
need tracking and not the quantum mechanical wavefunction at every grid-point.
Consequently, it is computationally less expensive and makes simulations of
larger systems possible. The structure of this paper is as follows: Section 2
describes how data for the out of equilibrium dynamics of a two-dimensional
quantum fluid is generated using a dissipative PVM. Section 3 details how a
two-point correlation function is calculated and rescaled and how the scaling
region is determined using a method based on the Mahalanobis distance. Section
4 summarises the main findings and suggests some possibilities for future work
using this method.

##  2 Generating Data Using the Point Vortex Model (PVM)

###  2.1 The PVM on an Infinite Plane

The PVM applies to two-dimensional incompressible fluids with zero viscosity
and assumes that vorticity is confined to delta-function-like points [8]. Each
vortex is advected by its local fluid velocity, which is made up of
contributions from each of the other vortices. Consequently, the equations of
motion for vortex $i$, of circulation, $\kappa_{i}$, at position
$(x_{i},y_{i})$, in a system of $N_{v}$ point vortices are:

$\frac{dx_{i}}{dt}=-\frac{1}{2\pi}\sum_{j=1,j\neq
i}^{N_{v}}\frac{\kappa_{j}y_{ij}}{|\mathbf{r}_{ij}|^{2}}\,,\,\,\frac{dy_{i}}{dt}=\frac{1}{2\pi}\sum_{j=1,j\neq
i}^{N_{v}}\frac{\kappa_{j}x_{ij}}{|\mathbf{r}_{ij}|^{2}}\,.$ (3)

Where $\kappa_{j}$ and $(x_{j},y_{j})$ are the circulations and positions of
all the vortices other than $i$ and $x_{ij}=x_{i}-x_{j}$, $y_{ij}=y_{i}-y_{j}$
and $\mathbf{r}_{ij}=(x_{ij},y_{ij})$. This result may also be derived from
the Hamiltonian [9]:

$H=-\frac{1}{4\pi}\sum_{i=1}^{N_{v}}\sum_{j=1,j\neq
i}^{N_{v}}\kappa_{i}\kappa_{j}\ln{|\mathbf{r}_{ij}|}\,,$ (4)

in which the double sum adds contributions from all pairs of vortices.

In the case of quantum fluids, vortex circulations are quantised and take only
the values of $\pm h/m$, where $h$ is the Planck constant and $m$ is the mass
of an atom. When modelling with the GPE, the natural units to use are the
healing length units. Distances are measured in terms of the healing length,
$\xi=\hbar/\sqrt{m\mu}$, and times in terms of, $\tau=\hbar/\mu$, where $\mu$
is the chemical potential. In these units, the circulations are $\pm 2\pi$. We
refer to a vortex of circulation $+2\pi$ as a vortex and a vortex of
circulation $-2\pi$ as an antivortex.

###  2.2 A Dissipative PVM in a Doubly-Periodic Square Box

To model the dynamics of the vortices we use a doubly-periodic square box of
length $L=2048\xi$. To account for interactions between vortices across the
periodic boundaries we use the following equations [15]:

$\frac{d{\tilde{x}_{i}}}{dt}=\sum_{j=1,j\neq
i}^{N_{v}}\kappa_{j}\sum_{m=-\infty}^{\infty}\frac{-\sin(\tilde{y}_{ij})}{\cosh(\tilde{x}_{ij}-2\pi
m)-\cos(\tilde{y}_{ij})}\,.$ (5) $\frac{d{\tilde{y}_{i}}}{dt}=\sum_{j=1,j\neq
i}^{N_{v}}\kappa_{j}\sum_{m=-\infty}^{\infty}\frac{\sin(\tilde{x}_{ij})}{\cosh(\tilde{y}_{ij}-2\pi
m)-\cos(\tilde{x}_{ij})}\,.$ (6)

As above, the circulation of vortex $i$ is $\kappa_{i}$ but its coordinates,
$(x_{i},y_{i})$, are rescaled by $2\pi/L$ to give
$(\tilde{x}_{i},\tilde{y}_{i})$, so that in these units the box has a length
of $2\pi$. $\tilde{x}_{i}-\tilde{x}_{j}$ and $\tilde{y}_{i}-\tilde{y}_{j}$ are
written as $\tilde{x}_{ij}$ and $\tilde{y}_{ij}$. The time units are also
scaled to give $\kappa_{i}=\pm 1$ The first sum is over all vortices except
for vortex $i$ and the second is over the integer $m$. Since the term in this
second sum falls with increasing $m$, an excellent level of precision is
achieved by curtailing the sum to be from $m=-10$ to $m=+10$.

We use a Dormand-Prince 8th Order Runge-Kutta algorithm [16] to solve these
equations of motion, with an absolute tolerance of $1\times 10^{-8}$. This
gives us a non-dissipative velocity,
$\mathbf{v}_{i}=(d\tilde{x}_{i}/dt,d\tilde{y}_{i}/dt)$, for each vortex. To
include the effects of dissipation we add a dissipative component at right
angles to this to give the overall velocity, including dissipation,
$\mathbf{v}_{i}^{\rm{diss}}=\mathbf{v}_{i}-\gamma\kappa_{i}\hat{\mathbf{z}}\times\mathbf{v}_{i}$,
[11]. The values of $\gamma$ we use are
$\\{2^{-8},2^{-7},\,.\,.\,.\,,2^{-1},1,\infty\\}$, where $\gamma=\infty$ is
the purely dissipative case, in which
$\mathbf{v}_{i}^{\rm{diss}}=-\kappa_{i}\hat{\mathbf{z}}\times\mathbf{v}_{i}$.

A final, key aspect of vortex dynamics that we add to the simulation is the
annihilation of vortices and antivortices that come into too close proximity.
This is observed in experiments and in GPE simulations and is the main
mechanism by which phase ordering takes place. To include this in our
dissipative PVM, at each time-step we remove any vortex-antivortex pair that
are within $1\xi$ of each other.

###  2.3 Generating Initial Conditions

The initial arrangement of vortices and antivortices affects the dynamics that
follow. For example, an isolated vortex-antivortex pair, known as a dipole,
will move in a straight line, perpendicular to the line joining the vortices.
Alternatively, a cluster made up only of vortices (or only of antivortices)
will tend to rotate about its centre, so long as no other vortices are near.
The effect of dissipation on dipoles is to draw the constituent vortex and
antivortex closer together over time, whereas it tends to scatter clusters
apart. Therefore, when generating multiple realizations of the initial
conditions for repeat runs of the point vortex simulations it is vital to
ensure they have simlar distributions of dipoles, clusters and free vortices.
To achieve this, we note that the Hamiltonian, $H_{\square}$, for a system of
point vortices in a doubly-periodic square box [15], gives an unambiguous
parameter for measuring vortex configurations:

$H_{\square}=-\sum_{i=1}^{N_{v}}\sum_{j=1,j\neq
i}^{N_{v}}\frac{\kappa_{i}\kappa_{j}}{2}\sum_{m=-\infty}^{\infty}\ln{\left(\frac{\cosh{\tilde{x}_{ij}-2\pi
m}-\cos{\tilde{y}_{ij}}}{\cosh{2\pi
m}}\right)}-\frac{\tilde{x}_{ij}^{2}}{2\pi}\,.$ (7)

In this equation all symbols have the same meanings as in Eq. 5 & 6. The
double sum at the beginning adds contributions from all pairs of vortices and,
as before, an excellent level of precision is achieved by curtailing the sum
over $m$ to be from $m=-10$ to $m=+10$. Calculating this value for the initial
configuration and dividing the result by the number of vortices, $N_{v}$,
gives the initial energy per vortex $E^{0}/N_{v}$. A completely random
configuration of vortices and antivortices would have $E^{0}/N_{v}=0$,
increasing the amount of clustering of like-signed vortices increases
$E^{0}/N_{v}$ and pairing vortices and antivortices into dipoles makes
$E^{0}/N_{v}$ more negative.

To create initial conditions with specified values of $E^{0}/N_{v}$ we first
scatter $2069$ vortices and $2069$ antivortices in the box. Then we employ a
biased random walk algorithm in which the positions of vortices are updated at
random and these updated positions accepted only if they bring $E^{0}/N_{v}$
closer to the desired value. Once the desired values is achieved to within a
$1\%$ tolerance, further steps are undertaken in order to thermalize the
configuration. In each step, two vortices are randomly selected and moved by
different, random amounts. If $E^{0}/N_{v}$ remains within the tolerance the
step is accepted. Throughout all of this process, changes that bring vortices
closer than $1\xi$ are rejected.

Using this approach we create $100$ stochastic realizations of vortex
coordinates for each value of $E^{0}/N_{v}=\\{-4,-3,-2,-1,0.001,1,2,3,4\\}\pm
1\%$. These provide the initial conditions for our dissipative PVM.

##  3 A Method for Determining the Scaling Region

###  3.1 Constructing the Velocity Field and Correlation Function

In order to demonstrate dynamical scaling, a two-point correlation function is
generally calculated using a relevant scalar order parameter. In GPE
simulations this may be the quantum mechanical phase or the fluid density, in
the Ising model it is the magnetization, in our case we use the speed of the
fluid. To find this we first construct the velocity field from the coordinates
of the point vortices. We consider a $100\times 100$ grid of points
$(\tilde{x}_{g},\tilde{y}_{g})$, filling the square box of length $2\pi$ and
find the velocity of the fluid at each grid point,
$\mathbf{v}=(v_{\tilde{x}},v_{\tilde{y}})$. This is done by summing the
contributions, at that location, of all of the point vortices, whilst
accounting for the periodic boundary conditions [15]:

$v_{\tilde{x}}(\tilde{x}_{g},\tilde{y}_{g})=\sum_{i}\kappa_{i}\sum_{m=-\infty}^{\infty}\frac{-\sin(\tilde{y}_{g}-\tilde{y}_{i})}{\cosh(\tilde{x}_{g}-\tilde{x}_{i}-2\pi
m)-\cos(\tilde{y}_{g}-\tilde{y}_{i})}\,,$ (8)
$v_{\tilde{y}}(\tilde{x}_{g},\tilde{y}_{g})=\sum_{i}\kappa_{i}\sum_{m=-\infty}^{\infty}\frac{\sin(\tilde{x}_{g}-\tilde{x}_{i})}{\cosh(\tilde{y}_{g}-\tilde{y}_{i}-2\pi
m)-\cos(\tilde{x}_{g}-\tilde{x}_{i})}\,.$ (9)

Eq. 8 & 9 are almost the same as the non-dissipative point vortex equations of
motion, Eq. 5 & 6, except for two differences; the sum is over all of the
vortices and the coordinates, $(\tilde{x}_{g},\tilde{y}_{g})$, are those of a
grid point and not those of a vortex. This is because the point vortex
equations arise from considering the contributions of all other vortices to
the fluid velocity at that point and asserting that, in the absence of
dissipation, the vortex moves at the local fluid velocity. See Fig. 2 (a)-(c)
for examples of the velocity field constructed from the point vortex
coordinates in this fashion.

Using the velocity, $\mathbf{v}=(v_{\tilde{x}},v_{\tilde{y}})$, we find the
corrected speed at each grid point, $v_{c}$, from the difference between the
speed, $|\mathbf{v}|$, and the mean speed taken across the whole grid,
$\langle|\mathbf{v}|\rangle$:

$v_{c}=|\mathbf{v}|-\langle|\mathbf{v}|\rangle\,.$ (10)

This is the scalar order parameter we use to find the two-point correlation
function and ultimately provides the basis for demonstrating dynamical
scaling.

Figure 2: Fluid velocity calculated from point vortex configurations and two-
point correlation functions for the initial conditons $\gamma=1$,
$E^{0}/N_{v}=4$. (a)-(c) Open (white) circles show the locations of vortices
and filled (blue) circles show the locations of antivortices. The pale red
streamlines show the fluid velocity at the times
$t=\\{4200,13200,32700\\}\tau$. (d)-(f) Are the two-point correlation function
$G(r)$ calculated from the corrected speed at the same times.

The correlation function, $G$, for the corrected speed, $v_{c}$, between two
points separated by a displacement $(\tilde{x},\tilde{y})$, is given by:

$G(\tilde{x},\tilde{y})=\frac{\sum_{i}\sum_{j}v_{c}(\tilde{x}_{i},\tilde{y}_{j})v_{c}(\tilde{x}_{i}-\tilde{x},\tilde{y}_{j}-\tilde{y})}{\sum_{i}\sum_{j}v_{c}(\tilde{x}_{i},\tilde{y}_{j})^{2}}\,.$
(11)

Where the sum over $i$ is across all grid points in the $\tilde{x}$ direction,
the sum over $j$ is for all those in the $\tilde{y}$ direction and the
denominator ensures normalization such that $G(0,0)=1$. The Wiener-Khinchin
theorem [17] provides a more succint, and computationally more efficient,
formula in terms of Fourier transforms:

$G(\tilde{x},\tilde{y})=\mathcal{F}^{-1}[|\mathcal{F}(v_{c})|^{2}]\,.$

To convert from $G(\tilde{x},\tilde{y})$ to a one-dimensional $G(\tilde{r})$
we consider nested annuli of thickness, $\delta\tilde{r}$, where
$\delta\tilde{r}$ is the distance between adjacent grid points and
$\tilde{r}^{2}=\tilde{x}^{2}+\tilde{y}^{2}$. We then find $G(\tilde{r})$ by
calculating the mean value of $G(\tilde{x},\tilde{y})$ across all the grid
points in the annulus from $\tilde{r}$ to $\tilde{r}+\delta\tilde{r}$.
Finally, we multiply the distances by $L/2\pi$ to return to the original
length scale used for our system and recover $G(r)$. We repeat this for each
time-step, for all $100$ stochastic realizations and for the whole ensemble of
initial conditions.

###  3.2 Rescaling the Correlation Function and Testing for Dynamical Scaling

The theory of dynamical scaling predicts that the correlation function,
$G(r,t)$ should vary only by a change in length scale as time progresses. If
this is the case, finding a suitable length at each time-step and rescaling
the distances by this length, should cause the plots of $G(r)$ at different
times to collapse onto a single curve.

The process we use to rescale the correlation function for a given set of
initial conditions is as follows:

  1. 1.

Consider one time-step and find the mean values of $G(r)$ across the $100$
stochastic realizations (from here on $G$ should be understood as the average
over these realizations).

  2. 2.

Use cubic spline interpolation to find the radius at which $G(r)=0.1$. We
define this as the correlation length at this time-step, $L_{c}(t)$. (We use
$0.1$ for our data as it gives sufficient points, at both smaller and larger
radii, to ensure a reliable interpolation).

  3. 3.

Rescale the values of radius by dividing them by $L_{c}(t)$.

  4. 4.

Repeat the steps 1-3 for all time-steps.

This gives arrays of data for $G$ and $r/L_{c}$ at each time-step. Simply
plotting $G$ against $r/L_{c}$ allows us to check visually if the curves for
different time-steps collapse. However, a more systematic approach is desired
to determine the start and end of any period of dynamical scaling in a robust
and reliable fashion. The method we implement to test for collapse of the
different curves is:

  1. 1.

For every time-step, use spline interpolation to find the values of $G$ on the
same uniform grid of values of $r/L_{c}$. Use $N_{g}$ points in the range from
the smallest value of $r/L_{c}(t)$ for any time-step up to $r/L_{c}(t)=5$. By
this point the values of $G$ should be very close to zero in all cases, so
continuing further is unnecessary. We use $N_{g}=40$.

  2. 2.

For each time-step there are now $N_{g}$ values of $G$. Therefore each time-
step can now be represented by a point in an $N_{g}$-dimensional Euclidean
space:

$\mathbf{G}^{i}=(G_{1}^{i},G_{2}^{i},...,G_{N_{g}}^{i})\,,$ (12)

where, $\mathbf{G}^{i}$ is the vector position representing the correlation
function at timestep $i$ and $(G_{1}^{i},G_{2}^{i},...,G_{N_{g}}^{i})$ are the
coordinates in each of the dimensions (i.e. the values of $G$ on the grid of
$N_{g}$ radii). The problem is now one of determining which time-steps are
outliers and which are not, see Fig. 3.

Figure 3: Representing curves as points in an $N_{g}$-dimensional Euclidean
space. (a) Five curves $G(r)$ plotted on a grid of three points
${r_{1},r_{2},r_{3}}$. (b) The same curves represented by points in a three-
dimensional Euclidean space. Similar curves, such as those shown with red
circles, green squares, blue diamonds and magenta triangles, become points
that are clustered close together; whilst a very different curve, like the one
shown with orange stars, becomea a point very far from the others. To identify
an outlier, we need to quantify the deviation, $\Delta$, of a suspected point
from the mean position of all the other points. We use the Mahalanobis
distance, in the $N_{g}$-dimensional Euclidean space, to do this.

  3. 3.

Find the mean, $\boldsymbol{\mu}$, and the covariance, $\mathbf{C}$, of
$\mathbf{G}^{i}$ across all of the time-steps. $\boldsymbol{\mu}$ is a vector
of length $N_{g}$ and $\mathbf{C}$ is a matrix of size $N_{g}\times N_{g}$.

  4. 4.

The squared deviation of a point, $\mathbf{G}^{i}$, from the mean,
$\boldsymbol{\mu}$, normalized by the covariance between the data from the
orginal curves, is given by the Mahalanobis distance squared, $MD_{i}^{2}$,
[4]:

$MD_{i}^{2}=(\mathbf{G}^{i}-\boldsymbol{\mu})\mathbf{C}^{-1}(\mathbf{G}^{i}-\boldsymbol{\mu})^{T}\,.$
(13)

Calculate $MD_{i}^{2}$ for each time-step and order them from largest to
smallest. Then, to identify and eliminate outliers:

  5. 5.

Choose the time-step with the largest value of $MD_{i}^{2}$ and recalculate
$\boldsymbol{\mu}$ and $\mathbf{C}$ without it (i.e. treate it as an outlier).

  6. 6.

Recalculate $MD_{i}^{2}$ for the suspected outlier using the new
$\boldsymbol{\mu}$ and $\mathbf{C}$. If the value is below a cutoff (see
discussion below), this timestep is not an outlier and the process is stopped.
If it is above the cutoff, this timestep is identified as an outlier and
removed from future consideration.

  7. 7.

Recalculate $MD_{i}^{2}$ for all the other points, with the previously
identified outlier(s) excluded and order these values from largest to
smallest.

  8. 8.

Identify the time-step with the new largest value of $MD_{i}^{2}$ and repeat
steps 5-7.

  9. 9.

Continue iteratively until a suspected outlier has $MD_{i}^{2}$ less than the
cutoff, at which point all outliers have been identified.

  10. 10.

All remaining timesteps are considered to overlap. The scaling region is then
the longest stretch of uninterrupted timesteps not containing any outliers.
See Fig. 4 for an example.

The choice of cutoff for determining outliers is a slighlty subtle one. If the
values of $\mathbf{G}^{i}$ at each radius are independent and normally
distributed then $MD_{i}^{2}$ would qualify as a $\chi^{2}$ statistic and the
critical $\chi^{2}$ value for a suitable significance level would be an
appropriate choice. However, our data are not normally distributed, so instead
we use a cutoff of $4\times N_{g}$, where $N_{g}$ is the number of grid-
points. The logic behind this decision is that for the case in which
$\mathbf{G}^{i}$ at each radius varies independently of every other radius,
the covariance matrix is diagonal and $MD_{i}^{2}$ becomes a sum of similar
terms:

$MD_{i}^{2}=(\mathbf{G}^{i}-\boldsymbol{\mu})\begin{pmatrix}1/\sigma_{1}^{2}&0&.&.&0\\\
0&1/\sigma_{2}^{2}&.&.&0\\\ .&.&.&.&.\\\ .&.&.&.&.\\\
0&0&.&.&1/\sigma_{N_{g}}^{2}\end{pmatrix}(\mathbf{G}^{i}-\boldsymbol{\mu})^{T}\,=\sum_{k=1}^{N_{g}}\left(\frac{G^{i}_{k}-\mu_{k}}{\sigma_{k}}\right)^{2}\,.$
(14)

So, if a hypothetical curve has every point one standard deviation from the
mean, then $MD^{2}=\sum_{k=1}^{N_{g}}1^{2}=N_{g}$, whereas if each point is
two standard deviations from the mean then
$MD^{2}=\sum_{k=1}^{N_{g}}2^{2}=4N_{g}$. We consider this a reasonable choice
of cutoff, whilst recognizing that our values of $\mathbf{G}^{i}$ at each
radius are not, in general, independent. However, a lower value could be used
if a stricter definition of collapse is preferred. The cutoff is the one free
parameter in this approach and plays a similar role to the significance level
in hypothesis testing. Fig. 5 (a) shows how the choice of cutoff affects the
percentage of timesteps identified as within the scaling region and
consequently the values of dynamical critical exponent determined from the
scaling region data.

Figure 4: Identifying the scaling region and determing $z$ for the initial
conditions $\gamma=1,\,E^{0}/N_{v}=4$. (a) Shows the two-point correlation
function $G(r)$ for each time-step. (b) Shows $G(r)$ with the distances
rescaled by $L_{c}$, where $L_{c}$ is the value of $r$ when the unscaled
$G(r)=0.1$. In both (a) and (b) those time-steps identified as within the
scaling region are plotted as solid red lines, getting darker as time
progresses, and those outside of the scaling region are plotted as gray dashed
lines. In the rescaled plots, (b), all curves in the scaling region overlap.
(c) Shows the data for correlation length against time plotted as grey
crosses, the bars show the standard error in $L_{c}$ calculated across the
$100$ stochastic realizations at each time-step. The red line is a power law
$L_{c}\sim t^{1/z}$, with $z=1.985\pm 0.003$ determined using a bootstrap
method [18] in which we resample $L_{c}(t)$ $100$ times, conduct a power law
fit each time and find the mean and standard error in $z$ across the
resamplings. (d) Illustrates the time-steps identified as outliers as black
blocks in the bottom line and those in the scaling region as red blocks on the
top line.

###  3.3 The Dynamical Critical Exponent

To calculate values for the dynamical critical exponent, $z$, we conduct a
least squares fit of the power law $L_{c}\sim t^{1/z}$ to the data $L_{c}(t)$
for just the times within the scaling region. By identifying the scaling
region in a robust fashion, with a suitable choice of cutoff, prior to
conducting the fit, we can be more confident that the predicted power law
relationship is valid and the value of $z$ found is an accurate measure of the
dynamical critical exponent of the system. Fig. 5 (b) shows how the value of
$z$, determined in this way, varies with the level of dissipation, $\gamma$.

Figure 5: Determining the dynamical critical exponent, $z$, for
$E^{0}/N_{v}=4$. (a) The black crosses (left-hand axis) show how the value of
$z$ varies with the choice of cutoff for $\gamma=1$ (errorbars are not shown
as they are either too large, for cutoff $<3N_{g}$, or too small, for cutoff
$\geq 3N_{g}$). The red diamonds (right-hand axis) show how the % of timesteps
identified as within the scaling region varies. Using a cutoff of less than
$3N_{g}$ is problematic as only a tiny fraction of timesteps are within the
scaling region and the fit is unreliable. (b) The black crosses (left-hand
axis) show how the value of $z$ varies with the dissipation, $\gamma$, for a
cutoff of $4N_{g}$ (errorbars would be too small to see on this scale). The
red diamonds (right-hand axis) show the corresponding % of timesteps in the
scaling region. A large scaling region is identified for all $\gamma$ and
$z\approx 2$ but dependent upon $\gamma$, a qualitatively similar result to
that reported in [3]. In all cases we determine $z$ using the same bootstrap
method described in Fig. 4.

##  4 Conclusions

We have described a systematic method, based on the concept of Mahalanobis
distance, for identifying outliers when comparing graphs of rescaled
correlation functions. By eliminating those times for which the curves do not
collapse onto one single curve, our approach gives a reliable identification
of the scaling region. We have also demonstrated its application to the phase-
ordering behaviour of a quantum fluid following a rapid quench, using data
generated using a dissipative PVM. This produced results that correspond
qualitatively and quantitatively with those from other studies that make use
of GPE simulations and visual identification of the scaling region (with
subsequent checks) [3, 5]. In particular, we found values of dynamical
critical constant that are close to two but depend on dissipation. By applying
this method we can be confident that the data used to determine the dynamical
critical exponent came only from times during which the system displays
scaling behaviour. Further work could include an application of this approach
to the correlation functions for a wider ensemble of initial conditions, to
the results of GPE simulations or to studies of non-equilibrium behaviours of
systems outside of the field of quantum fluids.

## Declarations

  * •

Funding This work was supported by a Lady Bertha Jeffreys PhD Studentship.

  * •

Competing interests The authors have no competing interests to declare that
are relevant to the content of this article.

  * •

For the purpose of open access, the author has applied a Creative Commons
Attribution (CC BY) licence to any Author Accepted Manuscript version arising
from this submission.

## References

  *   * Bray [1994] Bray, A.J.: Theory of phase-ordering kinetics. Advances in Physics 43, 357–459 (1994) https://doi.org/10.1080/00018739400101505
  * Nishimori [2011] Nishimori, H.: Elements of Phase Transitions and Critical Phenomena. Oxford University Press, Oxford (2011)
  * Groszek et al. [2021] Groszek, A.J., Comaron, P., Proukakis, N.P., Billam, T.P.: Crossover in the dynamical critical exponent of a quenched two-dimensional Bose gas. Physical Review Research 3, 013212 (2021) https://doi.org/10.1103/PhysRevResearch.3.013212
  * Mahalanobis [2018] Mahalanobis, P.C.: On the generalised distance in statistics. Sankhya A 80, 1–7 (2018) https://doi.org/10.1007/s13171-019-00164-5
  * Karl and Gasenzer [2017] Karl, M., Gasenzer, T.: Strongly anomalous non-thermal fixed point in a quenched two-dimensional Bose gas. New Journal of Physics 19(9), 093014 (2017) https://doi.org/10.1088/1367-2630/aa7eeb
  * Pitaevskii and Stringari [2016] Pitaevskii, L., Stringari, S.: Bose-Einstein Condensation and Superfluidity. Oxford University Press, Oxford (2016). https://doi.org/%****␣RT_QFS.bbl␣Line␣125␣****10.1093/acprof:oso/9780198758884.001.0001
  * Pethick and Smith [2008] Pethick, C.J., Smith, H.: Bose–Einstein Condensation in Dilute Gases. Cambridge University Press, Cambridge (2008). https://doi.org/10.1017/CBO9780511802850
  * Helmholtz [1858] Helmholtz, H.: Über intergrale der hydrodynamischen gleichungen, welche den wirbelbewegungen entsprechen. J. Reine Agnew. Math. 55, 22–55 (1858)
  * Kirchoff [1876] Kirchoff, G.: Vorlesungen Über Matematische Physik. Mechanik. Teuber, Leipzig (1876)
  * Fetter [1966] Fetter, A.L.: Vortices in an imperfect bose gas. IV. translational velocity. Phys. Rev. 151, 100–104 (1966) https://doi.org/%****␣RT_QFS.bbl␣Line␣175␣****10.1103/PhysRev.151.100
  * Billam et al. [2015] Billam, T.P., Reeves, M.T., Bradley, A.S.: Spectral energy transport in two-dimensional quantum vortex dynamics. Phys. Rev. A 91, 023615 (2015) https://doi.org/10.1103/PhysRevA.91.023615
  * Skipp et al. [2023] Skipp, J., Laurie, J., Nazarenko, S.: Hamiltonian derivation of the point vortex model from the two-dimensional nonlinear Schrödinger equation. Phys. Rev. E 107, 025107 (2023) https://doi.org/10.1103/PhysRevE.107.025107
  * Anderson et al. [1995] Anderson, M.H., Ensher, J.R., Matthews, M.R., Wieman, C.E., Cornell, E.A.: Observation of Bose-Einstein condensation in a dilute atomic vapor. Science 269, 198–201 (1995) https://doi.org/10.1126/science.269.5221.198
  * Abo-Shaeer et al. [2001] Abo-Shaeer, J.R., Raman, C., Vogels, J.M., Ketterle, W.: Observation of vortex lattices in Bose-Einstein condensates. Science 292, 476–479 (2001) https://doi.org/10.1126/science.1060182
  * Weiss and McWilliams [1991] Weiss, J.B., McWilliams, J.C.: Nonergodicity of point vortices. Physics of Fluids A: Fluid Dynamics 3, 835–844 (1991) https://doi.org/10.1063/1.858014
  * Prince and Dormand [1981] Prince, P.J., Dormand, J.R.: High order embedded Runge-Kutta formulae. Journal of Computational and Applied Mathematics 7(1), 67–75 (1981) https://doi.org/10.1016/0771-050X(81)90010-3
  * Champeney [1987] Champeney, D.: Power spectra and Wiener’s theorems. In: A Handbook of Fourier Theorems, pp. 102–117. Cambridge University Press, Cambridge (1987)
  * Efron [1993] Efron, B.: An Introduction to the Bootstrap. Monographs on statistics and applied probability (Series) ; 57. Chapman & Hall, New York (1993)