# Generating Steganographic Text with LSTMs

Tina Fang
University of Waterloo
<EMAIL_ADDRESS>
Martin Jaggi
EPFL
<EMAIL_ADDRESS>
Katerina Argyraki
EPFL
<EMAIL_ADDRESS>
###### Abstract

Motivated by concerns for user privacy, we design a steganographic system
(“stegosystem”) that enables two users to exchange encrypted messages without
an adversary detecting that such an exchange is taking place. We propose a new
linguistic stegosystem based on a Long Short-Term Memory (LSTM) neural
network. We demonstrate our approach on the Twitter and Enron email datasets
and show that it yields high-quality steganographic text while significantly
improving capacity (encrypted bits per word) relative to the state-of-the-art.

##  1 Introduction

The business model behind modern communication systems (email services or
messaging services provided by social networks) is incompatible with end-to-
end message encryption. The providers of these services can afford to offer
them free of charge because most of their users agree to receive “targeted
ads” (ads that are especially chosen to appeal to each user, based on the
needs the user has implied through their messages). This model works as long
as users communicate mostly in the clear, which enables service providers to
make informed guesses about user needs.

This situation does not prevent users from encrypting a few sensitive
messages, but it does take away some of the benefits of confidentiality. For
instance, imagine a scenario where two users want to exchange forbidden ideas
or organize forbidden events under an authoritarian regime; in a world where
most communication happens in the clear, encrypting a small fraction of
messages automatically makes these messages—and the users who exchange
them—suspicious.

With this motivation in mind, we want to design a system that enables two
users to exchange encrypted messages, such that a passive adversary that reads
the messages can determine neither the original content of the messages nor
the fact that the messages are encrypted.

We build on linguistic steganography, i.e., the science of encoding a secret
piece of information (“payload”) into a piece of text that looks like natural
language (“stegotext”). We propose a novel stegosystem, based on a neural
network, and demonstrate that it combines high quality of output (i.e., the
stegotext indeed looks like natural language) with the highest capacity
(number of bits encrypted per word) published in literature.

In the rest of the paper, we describe existing linguistic stegosystems along
with ours (§2), provide details on our system (§3), present preliminary
experimental results on Twitter and email messages (§4), and conclude with
future directions (§5).

##  2 Linguistic Steganography

In this section, we summarize related work (§2.1), then present out proposal
(§2.2).

###  2.1 Related Work

Traditional linguistic stegosystems are based on modification of an existing
cover text, e.g., using synonym substitution Topkara et al. (2006); Chang and
Clark (2014) and/or paraphrase substitution Chang and Clark (2010). The idea
is to encode the secret information in the transformation of the cover text,
ideally without affecting its meaning or grammatical correctness. Of these
systems, the most closely related to ours is CoverTweet Wilson et al. (2014),
a state-of-the-art cover modification stegosystem that uses Twitter as the
medium of cover; we compare to it in our preliminary evaluation (§4).

Cover modification can introduce syntactic and semantic unnaturalness Grosvald
and Orgun (2011); to address this, Grovsald and Orgun proposed an alternative
stegosystem where a human generates the stegotext manually, thus improving
linguistic naturalness at the cost of human effort Grosvald and Orgun (2011).

Matryoshka Safaka et al. (2016) takes this further: in step 1, it generates
candidate stegotext automatically based on an $n$-gram model of the English
language; in step 2, it presents the candidate stegotext to the human user for
polishing, i.e., ideally small edits that improve linguistic naturalness.
However, the cost of human effort is still high, because the (automatically
generated) candidate stegotext is far from natural language, and, as a result,
the human user has to spend significant time and effort manually editing and
augmenting it.

Volkhonskiy et al. have applied Generative Adversarial Networks Goodfellow et
al. (2014) to image steganography Volkhonskiy et al. (2017), but we are not
aware of any text stegosystem based on neural networks.

###  2.2 Our Proposal: Steganographic LSTM

Motivated by the fact that LSTMs Hochreiter and Schmidhuber (1997) constitute
the state of the art in text generation Jozefowicz et al. (2016), we propose
to automatically generate the stegotext from an LSTM (as opposed to an
$n$-gram model). The output of the LSTM can then be used either directly as
the stegotext, or Matryoshka-style, i.e., as a candidate stegotext to be
polished by a human user; in this paper, we explore only the former option,
i.e., we do not do any manual polishing. We describe the main components of
our system in the paragraphs below; for reference, Fig. 1 outlines the
building blocks of a stegosystem Salomon (2003).

Figure 1:  Stegosystem building blocks.

##### Secret data.

The secret data is the information we want to hide. First, we compress and/or
encrypt the secret data (e.g., in the simplest setting using the ASCII coding
map) into a _secret-containing bit string_ $S$. Second, we divide $S$ into
smaller _bit blocks_ of length $|B|$, resulting in a total of $|S|/|B|$111If
$|B|\centernot\mid|S|$, then we leave the remainder bit string out of
encryption. bit blocks. For example, if $S=100001$ and $|B|=2$, our bit-block
sequence is $10$, $00$, $01$. Based on this bit-block sequence, our
steganographic LSTM generates words.

##### Key.

The sender and receiver share a key that maps bit blocks to token sets and is
constructed as follows: We start from the _vocabulary_ , which is the set of
all possible tokens that may appear in the stegotext; the tokens are typically
words, but may also be punctuation marks. We partition the vocabulary into
$2^{|B|}$ _bins_ , i.e., disjoint token sets, randomly selected from the
vocabulary without replacement; each token appears in exactly one bin, and
each bin contains $|V|/2^{|B|}$ tokens. We map each bit block $B$ to a bin,
denoted by $W_{B}$. This mapping constitutes the shared key.

Bit Block | Tokens
---|---
00 | This, am, weather, …
01 | was, attaching, today, …
10 | I, better, an, Great, …
11 | great, than, NDA, ., …
Table 1:  Example shared key.

##### Embedding algorithm.

The embedding algorithm uses a modified word-level LSTM for language modeling
Mikolov et al. (2010). To encode the secret-containing bit string $S$, we
consider one bit block $B$ at a time and have our LSTM select one token from
bin $W_{B}$; hence, the candidate stegotext has as many tokens as the number
of bit blocks in $S$. Even though we restrict the LSTM to select a token from
a particular bin, each bin should offer sufficient variety of tokens, allowing
the LSTM to generate text that looks natural. For example, given the bit
string “1000011011” and the key in Table 1, the LSTM can form the partial
sentence in Table 2. We describe our LSTM model in more detail in the next
section.

Bit String | 10 | 00 | 01 | 10 | 11
---|---|---|---|---|---
Token | I | am | attaching | an | NDA
Table 2:  Example stegotext generation.

##### Decoder.

The decoder recovers the original data deterministically and in a
straightforward manner: it takes as input the generated stegotext, considers
one token at a time, finds the token’s bin in the shared key, and recovers the
original bit block.

##### Common-token variant.

We also explore a variant where we add a set of _common tokens_ , $C$, to all
bins. These common tokens do not carry any secret information; they serve only
to enhance stegotext naturalness. When the LSTM selects a common token from a
bin, we have it select an extra token from the same bin, until it selects a
non-common token. The decoder removes all common tokens before decoding. We
discuss the choice of common tokens and its implication on our system’s
performance in Section 4.

##  3 Steganographic LSTM Model

In this section, we provide more details on our system: how we modify the LSTM
(§3.1) and how we evaluate its output (§3.2).

###  3.1 LSTM Modification

##### Text generation in classic LSTM.

Classic LSTMs generate words as follows Sutskever et al. (2011): Given a word
sequence $(x_{1},x_{2},\dots,x_{T})$, the model has hidden states
$(h_{1},\dots,h_{T})$, and resulting output vectors $(o_{1},\dots,o_{T})$.
Each output vector $o_{t}$ has length $|V|$, and each output-vector element
$o_{t}^{(j)}$ is the unnormalized probability of word $j$ in the vocabulary.
Normalized probabilities for each candidate word are obtained by the following
softmax activation function:

$\mathop{softmax}(o_{t})_{j}:=\exp(o_{t}^{(j)})\Big{/}\sum_{k}\exp(o_{t}^{(k)}).$

The LSTM then selects the word with the highest probability
$P[x_{t+1}\,|\,x_{\leq t}]$ as its next word.

##### Text generation in our LSTM.

In our steganographic LSTM, word selection is restricted by the shared key.
That is, given bit block $B$, the LSTM has to select its next word from bin
$W_{B}$. We set $P[x=w_{j}]=0$ for $j\notin W_{B}$, so that the multinomial
softmax function selects the word with the highest probability within $W_{B}$.

##### Common tokens.

In the common-token variant, we restrict $P[x=w_{j}]=0$ only for
$j\notin(W_{B}\cup C)$, where $C$ is the set of common tokens added to all
bins.

###  3.2 Evaluation Metrics

We use perplexity to quantify stegotext quality; and capacity (i.e., encrypted
bits per output word) to quantify its efficiency in carrying secret
information. In Section 4, we also discuss our stegotext quality as
empirically perceived by us as human readers.

##### Perplexity.

Perplexity is a standard metric for the quality of language models Martin and
Jurafsky (2000), and it is defined as the average per-word log-probability on
the valid data set: $\exp(-1/N\sum_{i}\ln p[w_{i}])$ Jozefowicz et al. (2016).
Lower perplexity indicates a better model.

In our steganographic LSTM, we cannot use this metric as is: since we enforce
$p[w_{i}]=0$ for $w_{i}\notin W_{B}$, the corresponding $\ln p[w_{i}]$ becomes
undefined under this vocabulary.

Instead, we measure the probability of $w_{i}$ by taking the average of
$p[w_{i}]$ over all possible secret bit blocks $B$, under the assumption that
bit blocks are distributed uniformly. By the Law of Large Numbers Révész
(2014), if we perform many stegotext-generating trials using different random
secret data as input, the probability of each word will tend to the expected
value, $\sum p[w_{i},B]/2^{|B|}$, Hence, we set $p[w_{i}]:=\sum
p[w_{i},B]/2^{|B|}$ instead of $p[w_{i}]=0$ for $w_{i}\notin W_{B}$.

##### Capacity.

Our system’s capacity is the number of encrypted bits per output word. Without
common tokens, capacity is always $|B|$ bits/word (since each bit block of
size $|B|$ is always mapped to one output word). In the common-token variant,
capacity decreases because the output includes common tokens that do not carry
any secret information; in particular, if the fraction of common tokens is
$p$, then capacity is $(1-p)\cdot|B|$.

##  4 Experiments

In this section, we present our preliminary experimental evaluation: our
Twitter and email datasets (§4.1), details about the LSTMs used to produce our
results (§4.2), and finally a discussion of our results (§4.3).

###  4.1 Datasets

Tweets and emails are among the most popular media of open communication and
therefore provide very realistic environments for hiding information. We thus
trained our LSTMs on those two domains, Twitter messages and Enron emails
Klimt and Yang (2004), which vary greatly in message length and vocabulary
size.

For Twitter, we used the NLTK tokenizer to tokenize tweets Bird (2006) into
words and punctuation marks. We normalized the content by replacing usernames
and URLs with a username token ($<$user$>$) and a URL token ($<$url$>$),
respectively. We used $600$ thousand tweets with a total of $45$ million words
and a vocabulary of size $225$ thousand.

For Enron, we cleaned and extracted email message bodies Zhou et al. (2007)
from the Enron dataset, and we tokenized the messages into words and
punctuation marks. We took the first $100$MB of the resulting messages, with
$16.8$ million tokens and a vocabulary size of $406$ thousand.

###  4.2 Implementation Details

We implemented multi-layered LSTMs based on
PyTorch222https://github.com/pytorch in both experiments. We did not use pre-
trained word embeddings Mikolov et al. (2013); Pennington et al. (2014), and
instead trained word embeddings of dimension $200$ from scratch.

We optimized with Stochastic Gradient Descent and used a batch size of $20$.
The initial learning rate was $20$ and the decay factor per epoch was $4$. The
learning rate decay occurred only when the validation loss did not improve.
Model training was done on an NVIDIA GeForce GTX TITAN X.

For Twitter, we used a $2$-layer LSTM with $600$ units, unrolled for $25$
steps for back propagation. We clipped the norm of the gradients Pascanu et
al. (2013) at $0.25$ and applied $20$% dropout Srivastava et al. (2014). We
stopped the training after $12$ epochs ($10$ hours) based on validation loss
convergence.

For Enron, we used a $3$-layer LSTM with $600$ units and no regularization. We
unrolled the network for $20$ steps for back propagation. We stopped the
training after $6$ epochs ($2$ days).

###  4.3 Results and Discussion

####  4.3.1 Tweets

We evaluate resulting tweets generated by LSTMs of $1$ (non-steganographic),
$2$, $4$, $8$ bins. Furthermore, we found empirically that adding $10$ most
frequent tokens from the Twitter corpus was enough to significantly improve
the grammatical correctness and semantic reasonableness of tweets. Table 3
shows the relationship between capacity (_bits per word_), and quantitative
text quality (_perplexity_). It also compares models with and without adding
common tokens using perplexity and bits per word.

Table 4 shows example output texts of LSTMs with and without common tokens
added. To reflect the variation in the quality of the tweets, we represent
tweets that are good and poor in quality333For each category, we manually
evaluate $60$ randomly generated tweets based grammatical correctness,
semantic coherence, and resemblance to real tweets. We select tweets from the
$25$th, $50$th, and $75$th percentile, and call them “good”, “average”, and
“poor” respectively. We limit to tweets that are not offensive in language..

We replaced $<$user$>$ generated by the LSTM with mock usernames for a more
realistic presentation in Table 4. In practice, we can replace the $<$user$>$
tokens systematically, randomly selecting followers or followees of that tweet
sender, for example.

Re-tweet messages starting with “RT” can also be problematic, because it will
be easy to check whether the original message of the retweeted message exists.
A simple approach to deal with this is to eliminate “RT” messages from
training (or at generation). Finally, since we made all tweets lower case in
the pre-processing step, we can also post-process tweets to adhere to proper
English capitalization rules.

| Original | Common Tokens
---|---|---
# of Bins | bpw | ppl | bpw | ppl
$1$ | $0$ | $134.73$ | $0$ | $134.73$
$2$ | $1$ | $190.84$ | $0.65$ | $171.35$
$4$ | $2$ | $381.2$ | $1.17$ | $277.55$
$8$ | $3$ | $833.11$ | $1.53$ | $476.66$
Table 3:  An increase of of capacity correlates with an increase of perplexity, which implies that there is a negative correlation between capacity and text quality. After adding common tokens, there is a significant reduction in perplexity _(ppl)_ , at the expense of a lower capacity _(bits per word)_. # of Bins |  Tweets |  Tweets with Common Tokens
---|---|---
$2$ |  good: i was just looking for someone that i used have. poor: cry and speak! rt @user421: relatable personal hygiene for body and making bad things as a best friend in lifee |  good: i’m happy with you. i’ll take a pic poor: rt: cut your hair, the smallest things get to the body.
$4$ |  good: @user390 looool yeah she likes me then ;). you did? poor: “where else were u making?… i feel fine? - e? lol” * does a voice for me & take it to walmart? |  good: i just wanna move. collapses. poor: i hate being contemplating for something i want to.
$8$ |  good: @user239 hahah. sorry that my bf is amazing because i’m a bad influence ;). poor: so happy this to have been working my ass and they already took the perfect. but it’s just cause you’re too busy the slows out! love… * dancing on her face, holding two count out cold * ( a link with a roof on punishment… - please :) |  good: i hate the smell of my house. poor: a few simple i can’t. i need to make my specs jump surprisingly.
Table 4:  We observe that the model with common tokens produces tweets simpler
in style, and uses more words from the set of common tokens. There is a large
improvement in grammatical correctness and context coherence after adding
common tokens, especially in the “poor” examples. For example, adding the line
break token reduced the length of the tweet generated from the 8-bin LSTM.

####  4.3.2 Emails

We also tested email generation, and Table 5 shows sample email passages444We
only present passages ”average” in quality to conserve space. from each bin.
We post-processed the emails with untokenization of punctuations.

The biggest difference between emails and tweets is that emails have a much
longer range for context dependency, with context spanning sentences and
paragraphs. This is challenging to model even for the non-steganographic LSTM.
Once the long-range context dependency of the non-steganographic LSTM
improves, the context dependency of the steganographic LSTMs should also
improve.

# of Bins |  Sample Email
---|---
$1$ |  —–Original Message—– From: Nelson, Michelle Sent: Thursday, January 03, 2002 3:35 PM To: Maggi, Mike Subject: Hey, You are probably a list of people that are around asleep about the point of them and our wife. Rob
$2$ |  If you do like to comment on the above you will not contact me at the above time by 8:00 a.m. on Monday, March 13 and July 16 and Tuesday, May 13 and Tuesday, March 9 - Thursday, June 17, - 9:00 to 11:30 AM.
$4$ |  At a moment when my group was working for a few weeks, we were able to get more flexibility through in order that we would not be willing.
Table 5:  The issue of context inconsistency is present for all bins. However,
the resulting text remains syntactical even as the number of bins increases.

###  4.4 Comparison with Other Stegosystems

For all comparisons, we use our 4-bin model with no common tokens added.

Our model significantly improves the state-of-the-art capacity. Cover
modification based stegosystems hide $1$-$2$ bits per sentence Chang and Clark
(2012). The state-of-the-art Twitter stegosystem hides $2.8$ bits of per tweet
Wilson and Ker (2016). Assuming $16.04$ words per tweet555Based a random
sample of $2$ million tweets., our $4$-bin system hides 32 bits per tweet,
over $11$ times higher than Wilson and Ker (2016).

We hypothesize that the subjective quality of our generated tweets will be
comparable to tweets produced by CoverTweet Wilson et al. (2014). We present
some examples666Tweets selected for comparison are “average” in quality. in
Table 6 to show there is potential for a comparison. This contrasts the
previous conception that cover generation methods are fatally weak against
human judges Wilson et al. (2014). CoverTweet was tested to be secure against
human judges. Formal experiments will be necessary to establish that our
system is also secure against human judges.

CoverTweet Wilson et al. (2014) |  Steganographic LSTM
---|---
yall must have 11:11 set 1 minute early before yall tweet it, because soon as 11:11 hit yall don’t wastes no time. lol |  i wanna go to sleep in the gym, ny in peoples houses & i’m in the gym..! :((
you can tell when somebody hating on you! |  i would rather marry a regular sunday!!
most of the people who got mouth can’t beat you. |  my mom is going so hard to get his jam.
Table 6:  The tweets generated by the 4-bin LSTM ($32$ bits per tweet) are
reasonably comparable in quality to tweets produced by CoverTweet ($2.8$ bits
per tweet).

Our system also offers flexibility for the user to freely trade-off capacity
and text quality. Though we chose the 4-bin model with no common tokens for
comparison, user can choose to use more bins to achieve an even higher
capacity, or use less bins and add common tokens to increase text quality.
This is not the case with existing cover modification systems, where capacity
is bounded above by the number of transformation options Wilson et al. (2014).

##  5 Conclusion and Future Work

In this paper, we opened a new application of LSTMs, namely, steganographic
text generation. We presented our steganographic model based on existing
language modeling LSTMs, and demonstrated that our model produces realistic
tweets and emails while hiding information.

In comparison to the state-of-the-art steganographic systems, our system has
the advantage of encoding much more information (around $2$ bits per word).
This advantage makes the system more usable and scalable in practice.

In future work, we will formally evaluate our system’s security against human
judges and other steganography detection (steganalysis) methods Wilson et al.
(2015); Kodovsky et al. (2012). When evaluated against an automated
classifier, the setup becomes that of a Generative Adversarial Network
Goodfellow et al. (2014), though with additional conditions for the generator
(the secret bits) which are unknown to the discriminator, and not necessarily
employing joint training. Another line of future research is to generate
tweets which are personalized to a user type or interest group, instead of
reflecting all twitter users. Furthermore, we plan to explore if capacity can
be improved even more by using probabilistic encoders/decoders, as e.g. in
Matryoshka (Safaka et al., 2016, Section 4).

Ultimately, we aim to open-source our stegosystem so that users of open
communication systems (e.g. Twitter, emails) can use our stegosystem to
communicate private and sensitive information.

## References

  * Bird (2006) Steven Bird. 2006.  Nltk: the natural language toolkit.  In Proceedings of the COLING/ACL on Interactive presentation sessions. Association for Computational Linguistics, pages 69–72.
  * Chang and Clark (2010) Ching-Yun Chang and Stephen Clark. 2010.  Linguistic steganography using automatically generated paraphrases.  In Human Language Technologies: The 2010 Annual Conference of the North American Chapter of the Association for Computational Linguistics. Association for Computational Linguistics, pages 591–599.
  * Chang and Clark (2012) Ching-Yun Chang and Stephen Clark. 2012.  Adjective deletion for linguistic steganography and secret sharing.  In COLING. pages 493–510.
  * Chang and Clark (2014) Ching-Yun Chang and Stephen Clark. 2014.  Practical linguistic steganography using contextual synonym substitution and a novel vertex coding method.  Computational Linguistics 40(2):403–448.
  * Goodfellow et al. (2014) Ian Goodfellow, Jean Pouget-Abadie, Mehdi Mirza, Bing Xu, David Warde-Farley, Sherjil Ozair, Aaron Courville, and Yoshua Bengio. 2014.  Generative adversarial nets.  In Advances in neural information processing systems. pages 2672–2680.
  * Grosvald and Orgun (2011) Michael Grosvald and C Orhan Orgun. 2011.  Free from the cover text: a human-generated natural language approach to text-based steganography.  Journal of Information Hiding and Multimedia Signal Processing 2(2):133–141.
  * Hochreiter and Schmidhuber (1997) Sepp Hochreiter and Jürgen Schmidhuber. 1997.  Long short-term memory.  Neural computation 9(8):1735–1780.
  * Jozefowicz et al. (2016) Rafal Jozefowicz, Oriol Vinyals, Mike Schuster, Noam Shazeer, and Yonghui Wu. 2016\.  Exploring the limits of language modeling.  arXiv preprint arXiv:1602.02410 .
  * Klimt and Yang (2004) Bryan Klimt and Yiming Yang. 2004.  The enron corpus: A new dataset for email classification research.  In European Conference on Machine Learning. Springer, pages 217–226.
  * Kodovsky et al. (2012) Jan Kodovsky, Jessica Fridrich, and Vojtěch Holub. 2012.  Ensemble classifiers for steganalysis of digital media.  IEEE Transactions on Information Forensics and Security 7(2):432–444.
  * Martin and Jurafsky (2000) James H Martin and Daniel Jurafsky. 2000.  Speech and language processing.  International Edition 710\.
  * Mikolov et al. (2010) Tomas Mikolov, Martin Karafiát, Lukas Burget, Jan Cernockỳ, and Sanjeev Khudanpur. 2010.  Recurrent neural network based language model.  In Interspeech. volume 2, page 3.
  * Mikolov et al. (2013) Tomas Mikolov, Ilya Sutskever, Kai Chen, Greg S Corrado, and Jeff Dean. 2013.  Distributed representations of words and phrases and their compositionality.  In Advances in neural information processing systems. pages 3111–3119.
  * Pascanu et al. (2013) Razvan Pascanu, Tomas Mikolov, and Yoshua Bengio. 2013.  On the difficulty of training recurrent neural networks.  ICML (3) 28:1310–1318.
  * Pennington et al. (2014) Jeffrey Pennington, Richard Socher, and Christopher D Manning. 2014.  Glove: Global vectors for word representation.  In EMNLP. volume 14, pages 1532–1543.
  * Révész (2014) Pál Révész. 2014.  The laws of large numbers, volume 4.  Academic Press.
  * Safaka et al. (2016) Iris Safaka, Christina Fragouli, and Katerina Argyraki. 2016.  Matryoshka: Hiding secret communication in plain sight.  In 6th USENIX Workshop on Free and Open Communications on the Internet (FOCI 16). USENIX Association.
  * Salomon (2003) David Salomon. 2003.  Data privacy and security: encryption and information hiding.  Springer Science & Business Media.
  * Srivastava et al. (2014) Nitish Srivastava, Geoffrey E Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan Salakhutdinov. 2014.  Dropout: a simple way to prevent neural networks from overfitting.  Journal of Machine Learning Research 15(1):1929–1958.
  * Sutskever et al. (2011) Ilya Sutskever, James Martens, and Geoffrey E Hinton. 2011.  Generating text with recurrent neural networks.  In Proceedings of the 28th International Conference on Machine Learning (ICML-11). pages 1017–1024.
  * Topkara et al. (2006) Umut Topkara, Mercan Topkara, and Mikhail J Atallah. 2006.  The hiding virtues of ambiguity: quantifiably resilient watermarking of natural language text through synonym substitutions.  In Proceedings of the 8th workshop on Multimedia and security. ACM, pages 164–174.
  * Volkhonskiy et al. (2017) Denis Volkhonskiy, Ivan Nazarov, Boris Borisenko, and Evgeny Burnaev. 2017.  Steganographic generative adversarial networks.  arXiv preprint arXiv:1703.05502 .
  * Wilson et al. (2015) Alex Wilson, Phil Blunsom, and Andrew Ker. 2015.  Detection of steganographic techniques on twitter.  In EMNLP. pages 2564–2569.
  * Wilson et al. (2014) Alex Wilson, Phil Blunsom, and Andrew D Ker. 2014.  Linguistic steganography on twitter: hierarchical language modeling with manual interaction.  In IS&T/SPIE Electronic Imaging. International Society for Optics and Photonics, pages 902803–902803.
  * Wilson and Ker (2016) Alex Wilson and Andrew D Ker. 2016.  Avoiding detection on twitter: embedding strategies for linguistic steganography.  Electronic Imaging 2016(8):1–9.
  * Zhou et al. (2007) Yingjie Zhou, Mark Goldberg, Malik Magdon-Ismail, and Al Wallace. 2007.  Strategies for cleaning organizational emails with an application to enron email dataset.  In 5th Conf. of North American Association for Computational Social and Organizational Science.


# Statistical mechanics of one-dimensional quantum droplets

T. Mithun  Department of Mathematics and Statistics, University of
Massachusetts, Amherst MA 01003-4515, USA  S. I. Mistakidis  Center for
Optical Quantum Technologies, Department of Physics, University of Hamburg,
Luruper Chaussee 149, 22761 Hamburg Germany  P. Schmelcher  Center for Optical
Quantum Technologies, Department of Physics, University of Hamburg, Luruper
Chaussee 149, 22761 Hamburg Germany  The Hamburg Centre for Ultrafast Imaging,
University of Hamburg, Luruper Chaussee 149, 22761 Hamburg, Germany  P. G.
Kevrekidis  Department of Mathematics and Statistics, University of
Massachusetts, Amherst MA 01003-4515, USA

###### Abstract

We study the statistical mechanics and the dynamical relaxation process of
modulationally unstable one-dimensional quantum droplets described by a
modified Gross-Pitaevskii equation. To determine the classical partition
function thereof, we leverage the semi-analytical transfer integral operator
(TIO) technique. The latter predicts a distribution of the observed wave
function amplitudes and yields two-point correlation functions providing
insights into the emergent dynamics involving quantum droplets. We compare the
ensuing TIO results with the probability distributions obtained at large times
of the modulationally unstable dynamics as well as with the equilibrium
properties of a suitably constructed Langevin dynamics. We find that the
instability leads to the spontaneous formation of quantum droplets featuring
multiple collisions and by which are found to coalesce at large evolution
times. Our results from the distinct methodologies are in good agreement aside
from the case of low temperatures in the special limit where the droplet
widens. In this limit, the distribution acquires a pronounced bimodal
character, exhibiting a deviation between the TIO solution and the Langevin
dynamics.

##  I Introduction

The celebrated Gross-Pitaveskii model has proved particularly successful for
studying and describing a variety of macroscopic many-body phenomena in zero-
temperature Bose-Einstein condensates (BECs) Pethick and Smith (2008);
Stringari and Pitaevskii (2003); Kevrekidis  _et al._ (2015). On the other
hand, more recently, a new type of matter wave has emerged Petrov (2015);
Petrov and Astrakharchik (2016), namely the so-called quantum droplets. The
theoretical basis is a two-component (binary) BECs, with intra-component self-
repulsion, yet also inter-component attraction lying in the vicinity of the
self-repulsion. Here, the famous Lee-Huang-Yang quantum correction Lee  _et
al._ (1957) comes into play to account for the averaged effect of quantum
fluctuations beyond the mean-field description. It competes with the mean-
field effects Pethick and Smith (2008); Stringari and Pitaevskii (2003) and
can prevent the possibility of BEC collapse present in the mean-field realm.
Notice that this stabilization mechanism depends crucially on the system’s
dimensionality; namely beyond mean-field fluctuations are attractive in one-
dimension (1D) while being repulsive in higher spatial dimensions.
Importantly, these predictions led to a sequence of experimental realizations
of such quantum droplets Cabrera  _et al._ (2018); Cheiney  _et al._ (2018);
Semeghini  _et al._ (2018); Ferioli  _et al._ (2019); D’Errico  _et al._
(2019), firstly observed in dipolar gases Chomaz  _et al._ (2016); Ferrier-
Barbut  _et al._ (2016). Accordingly, these states constitute nowdays an
emerging topic that is under intensive research in the context of BECs.

Along this vein of experimental developments, it is remarkable that collisions
of such droplet patterns have been experimentally observed recently (leading
to mergers for slow collisions and quasi-elastic passage for fast ones)
Ferioli  _et al._ (2019). Additionally and while the above examples focused on
39K droplets, heteronuclear binary BECs of 87Rb and 41K showcased very stable
quantum droplets on time scales of the order of a second D’Errico  _et al._
(2019). In parallel to this steady stream of experimental demonstrations,
theoretical studies have spearheaded a number of parallel directions. These
include but are not limited to the exploration of self-evaporation dynamics of
droplets Ferioli  _et al._ (2020); Fort and Modugno (2021), quantum droplets
with intrinsic vorticity Li  _et al._ (2018), the impact of discreteness in
the form of semidiscrete droplets with or without topological charge Zhang
_et al._ (2019), or the investigation of 3D stable generalizations of such
states Kartashov  _et al._ (2018). Many of these developments have been
recently summarized in Ref. Luo  _et al._ (2021). Interestingly, the 1D
dynamical features of droplet states are far less explored and are currently
mainly restricted to inelastic collisional aspects of these configurations
especially for high momenta and flat-top droplets Astrakharchik and Malomed
(2018) or their spontaneous generation due to the modulation instability (MI)
Mithun  _et al._ (2020). Additionally, the phase diagram of quantum droplets
trapped in 1D optical lattices was studied very recently Morera  _et al._
(2020).

The above efforts have mainly focused on the dynamical aspects of the
droplets, and principally so at the zero temperature setting. However,
naturally, studying the finite temperature dynamics of BEC systems is a topic
of broad theoretical and experimental appeal Griffin  _et al._ (2009). It is
indeed challenging to completely eliminate thermal effects in current
experiments, a fact that further justifies the usage of models that can
operate at finite temperatures De Rosi  _et al._ (2021). This is particularly
intriguing in 1D where the role of quantum fluctuations, being inherently
related to droplet formation, is more prominent Parisi  _et al._ (2019);
Parisi and Giorgini (2020); Mistakidis  _et al._ (2021) and three-body losses
are suppressed Astrakharchik and Giorgini (2006) compared to higher dimensions
while the experimental probe of these configurations is still elusive. One of
the commonly used methods to include the temperature induced fluctuations is
the so-called truncated Wigner method Blakie  _et al._ (2008); Sreedharan  _et
al._ (2020). Recently, a considerable amount of attention has been drawn
towards the so-called positive P-method too Deuar and Drummond (2006a, b). In
general, the final task of these methods is solving a stochastic equation that
incorporates the thermal fluctuations to the corresponding Hamiltonian
dynamical system. In statistical mechanics, a principal task consists of the
evaluation of the partition function corresponding to the energy functional of
the system at hand. A central tool to this effect in 1D settings consists of
the so-called transfer integral operator (TIO) method Scalapino  _et al._
(1972). Indeed, this method renders the above problem equivalent to solving
the single particle Schrödinger equation. Moreover, it has been shown that a
stochastic method represented by a suitable Langevin equation is comparable
with the solution of TIO Khare  _et al._ (1997); Nunnenkamp  _et al._ (2007);
Habib and Lythe (2000). Further, recently it has been shown that the steady
state of molecular dynamics and the classical field method also reproduces the
TIO results Nunnenkamp  _et al._ (2007); Aarts  _et al._ (2000).

To the best of our knowledge, such techniques motivated by statistical
mechanics, while being successful at the level of single-component Gross-
Pitaevskii model Nunnenkamp  _et al._ (2007), have yet to be explored in the
quantum droplet realm of effectively multi-component systems. Indeed, this is
a central focus of the present work. More concretely, we study the statistical
mechanics of a quantum droplet in a 1D ring geometry. One of the peculiarities
of a quantum droplet is its incompressibility leading to a maximal critical
density for sufficiently large particle numbers.

This relevant system, as mentioned above, can be represented with a modified
binary Gross-Pitaveskii (MGP) equation, where the form of a Lee-Huang-Yang
term is determined by the dimension of the system Petrov and Astrakharchik
(2016); Astrakharchik and Malomed (2018); Tylutki  _et al._ (2020); Lavoine
and Bourdel (2021). For the case of particle-balanced components with equal
masses and equal intra-component interaction strengths, which will be of
primary interest herein, the binary MGP equation reduces to a single-component
one. We determine the classical partition function corresponding to our model
by mapping the functional integration of the partition function to a single-
particle Schrödinger equation via the TIO technique. Note that the classical
partition function corresponds to the high-temperature limit of the quantum
partition function describing an interacting Bose gas Pippard (1964). We then
numerically verify the predictions of TIO via a suitably crafted Langevin
dynamics Khare  _et al._ (1997); Nunnenkamp  _et al._ (2007). Our results
indicate that the equilibrium properties of the Langevin dynamics are well in
line with the TIO solution at intermediate and high temperatures for
$\mu\rightarrow\mu_{0}$, the limit where the droplets tend to disappear; the
agreement is less adequate at low temperatures. We also compare the TIO
results with the long-time dynamical evolution of the original droplet system
in the regime where the latter falls into MI Nguyen  _et al._ (2017) and
subsequently relaxes. The recent studies on the statistical properties of MI
motivate this analysis Gelash  _et al._ (2019); Kraych  _et al._ (2019).
Indeed, in Ref. Mithun  _et al._ (2020), it has been shown that the MI
resulting from the small perturbation of a plane-wave state leads to the
formation of quantum droplet structures that undergo inelastic collisions.
Interestingly, we find that, as a result of the inelastic collisions these
generated droplet structures coalesce and their equilibrium properties are
well matched by the TIO analysis. It is worthwhile to mention here that within
our study we neglect any possible temperature dependence of the Lee-Huang-Yang
term, a prospect that remains still an open question.

The work is organized as follows. We introduce the MGP model in Sec. II and
then determine the classical partition function by using the TIO method in
Sec. III. Section IV is devoted to developing and exploring the relevant
Langevin dynamics. We report the results of the MI dynamics within the MGP
framework in Sec. V and provide an outlook in Sec. VI. Appendix A elaborates
on the extraction of the single-particle Schrödinger picture through the
transfer integral method. Appendix B explicates the derivation of the Langevin
equations in interaction imbalanced two-component mixtures. Finally in
Appendix C it is shown that the inclusion of thermal effects in the mixture,
emulated by a dissipative term in the MGP, leads to a more pronounced
coalescence process of the droplets associated with the suppression of the
background fluctuations.

##  II The modified Gross-Pitaevskii framework

The starting point of our analysis will be the dimensionless MGP equation
describing 1D quantum droplets emerging in symmetric binary mixtures Petrov
and Astrakharchik (2016); Mithun  _et al._ (2020). In particular,

$i\frac{\partial\psi}{\partial t}=-\frac{1}{2}\frac{\partial^{2}\psi}{\partial
z^{2}}+|\psi|^{2}\psi-|\psi|\psi-\mu\psi,$ (1)

with the normalization condition

$\int_{-\infty}^{+\infty}n\leavevmode\nobreak\ dz=N,\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ \text{and}\leavevmode\nobreak\
\leavevmode\nobreak\ \leavevmode\nobreak\ n=|\psi(z)|^{2}.$ (2)

Here, $N$ is the number of atoms and $\mu$ represents the chemical potential.
The effective single-component Eq. (1) is obtained as a reduction of the
respective coupled set of two-component modified Gross-Pitaveskii equations
[see also Appendix B] under symmetry considerations. The latter refer to a
mixture with same particle number per species ($N_{1}=N_{1}\equiv N/2$), equal
repulsive intra-component interaction strengths ($g_{1}=g_{2}\equiv g>0$) and
equal masses ($m_{1}=m_{2}\equiv m$), see for details Petrov and Astrakharchik
(2016); Astrakharchik and Malomed (2018). Moreover according to Eq. (1), the
units of length, time and wave function are expressed in terms of the healing
length $\xi$, $\hbar/(m\xi^{2})$, and
$(2\sqrt{g})^{3/2}/\sqrt{\pi\xi}(2|\delta g|)^{3/4}$, respectively, with
$\xi=\frac{\pi\hbar^{2}}{m}\frac{\sqrt{2|\delta g|}}{g^{3/2}}$. Also, $\delta
g=g_{12}+g$ and $g=\sqrt{g_{1}g_{2}}$, where $g>0$ represent the repulsive
intra-component interaction strengths of the symmetric mixture, while
$g_{12}<0$ denotes the inter-component attractive interaction. Additionally,
for the existence of a quantum droplet $\delta g\ll g$ should hold. In a
corresponding experiment, $\delta g$ can be tuned using the Feshbach resonance
technique Cabrera  _et al._ (2018); Cheiney  _et al._ (2018).

For $\delta g/g>0$, Eq. (1) gives rise to an exact localized flat-top (FT)
solution Astrakharchik and Malomed (2018). The latter represents a quantum
droplet which originates from the balance between the effective cubic self-
repulsion and quadratic attraction characterized by the central density
$n_{0}$ and the chemical potential $\mu_{0}$, where

$n_{0}=\frac{4}{9},\leavevmode\nobreak\ \leavevmode\nobreak\
\mu_{0}=-\frac{2}{9},$ (3)

respectively. In what follows we shall focus on the $\delta g/g>0$ regime.
However, for reasons of completeness we remark that solutions to Eq. (1) exist
also for $\delta g/g<0$, see Ref. Mithun  _et al._ (2020) for a more
elaborated discussion.

##  III Statistical Mechanics

To study the statistical mechanics of the model under consideration, as
described by Eq. (1), we establish the corresponding classical partition
function via the TIO method, as indicated in the Introduction in line with
Refs. Nunnenkamp  _et al._ (2007); Scalapino  _et al._ (1972); Khare  _et al._
(1997). For completeness, we will present selected results for the two-
component case in Appendix B. In order to proceed, we consider the droplet on
a periodic domain (i.e., effectively a ring) of length $L$. Then, the form of
the classical partition function reads

$Z=\int D(\psi,\psi^{\ast})e^{-\beta F[\psi,\psi^{\ast}]},$ (4)

where

$\begin{split}F&=\int dz(H-\mu N)\end{split}$ (5)

is the free energy, $\beta$ is the inverse temperature, $H$ represents the
Hamiltonian corresponding to Eq. (1) and $N$ is the total number of particles.
In accordance with the TIO methodology, the functional integration can be
reduced to an eigenvalue problem Scalapino  _et al._ (1972); Trullinger and
DeLeonardis (1980). For the Eq. (1), we arrive at the corresponding eigenvalue
equation (see Appendix A for a detailed derivation)

$\Bigg{[}-\frac{1}{2\beta^{2}}\frac{\delta^{2}}{\delta\psi^{2}}+V_{d}(\psi)\Bigg{]}\phi_{n}(\psi)=E_{n}\phi_{n}(\psi).$
(6)

Evidently, Eq. (6) corresponds to a 1D single-particle Schrödinger equation.
In this expression, $\phi_{n}$ represent the corresponding single-particle
eigenfunctions, being functionals of the $\psi$ field, while $E_{n}$ denote
the eigenvalues. Importantly, the effective TIO potential appearing in Eq.
(6), possesses the form

$V_{d}(\psi)=\frac{1}{2}|\psi|^{4}-\frac{2}{3}|\psi|^{3}-\mu|\psi|^{2}.$ (7)

Notably, this represents the common mean-field potential of the Gross-
Pitaevskii theory Petrov and Astrakharchik (2016). The shape of these
effective potentials in the different parameter regions of the system is
crucial in order to understand where bound state solutions, and thus droplet-
like configurations, are prone to appear. For this reason, we next show
$V_{d}(\psi=u)$ in Fig. 1 since the one-component case will be our focus in
the following. According to Eq. (7) the minimum of $V_{d}(\psi=u)$ is given by
$u=(3b+\sqrt{9b^{2}+32a\mu})/(8a)$ for $\mu>0$, while
$u=(3b\pm\sqrt{9b^{2}+32a\mu})/(8a)$ for $\mu\leq 0$ where $a=\frac{1}{2}$ and
$b=\frac{2}{3}$. As a result for $\mu<0$ the effective potential features a
double-well structure, see Fig. 1. This feature is related to the presence of
a bimodal probability distribution as it will be argued later on [Eq. (8)].

Figure 1: Effective anharmonic potential $V_{d}(u)$, described by Eq. (7).
Shown are different values of the chemical potential $\mu$ (see legend) with
$\mu_{0}$ given by Eq. (3).

We solve the Eq. (6) numerically to determine the corresponding eigenvalues
and eigenvectors utilizing exact diagonalization. Notice that this is the only
aspect that renders the TIO approach semi-analytical. The main contribution to
the partition function in Eq. (4) will be from the lowest TIO eigenvalue. The
thermodynamic properties of an equilibrium system can be well understood by
resorting to the underlying probability distribution of amplitudes (PDA),
$P(|\psi|=u)$. The latter can be expressed in terms of the eigenfunction
corresponding to the lowest eigenvalue $E_{0}$ of Eq. (6) Nunnenkamp  _et al._
(2007). In particular $P(|\psi|=u)$ acquires the form

$P(|\psi|=u)=2u|\phi_{0}(u)|^{2}.$ (8)

Moreover, the steady state properties of a thermodynamic system can be further
characterized by employing the two-point spatial correlation function $C(z)$
Naraschewski and Glauber (1999); Katsimiga  _et al._ (2017). Here, $z$ denotes
the relative distance between two distinct spatial locations. This correlation
function can be written with respect to the eigenvalues and eigenfunctions of
Eq. (6) as

$\begin{split}C(z)&=\langle\psi(.)\psi(.+z)\rangle\\\ &=\sum_{n}\lvert\int
du\phi_{n}^{\ast}(u)u\phi_{0}(u)\rvert^{2}e^{-\beta|z|(E_{n}-E_{0})}.\end{split}$
(9)

It provides a measure of the coherence among the particles Naraschewski and
Glauber (1999); Mistakidis  _et al._ (2018). Concretely, it is bounded from
above and below taking as a maximal and minimal value the considered particle
number and zero respectively. If $C(z)$ is maximal for every $z$ then the
system is termed fully coherent and it is characterized by quasi long-range
order. Otherwise losses of coherence occur. In the thermodynamic equilibrium,
the long-range order is expected to vanish and accordingly $C(z)$ tends to
zero for increasing $z$ Nandkishore and Huse (2015). We finally remark that,
for the settings considered herein, the contributions of the (higher) excited
states of the effective potentials (described by Eqs. (8) and (9)) are
neglected, as being exponentially smaller than the dominant lowest-lying
states.

##  IV Langevin dynamics

It has been demonstrated that, for any positive temperature, the exact results
obtained from the TIO in a BEC system can be compared with that of a Langevin
equation with a Gaussian additive white noise term Nunnenkamp  _et al._
(2007). This type of stochastic dynamics contributes towards relaxing the
system configuration to the free-energy minimum, while at the same time
accounting for the thermal fluctuations arising around this minimum Parisi
_et al._ (1981). Indeed, the point of the Langevin dynamics is to allow the
fields to evolve in a way such that they eventually relax to the appropriate
equilibrium distribution and then to sample the suitable observables. I.e.,
this dynamics relaxes the field configuration via fictitious time-evolution to
the relevant minimum of the free energy surface, while the white noise
randomly drives the field. In the following, we consider the Langevin
equations corresponding to Eq. (1), in line also with the earlier work of
Nunnenkamp  _et al._ (2007) for regular one-component BECs

$\frac{\partial\psi}{\partial
t}=-\Big{[}-\frac{1}{2}\frac{\partial^{2}}{\partial
z^{2}}+|\psi|^{2}-|\psi|-\mu\Big{]}\psi+\xi_{1}(z,t).$ (10)

The term $\xi_{1}$ represents the white noise having a Gaussian distribution
with the correlation

$\langle\xi_{1}^{\ast}(z,t)\xi_{1}(z^{\prime},t^{\prime})\rangle=\frac{2}{\beta}\delta(z-z^{\prime})\delta(t-t^{\prime}).$
(11)

The Langevin equations are indeed found to relax to the equilibrium state at
large evolution times. This is caused by the fact that the time average of the
spatio-temporal correlation relaxes to its equilibrium spatial correlation or
in other words the time average of the noise approaches its equilibrium
distribution Hohenberg and Halperin (1977). We solve Eq. (10) numerically by
using the xmds package Dennis  _et al._ (2013).

Figure 2: Probability distribution function obtained from TIO (solid lines)
and within the Langevin dynamics for different values of the inverse
temperature $\beta$. Namely, $\beta=0.2$ (green crosses), $\beta=2$ (blue
circles) and $\beta=95$ (red squares). The chemical potential used is
$\mu=\mu_{0}+0.00001$. Figure 3: Probability distribution function
$P(|\psi|=u)$ at different total evolution times $t_{F}$ (see legend) for
$\beta=95$ obtained within the Langevin dynamics. A saturation of
$P(|\psi|=u)$ is observed for $t_{F}>200$. The chemical potential used is
$\mu=\mu_{0}+0.00001$.

Note that for the numerical simulations of Eq. (10), to be presented below, a
sample of 1000 trajectories is found to be sufficient for the convergence of
the relevant probability distributions, with the domain size being $L=30$. We
first fix the chemical potential $\mu=\mu_{0}+0.00001$, where the exact
solution of Eq. (1) is a FT shaped droplet. Figure 2 depicts the probability
distribution $P(|\psi|=u)$ results of the Langevin dynamics for (long) total
evolution time $t_{F}\geq 200$. The green (crosses), blue (circles) and red
(squares) symbols represent the numerically obtained steady states for
$\beta=0.2$, $\beta=2$ and $\beta=95$, respectively; the solid lines represent
the corresponding TIO solutions. Inspecting $P(|\psi|=u)$ [Fig. 2] we can
readily deduce that as the temperature decreases (from $\beta=0.2$ to
$\beta=95$), it changes from a single peak distribution (unimodal) to a
bimodal one in terms of the TIO prediction. The comparison between the TIO
results and the ones within the Langevin approach highlights that the steady
state of the Langevin dynamics is in very good agreement with the TIO
solutions both at higher and intermediate temperatures. On the other hand, at
the low temperature regime, the steady state of the Langevin dynamics deviates
from the TIO solution. The latter case is the more numerically delicate one
due to the nature of the relevant effective TIO potential, a feature to which
we attribute the observed discrepancy. In order to confirm that this
discrepancy is not due to the considered total evolution time $t_{F}$ that
allows to reach the steady state, we determine the $P(|\psi|=u)$ for different
$t_{F}$ at $\beta=95$, see in particular Fig. 3. Interestingly, we observe
that the probability distribution for evolution times in the range of
$t_{F}\in(160-240$) develops a bimodal structure, being proximal to the
expected picture from the TIO analysis. However, upon considering larger
timescales the relative difference in the amplitude between the humps
appearing in $P(|\psi|=u)$ tends to be suppressed.

We will further address this discrepancy when considering the dynamics of the
full model, Eq. (1), (rather than the modified Langevin one) in the following
section. To gain additional insight on the dependence of the steady state
distribution at low temperatures, we now vary the chemical potential $\mu$ by
fixing the inverse temperature parameter $\beta=95$. The obtained results for
the $P(|\psi|=u)$ are provided in Fig. 4. It shows that the TIO solution
develops a bimodal distribution only when $\mu\rightarrow\mu_{0}$. As $\mu$
deviates from this FT droplet limit, the TIO solution leads to a unimodal
distribution and the results of the Langevin dynamics match well with these
solutions. This once again suggests that presumably the Langevin dynamics is
not able to capture the delicate FT configuration of the droplet at the
intermediate temperature limit, as we will argue further below. Nevertheless,
in all other settings, the TIO semi-analytical prediction is adequately
captured by the modified Langevin dynamics.

Figure 4: Probability distribution $P(|\psi|=u)$ obtained from the Langevin
dynamics for $\mu=\mu_{0}+0.00001$ (red crosses), $\mu=\mu_{0}+0.005$ (green
circles), $\mu=\mu_{0}+0.01$ (blue squares), and $\mu=\mu_{0}+0.05$ (black
triangles). The solid lines represent the corresponding TIO solutions. The
inverse temperature is $\beta=95$.

To gain further insight into the emergent steady state of the droplets, we
additionally determine the two-point correlation function $C(|z|)$ [Eq. (9)].
The latter is provided in Fig. 5 for sufficiently large evolution times, i.e.,
$t_{F}\geq 200$. This observable allows to assess the underlying coherence
losses which are induced here by the temperature. It evinces that the
predictions of the TIO are in accordance with the ones of the Langevin
dynamics for both $\beta=0.2$ and $\beta=2$ inverse temperatures [see lower
and upper panels of Fig. 5]. Particularly, quasi long-range order is
suppressed and therefore coherence is lost since $C(|z|)\to 0$ for increasing
$z$. This fact also supports the appearance of a steady state Nandkishore and
Huse (2015) for the droplet. On the other hand, at the large inverse
temperature ($\beta=95$), in line with the mismatch in the probability
distributions, the behavior of the correlation function deviates from the
respective transfer integral prediction in the $\mu\rightarrow\mu_{0}$ limit
as seen in the lower panel of Fig. 5. Concretely, the TIO method shows a
vanishing coherence for larger values of $z$, whilst the Langevin dynamics
predicts that the coherence of the system is maintained.

Figure 5: Correlation function $C(z)$ obtained from the TIO method (solid
lines) and from the Langevin dynamics (dots) for different temperatures.
Namely, $\beta=0.2$ (top panel), $\beta=2$ (middle panel) and $\beta=95$
(bottom panel) at $\mu=\mu_{0}+0.00001$. Here the considered total evolution
time is $t_{F}\geq 200$.

##  V MI dynamics in the modified Gross-Pitaevskii approach

Recently, it has been shown that the dynamics of a classical field method
(CFM) adequately traces the features of the TIO solution for the equilibrium
statistical mechanics of the single-component BEC case Nunnenkamp  _et al._
(2007). Indeed, a key point of that work is the comparison of three different
methods: the Langevin approach detailed above with a molecular dynamics (MD)
one Aarts  _et al._ (2000) (based on Hamiltonian dynamics of Klein-Gordon
type) and the CFM. Considering the MD methodology as less directly related to
the system at hand, we focus here on the CFM as a complement of the Langevin
approach presented above. The core idea of the CFM is to obtain a steady state
where the dynamics is governed by the Gross-Pitaevskii equation and then
compare this steady state solution with the TIO solution. In this method, the
initial state is a non-equlibirium one where only a few modes in momentum
space are excited. Moreover, it is well known that in such systems the MI
phenomenon naturally creates non-equilibrium conditions when subjected to a
small perturbation (seeding the relevant instability) Nguyen  _et al._ (2017);
Everitt  _et al._ (2017). The Gross-Pitaevskii equation with an effective
attractive interaction is known to be modulationally unstable. The MI dynamics
of the model as described by Eq. (1) has been recently studied in detail
Mithun  _et al._ (2020) and argued to lead to droplet nucleation. Motivated by
this finding we shall subsequently focus on the dynamical response of the
symmetric mixture generated by its direct time-evolution within the MGP and
being associated with the unstable dynamics stemming from the MI. In
particular, we aim to unravel the spontaneous creation of droplets due to MI
and their fate in the long-time dynamics. Indeed, we allow the system to
evolve into its asymptotic (equilibrium) state from an initially dynamically
unstable initial condition and “observe” the resulting effective temperature
by comparison with the TIO solution results.

Figure 6: MI dynamics for the system parameter values $n\approx 0.11$ and
$\mu=\mu_{0}+0.00001$. (a) Spatiotemporal evolution of $|\psi(z,t)|$. (b)
Probability distribution function $P(|\psi|=u)$ and (c) correlation function
in momentum space $\langle|\psi(k)|^{2}\rangle$ at a final evolution time of
$t_{F}=4000$. The results shown in (b) and (c) are averaged ones over a sample
of twelve different initial conditions and the error bars shown in (b)
represent the respective standard deviation. The green solid line represents
the TIO solution at $\beta=95$, while, the blue dotted points are obtained
through the direct simulation of the MGP as described by Eq. (1).

As an initial condition we consider a homogeneous particle density
distribution subjected to a weak amplitude perturbation of the form
$\delta\psi=A_{p}e^{i\theta_{r}}$. The strength $A_{p}$ of the perturbation is
of the order of $10^{-8}$ and $\theta_{r}$ represents the uniform noise
distribution taking values in the interval $[0,\pi]$ Mithun  _et al._ (2020).
We fix $\mu=\mu_{0}+0.00001$ that corresponds to $n\approx 0.11$. The results
are presented in Fig. 6. The MI initially leads to the formation of small
droplets as discussed in Mithun  _et al._ (2020). In the course of the
evolution these small droplets undergo inelastic collisions (analogously to
what was observed in Ferioli  _et al._ (2019)) and feature a coarsening stage.
This results in the formation of “large” droplets as shown in Fig. 6(a). Also,
we noted that the amplitude of these waveforms is close to the central
density, $n_{0}$ of the FT shaped droplet. The corresponding PDA is depicted
in Fig. 6(b) (blue dotted points). Interestingly, the PDA exhibits a bimodal
distribution, where the peak at large amplitude ($u\approx 0.66$) represents
the large amplitude droplets, while peaks at small amplitude ($u\approx 0.33$)
are related to small-amplitude phononic excitations Tylutki  _et al._ (2020)
that are widespread within the spatial extension of the cloud.

To compare this equilibrium state with the TIO solution, we considered
different values of $\beta$ and obtained that, e.g., for $\beta=95$ the
probability distribution functions are in close correspondence with each other
as shown in Fig. 6(b) with the green solid line. Interestingly, the accordance
between the PDA and TIO solution predictions holds also for the corresponding
correlation function in momentum space Mistakidis  _et al._ (2015) i.e.
$\braket{}{\psi(k)}{{}^{2}}=(1/2\pi)\int dzC(z)e^{-ikz}$. The latter is
routinely accessible in current ultracold atom experiments via time-of-flight
imaging Bloch  _et al._ (2008). This observable is provided in Fig. 6(c) at
$t_{F}=4000$ and $\beta=95$. As can be seen, it features a central peak
structure around $k=0$ signaling the steady state of the droplet while its
shape is in very good agreement between the two approaches. To further
corroborate that the same argumentation holds also for other values of the
chemical potential, we exemplarily showcase in Fig. 7 the cases of
$\mu=\mu_{0}+0.01$ and $\mu=\mu_{0}+0.1$. It becomes apparent that even for
these values of the chemical potential, the final state resulting from the MI
dynamics is a large excitation on top of a fluctuating background bearing
small amplitudes. Hence, the corresponding PDA is a bimodal distribution.
However, the weight of the small amplitude excitations in the corresponding
PDA distribution decreases with increasing $\mu$, in line with the respective
TIO prediction shown in Fig. 4.

In our considerations so far, we have limited the study of an effective
temperature and chemical potential associated with the free energy and atom
number of our droplet forming system. In Appendix C, we also briefly consider
a variant of the MGP dynamics involving dissipative perturbations, as a first
glimpse towards the role of finite temperature effects in the spontaneous
generation of droplets. In particular, we explicate that the above-discussed
coalescence process of droplets becomes more prominent in the presence of
thermal effects of the background as captured by the direct inclusion of a
damping term in the modified Gross-Pitaevskii equation.

Figure 7: Time-evolution of $|\psi(z,t)|$ explicating the MI dynamics, the
subsequent generation of droplets and their coalescence while (a)
$\mu=\mu_{0}+0.01$, and (b) $\mu=\mu_{0}+0.1$.

##  VI Conclusions and Future Challenges

In the present work, we have sought to address the equilibrium statistical
mechanics of one-dimensional quantum droplets. This system is represented by a
modified two-component Gross-Pitaevskii equation and under the symmetry
considerations (equal wave functions and equal intra-component interaction
strengths and masses) it reduces to a single-component equation. To determine
the classical partition function, we first reduce the functional integration
of the partition function to a single-particle Schrödinger problem by using
the TIO technique. The probability distribution of the amplitudes obtained
from TIO shows unimodality in the limit $\mu\rightarrow\mu_{0}$ at
intermediate and large temperatures. However, at small temperatures the
relevant probability exhibits a bimodal shape. On the other hand, for the
values of $\mu$ far from $\mu_{0}$, the probability distribution is always
found to be unimodal.

We compared the results of the TIO analysis with the equilibrium properties of
a suitably crafted Langevin dynamics. The equilibrium properties of the
Langevin dynamics reproduce well the corresponding ones of the TIO solutions
except in the limit $\mu\rightarrow\mu_{0}$ for low temperatures. In this
latter scenario, a bimodal distribution is featured showing in turn a partial
discrepancy between the equilibrium properties of the Langevin dynamics and
the TIO solutions. Additionally, we compared the probability distribution of
the modulational instability dynamics at large times with the TIO findings,
obtaining good agreement in the cases under consideration for a suitable
choice of the temperature. Importantly, we have shown that the modulational
instability inherent to the system gives rise to spontaneous quantum droplet
nucleation. These in turn undergo a coarsening stage experiencing collision
events and leading to their coalescence.

Very recently, significant attention has been drawn to the thermodynamics of
quantum droplets De Rosi  _et al._ (2021). Beyond the realm of the present
study, several topics remain open for future studies. An extension of
thermodynamic considerations in higher dimensional settings would be a topic
of particular interest for a variety of reasons. On the one hand, semi-
analytical tools such as the TIO do not straightforwardly generalize to higher
dimensions, hence different theoretical approaches would need to be brought to
bear. Moreover, even numerically the nonlinearities of the system being
logarithmic in 2D and featuring a different form than the one considered
herein ($\propto|\psi|^{3}\psi$) in 3D could yield different outcomes as
regards the probability distributions and the correlation functions. On the
other hand, the genuinely multi-component case where the two components are
not forced to be equal and its systematic consideration and direct comparison
e.g. with a variational approach Mistakidis  _et al._ (2018) is also of
interest for further study.

## Acknowledgments

We thank S. Flach and K. Burnett for fruitful comments and discussions. S. I.
M. gratefully acknowledges financial support in the framework of the Lenz-
Ising Award of the University of Hamburg. The present paper is based on work
that was supported by the US National Science Foundation under DMS-1809074 and
PHY-2110030 (PGK).

##  Appendix A Derivation of the transfer integral problem

Here, we briefly discuss the derivation of the single-particle Schrödinger Eq.
(6), describing the equilibrium state of the symmetric mixture, from the
underlying classical partition function of Eq. (4) Scalapino  _et al._ (1972).
For a ring of length $L=M\Delta z$, where $M$ represents the number of grid
points and $\Delta z$ denotes the (quite fine) spatial discretization used in
the numerical implementation of the problem, we can discretize the free energy
functional introduced in Eq. (5) as

$\begin{split}F&=\Delta z\sum_{i=1}^{M}f(\psi_{i},\psi_{i+1})=\Delta
z\sum_{i=1}^{M}\Big{[}\frac{1}{2}|\frac{\psi_{i+1}-\psi_{i}}{\Delta z}|^{2}\\\
&+\frac{1}{2}n_{i}^{2}-\frac{2}{3}n_{i}^{3/2}-\mu n_{i}\Big{]},\end{split}$
(12)

and we obtain the partition function Scalapino  _et al._ (1972)

$Z=\prod_{i=1}^{M}\int_{-\infty}^{\infty}d\tilde{\psi}_{1}^{{}^{\prime}}d\tilde{\psi}_{i}e^{-\beta\Delta
zf(\psi_{i},\psi_{i+1})}\delta(\tilde{\psi}_{1}-\tilde{\psi}_{1}^{{}^{\prime}}),$
(13)

where $d\tilde{\psi}_{i}=\sqrt{\frac{\beta}{2\pi\Delta z}}d\psi_{i}$ and $M$
is the number of segments within the ring of width $\Delta z$. We now expand
the $\delta$ function in terms of a normalized set of (complete)
eigenfunctions and obtain

$Z=\sum_{n}\prod_{i=1}^{M}\int_{-\infty}^{\infty}d\tilde{\psi}_{1}^{{}^{\prime}}d\tilde{\psi}_{i}\phi_{n}(\tilde{\psi}_{1}^{{}^{\prime}})e^{-\beta\Delta
zf(\psi_{i},\psi_{i+1})}\phi_{n}(\tilde{\psi}_{1}).$ (14)

We can calculate this partition function from the eigenvalue problem

$\int_{-\infty}^{\infty}d\tilde{\psi}_{i}e^{-\beta\Delta
zf(\psi_{i},\psi_{i+1})}\phi_{n}(\tilde{\psi}_{i})=e^{-\beta\Delta
zE_{n}}\phi_{n}(\tilde{\psi}_{i+1}),$ (15)

where $\phi_{n}$ and $E_{n}$ are the eigenfunctions and eigenvalues of the
transfer matrix equation. We then reduce this equation to the single-particle
Schrödinger problem as follows. Performing the Taylor series expansion of the
eigenfunction $\phi_{n}(\tilde{\psi}_{i})$ around $(\tilde{\psi}_{i+1})$ we
arrive at

$\begin{split}\int_{-\infty}^{\infty}d\tilde{\psi}_{i}e^{-\beta\Delta
zf(\psi_{i},\psi_{i+1})}\phi_{n}(\tilde{\psi}_{i})=e^{-\beta\Delta
z\Big{[}\frac{1}{2}\frac{\delta
g}{g}n^{2}-\frac{1}{2}\frac{2^{5/2}}{3\pi}n^{3/2}-\mu
n\Big{]}}&\times\Big{(}1+\frac{\Delta
z}{2\beta}\frac{\partial^{2}}{\partial\psi_{i+1}^{2}}\Big{)}\phi_{n}(\tilde{\psi}_{i+1})\\\
&=e^{-\beta\Delta zH}\phi_{n}(\tilde{\psi}_{i+1})\\\ &\equiv e^{-\beta\Delta
zE_{n}}\phi_{n}(\tilde{\psi}_{i+1})\end{split}$ (16)

where $\Big{(}1+\frac{\Delta
z}{2\beta}\frac{\partial^{2}}{\partial\psi_{i+1}^{2}}\Big{)}\approx
e^{\frac{\Delta z}{2\beta}\frac{\partial^{2}}{\partial\psi_{i+1}^{2}}}$.

The effective Hamiltonian reads

$\begin{split}H=-\frac{1}{2\beta^{2}}\frac{\partial^{2}}{\partial\psi^{2}}+\frac{1}{2}n^{2}-\frac{2}{3}n^{3/2}-\mu
n.\end{split}$ (17)

Here, the second term corresponds to the mean-field contribution and the third
one is the leading-order beyond mean-field correction accounting for quantum
fluctuations.

Figure 8: Two-dimensional probability distribution function
$P[{\textrm{R}e}(\psi_{1})=u_{1},{\textrm{R}e}(\psi_{2})=u_{2}]$ obtained from
the Langevin dynamics (left panels) and the TIO approach (right panels) for
$\beta=0.2$ (top panels), $\beta=2$ (middle panels) and $\beta=95$ (bottom
panels). Other parameters used are $\mu=\mu_{0}+0.00001$, $\delta g/g=0.05$,
and $\mathcal{P}=1$. The total evolution time is $t_{F}=240$. Figure 9: Two-
dimensional probability distribution function
$P[Re(\psi_{1})=u_{1},Re(\psi_{2})=u_{2}]$ as predicted from the Langevin
dynamics (left panel) and the TIO method (right panel) for $\beta=2$,
$\mathcal{P}=1.25$ as well as $\mu=\mu_{0}+0.00001$. The considered total
evolution time corresponds to $t_{F}=240$.

##  Appendix B Two-component Langevin equations

Having analyzed in the main text the case of symmetric mixtures, here we
briefly discuss the generalization of our method to interaction imbalanced
two-component settings. For this reason we resort to the dimensionless MGP
equation for a 1D binary quantum droplet Petrov and Astrakharchik (2016);
Mithun  _et al._ (2020),

$\begin{split}i\frac{\partial\psi_{1}}{\partial
t}&=-\frac{1}{2}\frac{\partial^{2}\psi_{1}}{\partial
z^{2}}+\mathcal{Q}(\mathcal{P}+\mathcal{G}\mathcal{P}^{-1})|\psi_{1}|^{2}\psi_{1}-\mathcal{Q}(1-\mathcal{G})|\psi_{2}|^{2}\psi_{1}\\\
&-\mathcal{P}\sqrt{\mathcal{P}|\psi_{1}|^{2}+\mathcal{P}^{-1}|\psi_{2}|^{2}}\psi_{1}-\mu_{1}\psi_{1},\\\
i\frac{\partial\psi_{2}}{\partial
t}&=-\frac{1}{2}\frac{\partial^{2}\psi_{2}}{\partial
z^{2}}+\mathcal{Q}(\mathcal{P}^{-1}+\mathcal{G}\mathcal{P})|\psi_{2}|^{2}\psi_{2}-\mathcal{Q}(1-\mathcal{G})|\psi_{1}|^{2}\psi_{2}\\\
&-\frac{1}{\mathcal{P}}\sqrt{\mathcal{P}^{-1}|\psi_{2}|^{2}+\mathcal{P}|\psi_{1}|^{2}}\psi_{2}-\mu_{2}\psi_{2},\end{split}$
(18)

where the involved parameters stand for

$\begin{split}\mathcal{P}\equiv\sqrt{\frac{g_{1}}{g_{2}}},&\leavevmode\nobreak\
\mathcal{Q}=\frac{g(1+\mathcal{P})^{2}}{2\mathcal{P}\delta
g},\leavevmode\nobreak\ \mathcal{G}=\frac{2\mathcal{P}^{2}\delta
g}{g(1+\mathcal{P}^{2})^{2}}\\\ &\leavevmode\nobreak\
\text{and}\leavevmode\nobreak\ \delta g=g_{12}+g.\end{split}$ (19)

In the above expressions, $g_{1}>0$ and $g_{2}>0$ correspond to the repulsive
intra-component interaction strengths, while $g_{12}<0$ is the inter-component
attractive interaction. Also, $\mu_{1}$ and $\mu_{2}$ are the chemical
potentials of the individual subsystems. The parameter $\mathcal{P}$
quantifies the intra-component interaction imbalance and $\mathcal{G}$
measures the deviation from the balance point of the mean-field repulsion and
attraction where $\delta g=0$. In the coupled system of Eqs. (18) the units of
length, time and the wave function are expressed in terms of $\xi$,
$\hbar/(m\xi^{2})$, and
$\frac{(\sqrt{g_{1}}+\sqrt{g_{2}})^{3/2}}{\sqrt{\pi\xi}(2|\delta g|)^{3/4}}$,
with the healing length being $\xi=\frac{\pi\hbar^{2}}{m}\frac{\sqrt{2|\delta
g|}}{g(\sqrt{g_{1}}+\sqrt{g_{2}})}$.

For the Eq. (18), we arrive at the corresponding eigenvalue equation

$\begin{split}\Bigg{[}-\frac{1}{2\beta^{2}}\Big{(}\frac{\delta^{2}}{\delta\psi_{1}^{2}}+\frac{\delta^{2}}{\delta\psi_{2}^{2}}\Big{)}&+V_{2d}(\psi_{1},\psi_{2})\Bigg{]}\phi_{n}(\psi_{1},\psi_{2})\\\
&=E_{n}\phi_{n}(\psi_{1},\psi_{2}),\end{split}$ (20)

As it can be seen, Eq. (20) is a two-dimensional single-particle Schrödinger
equation. Also, $\phi_{n}$ refer to the single-particle eigenfunctions, being
functionals of the $\psi$ field, and $E_{n}$ are the respective eigenvalues.
The effective TIO potential of Eq. (20) reads

$\begin{split}V_{2d}(\psi_{1},\psi_{2})&=\mathcal{Q}\frac{\mathcal{P}+\mathcal{G}\mathcal{P}^{-1}}{2}|\psi_{1}|^{4}+\mathcal{Q}\frac{\mathcal{P}^{-1}+\mathcal{G}\mathcal{P}}{2}|\psi_{2}|^{4}\\\
&+\mathcal{Q}(\mathcal{G}-1)|\psi_{1}|^{2}|\psi_{2}|^{2}-\frac{2}{3}\left(\mathcal{P}|\psi_{1}|^{2}+\frac{|\psi_{2}|^{2}}{\mathcal{P}}\right)^{3/2}\\\
&-\mu(|\psi_{1}|^{2}+|\psi_{2}|^{2}).\end{split}$ (21)

Below, we consider the Langevin equations corresponding to Eqs. (18), namely

$\begin{split}\frac{\partial\psi_{1}}{\partial
t}&=-\Big{[}-\frac{1}{2}\frac{\partial^{2}\psi_{1}}{\partial
z^{2}}+\mathcal{Q}(\mathcal{P}+\mathcal{G}\mathcal{P}^{-1})|\psi_{1}|^{2}\psi_{1}\\\
&-\mathcal{Q}(1-\mathcal{G})|\psi_{2}|^{2}\psi_{1}-\mathcal{P}\sqrt{\mathcal{P}|\psi_{1}|^{2}+\mathcal{P}^{-1}|\psi_{2}|^{2}}\psi_{1}\\\
&-\mu_{1}\psi_{1}\Big{]}+\xi_{1}(z,t),\\\ \frac{\partial\psi_{2}}{\partial
t}&=-\Big{[}-\frac{1}{2}\frac{\partial^{2}\psi_{2}}{\partial
z^{2}}+\mathcal{Q}(\mathcal{P}^{-1}+\mathcal{G}\mathcal{P})|\psi_{2}|^{2}\psi_{2}\\\
&-\mathcal{Q}(1-\mathcal{G})|\psi_{1}|^{2}\psi_{2}-\frac{1}{\mathcal{P}}\sqrt{\mathcal{P}^{-1}|\psi_{2}|^{2}+\mathcal{P}|\psi_{1}|^{2}}\psi_{2}\\\
&-\mu_{2}\psi_{2}\Big{]}+\xi_{2}(z,t).\end{split}$ (22)

For the numerical simulation of Eq. (22) we utilize 1000 trajectories with the
domain size $L=8\times 8$. We first set
$\mu_{1}\approx\mu_{2}\equiv\mu_{0}+0.00001$ and $\mathcal{P}=1$. Further, by
considering two different noise distributions ($\xi_{1}$ and $\xi_{2}$) for
the two individual components, we break the symmetry between them. Figure 8
depicts the results of the Langevin dynamics (left panels) and the TIO
solutions (right panels). Comparing the corresponding 2D probability
distribution functions
$P[{\textrm{R}e}(\psi_{1})=u_{1},\textrm{Re}(\psi_{2})=u_{2}]$ between the
Langevin dynamics and the TIO findings for various inverse temperatures we
deduce that they are in a good agreement besides the region of small
temperatures. Particularly, as the temperature increases the distribution
widens similarly to the single-component case. However, at $\beta=95$, we
observe a mismatch among the predicted distributions. Indeed, the center of
the distribution within the TIO case is shifted from the origin $(0,0)$ in
contrast to the Langevin scenario where it is strongly localized around the
center. Furthermore, in order to estimate the impact of an increasing
population asymmetry we change $\mathcal{P}=1$ to $\mathcal{P}=1.25$. The
corresponding PDAs for $\beta=2$ are provided in Fig. 9 both within the
Langevin method and the TIO solution. Evidently, the predictions of these
approaches are quite proximal.

Figure 10: Time-evolution of $|\psi(z,t)|$ in the presence of temperature
showcasing the spontaneous creation of droplets due to the MI dynamics and
their consequent coalescence. We use $\mu=\mu_{0}+0.00001$. The
phenomenological factor accounting for the presence of temperature corresponds
here to $\gamma=0.03$.

##  Appendix C Phenomenological consideration of finite temperature effects
towards droplet formation

Next, we briefly discuss the role of finite (condensate) temperature effects
in the spontaneous generation of the droplets, relying on certain
approximations. It is important to recall that so far, we have been referring
to effective temperatures and chemical potentials, associated with the (free)
energy and atom number of the system in the grand-canonical ensemble
formulation herein. To emulate the role of thermal processes Proukakis and
Jackson (2008), we introduce a phenomenological damping factor $\gamma$ in the
MGP equation by replacing $i\rightarrow(i-\gamma)$ in Eq. (1). Recall that
this is a widely used approximation in the context of the common Gross-
Pitaevskii framework, i.e. in the absence of Lee-Huang-Yang contributions
Kasamatsu  _et al._ (2003); Choi  _et al._ (1998); Mithun  _et al._ (2016);
Katsimiga  _et al._ (2021). This factor accounts for the inclusion of thermal
components in the condensate being usually present in cold atom experiments
but it ignores any possible temperature dependence stemming from the Lee-
Huang-Yang contribution. The temperature dependence of the latter term and its
competition with the thermal condensate fraction are still open problems and
certainly worth to be pursued in future studies. Since the total number of
atoms is not conserved in the presence of $\gamma$, we treat the chemical
potential as a time-dependent parameter and adjust it at each evolution time
Mithun  _et al._ (2018).

Figure 10 shows the time-evolution of $|\psi(z,t)|$ in the case of
$\gamma=0.03$. Notice that in line with the discussion of Ref. Yan  _et al._
(2014), values of $\gamma$ within $0.00023--0.0023$ are relevant for
temperatures in the order of $10--100$nK. Nevertheless, we utilize a larger
value here to render apparent the fact that the nucleation of droplets takes
place within the same timescale as compared to the $\gamma=0$ scenario [Fig.
6(a)] and also their original number is the same. However, it should also be
noted that the presence of temperature enhances the formation of the ”large”
droplets by minimizing the background fluctuations, compare in particular Fig.
10 and Fig. 6(a). As mentioned above, a further systematic exploration of
finite temperature condensates, including the role of such corrections to LHY
terms is certainly a topic of particular interest for further study.

## References

  * Pethick and Smith (2008) C. J. Pethick and H. Smith,  _Bose–Einstein condensation in dilute gases_ (Cambridge university press, Cambridge, 2008).
  * Stringari and Pitaevskii (2003) S. Stringari and L. Pitaevskii,  _Bose-Einstein Condensation_ (Oxford University Press, Oxford, United Kingdom, 2003).
  * Kevrekidis  _et al._ (2015) P. G. Kevrekidis, D. J. Frantzeskakis, and R. Carretero-González,  _The defocusing nonlinear Schrödinger equation: from dark solitons to vortices and vortex rings_ (Society for Industrial and Applied Mathematics, Philadelphia, PA, 2015).
  * Petrov (2015) D. S. Petrov, “Quantum mechanical stabilization of a collapsing bose-bose mixture,” Phys. Rev. Lett. 115, 155302 (2015).
  * Petrov and Astrakharchik (2016) D. S. Petrov and G. E. Astrakharchik, “Ultradilute low-dimensional liquids,” Phys. Rev. Lett. 117, 100401 (2016).
  * Lee  _et al._ (1957) T. D. Lee, K. Huang, and C. N. Yang, “Eigenvalues and eigenfunctions of a bose system of hard spheres and its low-temperature properties,” Phys. Rev. 106, 1135 (1957).
  * Cabrera  _et al._ (2018) C. R. Cabrera, L. Tanzi, J. Sanz, B. Naylor, P. Thomas, P. Cheiney, and L. Tarruell, “Quantum liquid droplets in a mixture of bose-einstein condensates,” Science 359, 301 (2018).
  * Cheiney  _et al._ (2018) P. Cheiney, C. R. Cabrera, J. Sanz, B. Naylor, L. Tanzi, and L. Tarruell, “Bright soliton to quantum droplet transition in a mixture of bose-einstein condensates,” Phys. Rev. Lett. 120, 135301 (2018).
  * Semeghini  _et al._ (2018) G. Semeghini, G. Ferioli, L. Masi, C. Mazzinghi, L. Wolswijk, F. Minardi, M. Modugno, G. Modugno, M. Inguscio, and M. Fattori, “Self-bound quantum droplets of atomic mixtures in free space,” Phys. Rev. Lett. 120, 235301 (2018).
  * Ferioli  _et al._ (2019) G. Ferioli, G. Semeghini, L. Masi, G. Giusti, G. Modugno, M. Inguscio, A. Gallemí, A. Recati, and M. Fattori, “Collisions of self-bound quantum droplets,” Phys. Rev. Lett. 122, 090401 (2019).
  * D’Errico  _et al._ (2019) C. D’Errico, A. Burchianti, M. Prevedelli, L. Salasnich, F. Ancilotto, M. Modugno, F. Minardi, and C. Fort, “Observation of quantum droplets in a heteronuclear bosonic mixture,” Phys. Rev. Research 1, 033155 (2019).
  * Chomaz  _et al._ (2016) L. Chomaz, S. Baier, D. Petter, M. J. Mark, F. Wächtler, L. Santos, and F. Ferlaino, “Quantum-fluctuation-driven crossover from a dilute bose-einstein condensate to a macrodroplet in a dipolar quantum fluid,” Phys. Rev. X 6, 041039 (2016).
  * Ferrier-Barbut  _et al._ (2016) I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel, and T. Pfau, “Observation of quantum droplets in a strongly dipolar bose gas,” Phys. Rev. Lett. 116, 215301 (2016).
  * Ferioli  _et al._ (2020) G. Ferioli, G. Semeghini, S. Terradas-Briansó, L. Masi, M. Fattori, and M. Modugno, “Dynamical formation of quantum droplets in a ${}^{39}\mathrm{K}$ mixture,” Phys. Rev. Research 2, 013269 (2020).
  * Fort and Modugno (2021) C. Fort and M. Modugno, “Self-evaporation dynamics of quantum droplets in a 41k-87rb mixture,” Applied Sciences 11 (2021), 10.3390/app11020866.
  * Li  _et al._ (2018) Y. Li, Z. Chen, Z. Luo, C. Huang, H. Tan, W. Pang, and B. A. Malomed, “Two-dimensional vortex quantum droplets,” Phys. Rev. A 98, 063602 (2018).
  * Zhang  _et al._ (2019) X. Zhang, X. Xu, Y. Zheng, Z. Chen, B. Liu, C. Huang, B. A. Malomed, and Y. Li, “Semidiscrete quantum droplets and vortices,” Phys. Rev. Lett. 123, 133901 (2019).
  * Kartashov  _et al._ (2018) Y. V. Kartashov, B. A. Malomed, L. Tarruell, and L. Torner, “Three-dimensional droplets of swirling superfluids,” Phys. Rev. A 98, 013612 (2018).
  * Luo  _et al._ (2021) Z. Luo, W. Pang, B. Liu, Y. Li, and B. A. Malomed, “A new form of liquid matter: quantum droplets,” Frontiers of Physics 16, 32201 (2021).
  * Astrakharchik and Malomed (2018) G. E. Astrakharchik and B. A. Malomed, “Dynamics of one-dimensional quantum droplets,” Phys. Rev. A 98, 013631 (2018).
  * Mithun  _et al._ (2020) T. Mithun, A. Maluckov, K. Kasamatsu, B. A. Malomed, and A. Khare, “Modulational instability, inter-component asymmetry, and formation of quantum droplets in one-dimensional binary bose gases,” Symmetry 12, 174 (2020).
  * Morera  _et al._ (2020) I. Morera, G. E. Astrakharchik, A. Polls, and B. Juliá-Díaz, “Quantum droplets of bosonic mixtures in a one-dimensional optical lattice,” Phys. Rev. Research 2, 022008 (2020).
  * Griffin  _et al._ (2009) A. Griffin, T. Nikuni, and E. Zaremba,  _Bose-condensed gases at finite temperatures_ (Cambridge University Press, Cambridge, 2009).
  * De Rosi  _et al._ (2021) G. De Rosi, G. E. Astrakharchik, and P. Massignan, “Thermal instability, evaporation, and thermodynamics of one-dimensional liquids in weakly interacting bose-bose mixtures,” Phys. Rev. A 103, 043316 (2021).
  * Parisi  _et al._ (2019) L. Parisi, G. E. Astrakharchik, and S. Giorgini, “Liquid state of one-dimensional bose mixtures: A quantum monte carlo study,” Phys. Rev. Lett. 122, 105302 (2019).
  * Parisi and Giorgini (2020) L. Parisi and S. Giorgini, “Quantum droplets in one-dimensional bose mixtures: A quantum monte carlo study,” Phys. Rev. A 102, 023318 (2020).
  * Mistakidis  _et al._ (2021) S. I. Mistakidis, T. Mithun, P. G. Kevrekidis, H. R. Sadeghpour, and P. Schmelcher, “Formation and quench of homo-and hetero-nuclear quantum droplets in one-dimension,” arXiv:2108.00727 (2021).
  * Astrakharchik and Giorgini (2006) G. E. Astrakharchik and S. Giorgini, “Correlation functions of a lieb–liniger bose gas,” J. Phys. B: At. Mol. and Opt. Phys. 39, S1 (2006).
  * Blakie  _et al._ (2008) P. B. Blakie, A. Bradley, M. Davis, R. Ballagh, and C. Gardiner, “Dynamics and statistical mechanics of ultra-cold bose gases using c-field techniques,” Advances in Physics 57, 363 (2008).
  * Sreedharan  _et al._ (2020) A. Sreedharan, S. Choudhury, R. Mukherjee, A. Streltsov, and S. Wüster, “Collisions of solitary waves in condensates beyond mean-field theory,” Phys. Rev. A 101, 043604 (2020).
  * Deuar and Drummond (2006a) P. Deuar and P. Drummond, “First-principles quantum dynamics in interacting bose gases: I. the positive p representation,” J. Phys. A: Mathematical and General 39, 1163 (2006a).
  * Deuar and Drummond (2006b) P. Deuar and P. Drummond, “First-principles quantum dynamics in interacting bose gases ii: stochastic gauges,” J. Phys. A: Mathematical and General 39, 2723 (2006b).
  * Scalapino  _et al._ (1972) D. J. Scalapino, M. Sears, and R. A. Ferrell, “Statistical mechanics of one-dimensional ginzburg-landau fields,” Phys. Rev. B 6, 3409 (1972).
  * Khare  _et al._ (1997) A. Khare, S. Habib, and A. Saxena, “Exact thermodynamics of the double sinh-gordon theory in 1+ 1 dimensions,” Phys. Rev. Lett. 79, 3797 (1997).
  * Nunnenkamp  _et al._ (2007) A. Nunnenkamp, J. N. Milstein, and K. Burnett, “Classical field techniques for condensates in one-dimensional rings at finite temperatures,” Phys. Rev. A 75, 033604 (2007).
  * Habib and Lythe (2000) S. Habib and G. Lythe, “Dynamics of kinks: Nucleation, diffusion, and annihilation,” Phys. Rev. Lett. 84, 1070 (2000).
  * Aarts  _et al._ (2000) G. Aarts, G. F. Bonini, and C. Wetterich, “On thermalization in classical scalar field theory,” Nuclear Physics B 587, 403 (2000).
  * Tylutki  _et al._ (2020) M. Tylutki, G. E. Astrakharchik, B. A. Malomed, and D. S. Petrov, “Collective excitations of a one-dimensional quantum droplet,” Phys. Rev. A 101, 051601(R) (2020).
  * Lavoine and Bourdel (2021) L. Lavoine and T. Bourdel, “Beyond-mean-field crossover from one dimension to three dimensions in quantum droplets of binary mixtures,” Phys. Rev. A 103, 033312 (2021).
  * Pippard (1964) A. B. Pippard,  _Elements of classical thermodynamics: for advanced students of physics_ (Cambridge University Press, 1964).
  * Nguyen  _et al._ (2017) J. H. Nguyen, D. Luo, and R. G. Hulet, “Formation of matter-wave soliton trains by modulational instability,” Science 356, 422 (2017).
  * Gelash  _et al._ (2019) A. Gelash, D. Agafontsev, V. Zakharov, G. El, S. Randoux, and P. Suret, “Bound state soliton gas dynamics underlying the spontaneous modulational instability,” Phys. Rev. Lett. 123, 234102 (2019).
  * Kraych  _et al._ (2019) A. E. Kraych, D. Agafontsev, S. Randoux, and P. Suret, “Statistical properties of the nonlinear stage of modulation instability in fiber optics,” Phys. Rev. Lett. 123, 093902 (2019).
  * Trullinger and DeLeonardis (1980) S. E. Trullinger and R. M. DeLeonardis, “Analytic statistical mechanics for a two-component-kink system,” Phys. Rev. B 22, 5522 (1980).
  * Naraschewski and Glauber (1999) M. Naraschewski and R. J. Glauber, “Spatial coherence and density correlations of trapped bose gases,” Phys. Rev. A 59, 4595 (1999).
  * Katsimiga  _et al._ (2017) G. C. Katsimiga, G. M. Koutentakis, S. I. Mistakidis, P. G. Kevrekidis, and P. Schmelcher, “Dark–bright soliton dynamics beyond the mean-field approximation,” New J. Phys. 19, 073004 (2017).
  * Mistakidis  _et al._ (2018) S. I. Mistakidis, G. C. Katsimiga, P. G. Kevrekidis, and P. Schmelcher, “Correlation effects in the quench-induced phase separation dynamics of a two species ultracold quantum gas,” New J. Phys. 20, 043052 (2018).
  * Nandkishore and Huse (2015) R. Nandkishore and D. A. Huse, “Many-body localization and thermalization in quantum statistical mechanics,” Annu. Rev. Condens. Matter Phys. 6, 15 (2015).
  * Parisi  _et al._ (1981) G. Parisi, Y. S. Wu, _et al._ , “Perturbation theory without gauge fixing,” Scientia Sinica 24, 483 (1981).
  * Hohenberg and Halperin (1977) P. C. Hohenberg and B. I. Halperin, “Theory of dynamic critical phenomena,” Rev. Mod. Phys. 49, 435 (1977).
  * Dennis  _et al._ (2013) G. R. Dennis, J. J. Hope, and M. T. Johnsson, “Xmds2: Fast, scalable simulation of coupled stochastic partial differential equations,” Computer Physics Communications 184, 201 (2013).
  * Everitt  _et al._ (2017) P. J. Everitt, M. A. Sooriyabandara, M. Guasoni, P. B. Wigley, C. H. Wei, G. D. McDonald, K. S. Hardman, P. Manju, J. D. Close, C. C. N. Kuhn, S. S. Szigeti, Y. S. Kivshar, and N. P. Robins, “Observation of a modulational instability in bose-einstein condensates,” Phys. Rev. A 96, 041601(R) (2017).
  * Mistakidis  _et al._ (2015) S. I. Mistakidis, T. Wulf, A. Negretti, and P. Schmelcher, “Resonant quantum dynamics of few ultracold bosons in periodically driven finite lattices,” J. Phys. B: At. Mol. and Opt. Phys. 48, 244004 (2015).
  * Bloch  _et al._ (2008) I. Bloch, J. Dalibard, and W. Zwerger, “Many-body physics with ultracold gases,” Rev. Mod. Phys. 80, 885 (2008).
  * Proukakis and Jackson (2008) N. P. Proukakis and B. Jackson, “Finite-temperature models of bose–einstein condensation,” J. Phys. B: At. Mol. and Opt. Phys. 41, 203002 (2008).
  * Kasamatsu  _et al._ (2003) K. Kasamatsu, M. Tsubota, and M. Ueda, “Nonlinear dynamics of vortex lattice formation in a rotating bose-einstein condensate,” Phys. Rev. A 67, 033610 (2003).
  * Choi  _et al._ (1998) S. Choi, S. A. Morgan, and K. Burnett, “Phenomenological damping in trapped atomic bose-einstein condensates,” Phys. Rev. A 57, 4057 (1998).
  * Mithun  _et al._ (2016) T. Mithun, K. Porsezian, and B. Dey, “Disorder-induced vortex lattice melting in a bose-einstein condensate,” Phys. Rev. A 93, 013620 (2016).
  * Katsimiga  _et al._ (2021) G. C. Katsimiga, S. I. Mistakidis, P. Schmelcher, and P. G. Kevrekidis, “Phase diagram, stability and magnetic properties of nonlinear excitations in spinor bose–einstein condensates,” New J. Phys. 23, 013015 (2021).
  * Mithun  _et al._ (2018) T. Mithun, S. C. Ganguli, P. Raychaudhuri, and B. Dey, “Signatures of two-step impurity-mediated vortex lattice melting in bose-einstein condensates,” Europhys. Lett.) 123, 20004 (2018).
  * Yan  _et al._ (2014) D. Yan, R. Carretero-González, D. J. Frantzeskakis, P. G. Kevrekidis, N. P. Proukakis, and D. Spirn, “Exploring vortex dynamics in the presence of dissipation: Analytical and numerical results,” Phys. Rev. A 89, 043613 (2014).