# Decode Neural signal as Speech

Yiqian Yang  Yiqun Duan  Qiang Zhang  Renjing Xu  Hui Xiong

###### Abstract

Decoding language from brain dynamics is an important open direction in the
realm of brain-computer interface (BCI), especially considering the rapid
growth of large language models. Compared to invasive-based signals which
require electrode implantation surgery, non-invasive neural signals (e.g. EEG,
MEG) have attracted increasing attention considering their safety and
generality. However, the exploration is not adequate in three aspects: 1)
previous methods mainly focus on EEG but none of the previous works address
this problem on MEG with better signal quality; 2) prior works have
predominantly used “teacher-forcing” during generative decoding, which is
impractical; 3) prior works are mostly “BART-based” not fully auto-regressive,
which performs better in other sequence tasks. In this paper, we explore the
brain-to-text translation of MEG signals in a speech-decoding formation. Here
we are the first to investigate a cross-attention-based “whisper” model for
generating text directly from MEG signals without teacher forcing. Our model
achieves impressive BLEU-1 scores of 60.30 and 52.89 without pretraining &
teacher-forcing on two major datasets (GWilliams and Schoffelen). This paper
conducts a comprehensive review to understand how speech decoding formation
performs on the neural decoding tasks, including pretraining initialization,
training & evaluation set splitting, augmentation, and scaling law. Code is
available at https://github.com/NeuSpeech/NeuSpeech1/tree/main.

##  1 Introduction

Decoding language from brain activity is an important and rapidly advancing
area of neurotechnology. The ability to infer a person’s intended words,
phrases, or even full sentences directly from their brain signals has profound
applications for improving communication and control capabilities for those
suffering from severe speech and motor impairments. Conditions like high-level
spinal cord injury or late-stage ALS can trap a person inside their own body,
unable to fully express themselves or interact with the digital world around
them. By noninvasively measuring and modeling the complex neural
representations underlying language production in the brain, we aim to restore
digital communication simply through thought. Beyond assisting those with
disabilities, proficient decoding of full linguistic expressions from brain
signals could also enable entirely new modes of seamless human-machine
interaction. One could silently yet precisely command prosthetic devices,
control software applications, browse the internet, or navigate virtual and
augmented reality environments entirely using their brain. The potential
applications seem nearly limitless and could vastly expand what is possible
for both able-bodied individuals and those with severe impairments. It is an
enormously challenging task from an artificial intelligence and neuroscience
perspective but one with immense promise to transform lives if achieved.

Notably, using decoding speech on invasive signals achieved high-performance
recognition accuracy for real-time translation (Metzger et al., 2023; Willett
et al., 2023). However, the invasive recordings pose medical risks and issues
with maintaining over long periods. Thus, for non-invasive brain recordings,
Wang et al. (Wang & Ji, 2022) demonstrated open-vocabulary EEG-to-text
translation by employing pre-trained language models on word-level EEG
features. Duan et al. (Duan et al., 2023) improved this approach by decoding
from raw EEG waves without using time markers. However, limitations remain as
evaluation relies on teacher forcing and the model cannot generate meaningful
sentences without it. On MEG, previous works can merely decode limited classes
from MEG signal. (Dash et al., 2020; Csaky et al., 2023; Ghazaryan et al.,
2023). Recently, Defossez et al.(Défossez et al., 2023) decoded speech
perception of 3 second segments from MEG signals, and was capable of
classifying sentences by matching MEG signal with speech. However, it cannot
generate sentences from the MEG signal.

To address the issue of teacher forcing evaluation and generating text from
MEG signal. In this paper, we present NeuSpeech, a pioneering framework that
decodes neural as speech. NeuSpeech uses an encoder-decoder model (whisper
(Radford et al., 2023)) rather than BART-based models that are widely used in
previous works (Wang & Ji, 2022; Duan et al., 2023; Xi et al., 2023) to
directly translate raw MEG waves into text, without pretraining or
transforming neural signals to discrete codes. This approach 1) lets the model
efficiently learn text-related information from neural signals, and 2) keeps
the general meaning of the whole sentence when generating long texts.

Experiments employ non-invasive MEG signals from the large-scale public
GWilliams (Gwilliams et al., 2023) and Schoffelen (Schoffelen et al., 2019)
datasets, which recorded MEG during speech listening tasks. Notably, NeuSpeech
can generalize to all languages and equipment types. Performance is assessed
using standard translation metrics (Papineni et al., 2002; Lin, 2004). For raw
MEG waves without event markers, NeuSpeech achieves 60.30 BLEU-1 and 55.26
Rouge-1 F-measure on GWilliams, and 53.16 BLEU-1 on Schoffelen without teacher
forcing evaluation. The main contributions of this work are as follows:

  * •

NeuSpeech adoptes an end-to-end neural framework for open-vocabulary MEG-to-
text translation without markers without pre-training or teacher forcing, and
is the first to achieve high-performance raw MEG wave-to-text translation.

  * •

We explore multi-layout, multi-language, multi-subject MEG-to-text
transcription on multiple datasets under various training settings. This
provides a more comprehensive evaluation than prior work.

  * •

We exploit cross-attention to get global semantic information for sentences
rather than previous token-based methods.

##  2 Related Works

Significant progress has been made in the field of brain-to-text decoding in
recent years. In 2019, Anumanchipalli et al. (Anumanchipalli et al., 2019)
proposed a model that could decode ECoG signals into the kinematics of
articulation, followed by acoustic features (such as MFCCs), and ultimately
generate the speech waveform, which inspired many researchers to explore this
area further. In 2020, Wang et al. (Wang et al., 2020) introduced the use of
generative adversarial networks (GANs) to decode ECoG signals into speech. In
2021, Willett et al. (Willett et al., 2021) proposed a model that utilized an
RNN and a probabilistic language model to decode letters while a subject was
performing handwriting. Metzger et al. (Metzger et al., 2022) developed a
pipeline that decoded ECoG signals into text using an RNN and then refined the
results using the GPT-2 model. In the domain of open-vocabulary decoding,
Metzger et al. (Metzger et al., 2023) proposed an RNN framework that could
decode speech, text, emotion, and facial expressions from ECoG signals in
real-time. Simultaneously, Willett et al. (Willett et al., 2023) successfully
decoded text from brain signals. Liu et al. (Liu et al., 2023) developed a
three-part model to decode logosyllabic languages (e.g., Chinese) by first
decoding ECoG signals into Chinese pinyin with tones and syllables, and then
synthesizing speech. Feng et al. (Feng et al., 2023) also achieved text
decoding from SEEG signals. However, it is important to note that these
practical systems predominantly rely on invasive brain recordings.

For non-invasive recordings, Meta proposed a brain-to-speech framework using
contrastive learning with MEG and EEG signals (Défossez et al., 2023), which
is capable of classifying limited sentences but cannot be used for open-
vocabulary text decoding. Ghazaryan et al. (Ghazaryan et al., 2023) attempted
to decode limited words using MEG responses. Wang et al. (Wang & Ji, 2022)
developed a framework to decode word-level EEG signal features into text using
a pre-trained BART model (Lewis et al., 2020). Many subsequent works have
followed in the footsteps of Wang et al. (Wang & Ji, 2022), including Dewave
(Duan et al., 2023), which proposed a framework that employed wave2vec
(Baevski et al., 2020), and discrete codex to obtain effective representations
and fed them into a BART (Lewis et al., 2020) model to generate text. However,
these methods rely on teacher forcing, which is not feasible in real-world
settings. Different from these methods, our work uses an end-to-end paradigm
to generate open-vocabulary text and evaluate it without teacher forcing.

Table 1:  Evaluation metrics of MEG-to-Text translation under raw waves input.
For a fair comparison, these results has w/ and w/o the teacher-forcing
evaluation setting since EEG-to-Text (Wang & Ji, 2022) only has teacher-
forcing originally. Default is w/o teacher forcing and ‘w/ tf’ means with
teacher-forcing.

| BLEU-N (%) $\uparrow$ | ROUGE-1 (%)$\uparrow$ | WER (%) $\uparrow$
---|---|---|---
Modality | Method | N=1 | N=2 | N=3 | N=4 | F | P | R |
Speech | Whisper (Radford et al., 2023) | 67.79 | 59.15 | 51.97 | 45.78 | 75.90 | 73.53 | 81.24 | 41.21
Noise | NeuSpeech | 0.0 | 0.0 | 0.0 | 0.0 | 2.84 | 23.9 | 1.5 | 99.98
- | random | 4.99 | 1.18 | 0.59 | 0.41 | 4.79 | 5.95 | 5.87 | 132.30
MEG |  eeg-to-text (Wang & Ji, 2022) | 9.21 | 2.13 | 0.57 | 0.14 | 9.74 | 10.73 | 11.38 | 118.25
MEG | NeuSpeech | 57.68 | 52.76 | 48.84 | 45.42 | 60.24 | 61.86 | 61.92 | 58.89
MEG |  NeuSpeech w/ pp | 60.3 | 55.26 | 51.24 | 47.78 | 58.73 | 60.88 | 59.76 | 56.63
MEG |  eeg-to-text (Wang & Ji, 2022) w/ tf | 21.66 | 9.27 | 3.93 | 1.73 | 20.11 | 22.48 | 18.71 | 95.69
MEG |  NeuSpeech w/ tf | 78.13 | 71.25 | 64.79 | 58.86 | 81.91 | 80.52 | 83.68 | 25.38
MEG |  NeuSpeech w/ pp tf | 80.20 | 73.27 | 66.79 | 60.81 | 78.62 | 78.14 | 79.54 | 23.42

##  3 Method

NeuSpeech is a modification of whisper (Radford et al., 2023) model, where the
raw MEG signals are sent into the model and generate text. Whisper (Radford et
al., 2023) is the most powerful and noise-robust model in speech recognition
which is trained on hundreds of thousands of audio data and can transcribe
speech from most of the languages on earth.

###  3.1 Task Definition

Given a raw MEG waves $\mathcal{X}$, the aim is to decode the corresponding
open-vocabulary text tokens $\mathcal{W}$. These MEG-Text pairs
$\langle\mathcal{X},\mathcal{W}\rangle$ are collected during hearing, as
defined in Section 4.1. We decode neural signals as speech, and we only have
one setting which is: Raw MEG Waves to Text Translation, where MEG waves
$\mathcal{X}$ are directly sent into the whisper model to generate continuous
words without teacher forcing, a more challenging but practical setting.
NeuSpeech is the pioneering work in this task. Different from previous
works(Duan et al., 2023; Wang & Ji, 2022), we are the first to evaluate
generated text without teacher forcing using MEG.

###  3.2 Preprocessing

We only implement basic pre-processing here, the MEG signals are notched at
power line frequency and filtered with (1,60) Hz. Then, both of them are
resampled to 200Hz. After that, the signals are normalized using robust scaler
mentioned in (Défossez et al., 2023), then clipped to $\pm 10$ and scaled into
a value range of $(-1,1)$.

###  3.3 Model

Different from previous works (Duan et al., 2023; Tang et al., 2022) which are
usually trained through multi-steps. Our model is trained end-to-end from
scratch without bells and whistles. We freeze parameters in the decoder and
only train the encoder. Since the original input channel of the whisper model
is incompatible with the MEG channels, we replace the first convolution layer
of the model with two convolution layers to fit in MEG channels and align time
width respectively. For training data, we do not use any data augmentation. In
terms of efficiency, we use Adalora (Zhang et al., 2023) technique to tune our
model, and only the encoder including convolution layers are in the Lora
model.

###  3.4 Evaluation

We evaluate transcribing performance on GWilliams dataset (Gwilliams et al.,
2023) and Schoffelen (Schoffelen et al., 2019) using NLP metrics, BLEU, ROUGE
and WER, which are used in previous works (Wang & Ji, 2022; Duan et al.,
2023). We add three baselines to give a rough understanding of up and down
limits. 1) speech recognition, we use a pre-trained whisper base model to
demonstrate the high performance on speech modality. 2) input noise, we input
Gaussian noise into our model to check if our model can distinguish between
notoriously noisy MEG and noise, and the well-trained model should not
generate text with pure noise as input. 3) random selection, we also showed
the performance of randomly selecting sentences from training sets as
predictions to give the down limit of the dataset as learned models should
pass this limit to prove its effectiveness.

##  4 Experiments

###  4.1 Dataset

We utilize GWilliams (Gwilliams et al., 2023) and Schoffelen (Schoffelen et
al., 2019) datasets which are auditory MEG task for experiments. Both have a
detailed timetable for each sentence, and they are divided by sentence-
recording pairs into train:validation:test 8:1:1, which results in
23339:2917:2918 and 8570:1112:1076 MEG-text pairs respectively. However, we
should note that the unique sentences are 661 and 360 for two datasets, which
are very limited compared to speech datasets.

Table 2: Transcription results on the unseen MEG waves, where bold denotes a
correct match between ground truth and our prediction. Underline denotes a
fuzzy match with similar semantic meanings.  For a fair comparison, these
results do not keep the teacher-forcing evaluation setting as EEG-to-Text
(Wang & Ji, 2022) and Dewave (Duan et al., 2023). This means the decoding
process can hold on accumulated errors and predicts the next token with model-
generated tokens from previous steps.

Decoding Results on GWilliams (Gwilliams et al., 2023)
---
Ground Truth: At best they tales at worst expeditions into macabre realms no
healthy mind need ever see
Prediction: At best it is macred at worst expeditions into macabre realms no
healthy mind need ever see that
Ground Truth: He turned to regard the pillar still rising behind his house and
ran to see what had been done
Prediction: He-turned to regard a pillar still rising behind his house and ran
away from…to see what had been done.
Decoding Results on Schoffelen (Schoffelen et al., 2019)
Ground Truth: Tijdens de rumoerige bespreking beslisten de leden de staking
voort te zetten
Prediction: Tijdens de rumoerige bespreking beslisten de leden of de staking
voort.
Ground Truth: Op die goede prestatie was helemaal niets af te dingen
Prediction: Op die goede prestatie was naar een varken afkomen van het nieuwe
afguk…

###  4.2 Implementation Details

Model modification. We adopt whisper base model (Radford et al., 2023) as our
base model. To adapt the input size of the neural signal, we use two
convolution1d layers, the kernel sizes are 3, and stride is 1,2 respectively
to replace the first convolution layer in the whisper-base model of dimension
512.

Training. To implement eeg-to-text (Wang & Ji, 2022), we use the original word
timestamps provided by GWilliams (Gwilliams et al., 2023) to segment MEG waves
for each word feature. The training settings are the same as the origin code.

All models are trained on Nvidia H800 GPUs. We use a learning rate of 1e-3 and
batch size of 64 for training 500 epochs and select the best model on
evaluation loss. We use the AdamW as the optimizer for training all the
models.

Evaluation. For non-teacher forcing generation configuration, we use the beam
search with a beam size of 5 and set the repetition penalty to 5, with no
repeat n-gram size to 2, which is very important to avoid repeating. In
addition, we implement post-processing composed of filtering and lowering
alphabets to improve performance.

###  4.3 Metrics

As shown in Table 1, we compare our results to input Gaussian noise to the
NeuSpeech model and randomly select sentences from the test set and whisper
base model speech transcription to give the understanding of up and down
limits. We also compare previous method eeg-to-text (Wang & Ji, 2022).

Based on the information provided in Tab. 1, it is evident that our method
outperforms randomly selected sentences, noise input, and eeg-to-text (Wang &
Ji, 2022). NeuSpeech has demonstrated the feasibility of decoding text from
non-invasive brain recordings. Prior to NeuSpeech, many previous works (Duan
et al., 2023; Wang & Ji, 2022) in this field relied on teacher forcing to
evaluate their models, Defossez et al. (Défossez et al., 2023) used
classification accuracy. However, NeuSpeech introduces a more fair and real-
world setting by employing end-to-end training, which directly ties neural
signals to text. In comparison to previous works evaluated using teacher
forcing, NeuSpeech has shown remarkable improvement, achieving BLEU-1, BLEU-2,
BLEU-3, and BLEU-4 scores of 60.3, 55.26, 51.24, and 47.78, respectively,
which is close to the performance of speech recognition.

To ensure a fair comparison, we also implement teacher forcing for both the
eeg-to-text model (Wang & Ji, 2022) and NeuSpeech. It is evident that our
BLEU-1 score of 80.20 is almost four times that of the eeg-to-text (Wang & Ji,
2022) model, highlighting the superiority of our approach. Furthermore, we
have surpassed speech recognition methods with teacher forcing, emphasizing
that our model can significantly improve performance by eliminating
accumulated errors. It is worth noting that, at this stage, BLEU is the most
suitable metric for evaluating the MEG-to-text system, as ROUGE and WER tend
to yield high scores for noise and random models.

###  4.4 Generated Samples

Table 3:  Transcription results comparison with different methods. Bold font
means exact match of words. Underline for similar but not exact matching. w/
means with. tf means teacher-forcing, pt means pretraining, jt means joint-
training. eeg-to-text (Wang & Ji, 2022) here is using beam search as the same
setting of ours which is described in Sec. 4.2.

(1) | GT |  Arthur replaced the tray onto the lacquered cherry table that stood in the center of his ordered chamber
---|---|---
NeuSpeech |  The rental car records pointed to lacquered cherry-schedung and stood in the center of his ordered chamber in motel from Acres comb in?
NeuSpeech w/ tf |  Arthur replaced thethe cartray onto thethe lacquered cherry-table inthat stood in thethe center of his ordered chamber in
NeuSpeech w/ pt |  Arthur replaced the tray onto the lacquered cherry table to match from Acres.
NeuSpeech w/ jt |  The rental car records pointed to Acres and stay in the center of his ordered chamber
eeg-to-text (Wang & Ji, 2022) |  He said he had no idea how long it would take him to get there
eeg-to-text (Wang & Ji, 2022) w/ tf |  was the old with the tablekeyered table table and had in the middle of the room table
(2) | GT |  My patience was long gone and I was back in the car to warming up when Acres tapped on the window and told me he had found whatever he was looking for
NeuSpeech |  My patience was long gone and I was back in the car. But when I heard that many of you were looking for whatever it was, but what about this?
NeuSpeech w/ tf |  My patience was long gone and I was back in athe car.to warming up when Acres wastapped on thethe window and openedtold me he had found whatever he was looking for.
NeuSpeech w/ pt |  My patience was long gone and I was back in the car.
NeuSpeech w/ jt |  My patience was long gone and I was back in the car. The window was looking good when Acres is ongoing on her car and she now has…
eeg-to-text (Wang & Ji, 2022) |  He said he had no idea how long it would take him to get back home
eeg-to-text (Wang & Ji, 2022) w/ tf |  father was running gone I was ready to the dark drive up Iacia me the window said me to was a a he was looking for
(3) | GT |  Very well said Arthur then take your story to that and let him defend you when no civilized publisher will approach you
NeuSpeech |  Very well said Arthur! Tell me your story, and I’ll defend it from you when no civilized publisher will approach you.
NeuSpeech w/ tf |  Very well said Arthur!then thattake your story to andthat and let him defend you and no civilized publisher will approach you when
NeuSpeech w/ pt |  Very well said Arthur then take your story to that and let him defend you when no civilized publisher will approach you very well
NeuSpeech w/ jt |  Very well said! It was that and we can make sure that it and not let him defend you when no civilized publisher will approach you
eeg-to-text (Wang & Ji, 2022) |  It was the first time I had seen him in a long time and it made me very nervous
eeg-to-text (Wang & Ji, 2022) w/ tf |  little then he the seat and the place go it take it he one man would you

In Tab. 2, we present visual examples of text that have been generated from
MEG signals. The training of Schoffelen (Schoffelen et al., 2019) is described
in sec. 4.5.1 These examples showcase the capabilities of our model in
generating meaningful text, despite the difficulties associated with
translating thoughts and the limited prior research in this area. While our
model may not yet achieve the same level of performance as nowadays speech
recognition tasks, it demonstrates the ability to align keywords and construct
sentences with similar structures. This is a promising indication of the
potential of our approach in the field of brain-to-text translation.

The model’s capability to match long sequences of words continuously is indeed
impressive. As evidenced by the examples provided in Tab. 2, such as the first
sentence with 12 identical words, the model demonstrates the ability to
accurately predict groups of words throughout a sentence. It is noteworthy
that the model can leverage information from the MEG signals to make
predictions, allowing it to generate correct groups of words even when there
are occasional errors in the middle of a sentence. This highlights the model’s
capacity to capture and utilize contextual information from the brain signals,
enabling it to generate meaningful and coherent text. Also, the model is good
at transcribing Dutch which also shows the capability of generating continuous
words.

We have provided Tab. 3 below, which compares the results of different
approaches, namely NeuSpeech, NeuSpeech with teacher forcing, NeuSpeech with
pre-training, and joint training. We observe that using NeuSpeech alone often
result in semantically coherent sentences, even if not all the words were
correct. Particularly in example (3), there was even a hint of expressive
tone. However, when NeuSpeech was generating text with teacher forcing,
although there were more correct word matches, there were often consecutive
word repetitions, such as ”thethe.” The reason for this occurrence is not yet
clear, but it significantly affects the semantic coherence of the sentences.

The performance of NeuSpeech improves to some extent when pre-training is
employed. This improvement may be attributed, in part, to the model becoming
more cautious and refraining from generating further output after a few
errors. Among all the examples, only the sentences generated by the pre-
trained model had word sequences without any incorrect words. This has both
advantages and disadvantages. The advantage is that it improves the precision
of the model and reduces the occurrence of incorrect words. However, the
disadvantage is that once an error is predicted, the model becomes
ineffective.

In contrast, the joint training model seems to tolerate longer sequences of
word errors. For example, in (1), the model does not predict a single correct
word for most of the sentence but accurately predicts the last seven words,
which is neglected by the higher-scoring pre-trained model. This indicates
that, even for a model that reaches a certain performance level, we cannot
solely rely on scoring metrics. Each model may have unique attributes that we
can exploit. In addition, previous approaches for EEG-to-text translation
seems to be less effective, as even with the use of teacher forcing, they
produced only a few words that match the true labels. This further highlights
the superiority of our approach.

In summary, the analysis of generated samples from different approaches
reveals valuable insights. NeuSpeech demonstrates the ability to generate
semantically coherent sentences, while NeuSpeech with teacher forcing and pre-
training offers improvements in word accuracy, albeit with some drawbacks. The
joint training model shows resilience to long sequences of word errors. These
findings emphasize the importance of considering both quantitative metrics and
understanding the unique properties and strengths of different models.
Furthermore, the comparison with previous EEG-to-text (Wang & Ji, 2022)
methods underscores the advancements achieved by our approach.

###  4.5 Ablation Study

To demonstrate the effectiveness of our model and provide experimental and
theoretical foundations for future improvements in model design, training
methods, data augmentation techniques, and modifications, we conducted an
ablation study. This study provided a comprehensive analysis of various
factors influencing the performance of our model. The experimental results and
theoretical insights obtained from these experiments not only validated the
effectiveness of our model but also laid the groundwork for designing improved
models, selecting appropriate training methods, devising better data
augmentation techniques, and making informed model modifications.

####  4.5.1 Different Dataset

We evaluated our model in Schoffelen dataset (Schoffelen et al., 2019) to
illustrate our effectiveness. As shown in Tab. 4, our model has comparable
performance on two datasets. This illustrates that our model is applicable
under many scenarios, including English and Dutch, and can be adapted on
different layouts.

Table 4: Different datasets. Bold for higher socres.

Dataset | BLEU-N (%) $\uparrow$ | ROUGE-1 (%) $\uparrow$ |
---|---|---|---
N=1 | N=2 | N=3 | N=4 | F | P | R |
GWilliams | 60.3 | 55.26 | 51.24 | 47.78 | 58.73 | 60.88 | 59.76 |
Schoffelen | 52.89 | 48.53 | 46.02 | 44.18 | 54.28 | 53.87 | 56.14 |

####  4.5.2 Splits

We conducted further experiments using the GWilliams (Gwilliams et al., 2023)
dataset to evaluate NeuSpeech’s generalizability. To test the model’s ability
to generalize across recording sessions, subjects, and unique sentences, we
developed three additional training/validation splits beyond the standard 8:1
split. Specifically, we withheld: 1) all data from the second recording
session, 2) all instances of the story “lw1”, and 3) a random 10% sample of
unique sentences. For the latter two splits, we ensured no overlap sentences
between the training and test sets for the held-out data.

Table 5: Dataset splits.Noise means input model with Gaussian noise. Random
means randomly choosing sentences from training set as predictions for testing
set. Retrieve means find best match in training set for each test sentence as
predictions using sentence similarity using “all-MiniLM-L6-v2” model.
Generation means our model’s generated text. Here, we only shows the BLEU-1
score. Note, for hold-out story and sentences, we used nucleus sampling p=0.25
for generation.

Hold-out | Noise | Random | Retrieve | Generation
---|---|---|---|---
Session | 0.0 | 4.94 | 100 | 53.16
Story | 0.0 | 4.97 | 11.01 | 6.91
Sentences | 0.0 | 4.79 | 16.38 | 6.54

Table 6: Transcription results on different splits, where bold denotes a
correct match between ground truth and our prediction. Underline denotes a
fuzzy match with similar semantic meanings. Hold-out sentences
---
Ground Truth: I looked down and saw he was pointing to another marker just
slightly different from the one we had just left
Prediction: It was different from now and it just wasn’t working well.
Hold-out story
Ground Truth: Mahogany black curly hair cut short in the typical military
style framed a dark skinned face
Prediction: I must admit, this is a great black or cold weather. It had things
such as

As shown in Tab. 5, the model maintained good performance when trained and
tested on separate recording sessions, demonstrating its ability to generalize
across sessions. NeuSpeech significantly outperformed random baselines on
splits with no training-test overlap, validating its generalizability despite
low text similarity between train and test sentences.Tab. 6 presents sample
outputs, revealing the model can still capture some semantic meaning and words
when generalizing, though it struggles to produce fully correct novel
sentences.

####  4.5.3 Pretrain

To investigate the effectiveness of pretraining, we conducted experiments
involving pretraining on one dataset and fine-tuning on the other.
Specifically, we pretrained the model on the Schoffelen dataset (Schoffelen et
al., 2019) and fine-tuned on GWilliams (Gwilliams et al., 2023), and vice
versa. As shown in Table 7, pretraining provided benefits in some settings.

Table 7: The performance of pretraining. The column Pretrain Dataset means
trained on GWilliams, Schoffelen and none. Bold for highest socres.

Pretrain Dataset |  Finetune Dataset | BLEU-N ($\%$) $\uparrow$
---|---|---
N=1 | N=2 | N=3 | N=4
GWilliams | Schoffelen | 48.09 | 42.23 | 38.5 | 35.60
Schoffelen | GWilliams | 63.68 | 61.06 | 59.26 | 57.76
- | Schoffelen | 52.89 | 48.53 | 46.02 | 44.18
- | GWilliams | 60.30 | 55.26 | 51.24 | 47.78

Interestingly, fine-tuning on GWilliams after pretraining achieved higher
performance than training from scratch, while pretraining on GWilliams and
fine-tuning on Schoffelen yielded much lower scores. One potential reason for
this asymmetry is that GWilliams contains over twice as many samples as
Schoffelen. Thus, after pretraining on the larger GWilliams set, the model is
already highly optimized, making it difficult for the training loss to further
descend during fine-tuning on the smaller Schoffelen data.

In contrast, fine-tuning the model pretrained on Schoffelen Still allowed for
improvements to be made on the larger GWilliams distributions. This suggests
pretraining can help in low-data regimes but may inhibit learning when moving
to a dataset with a broader sample scope. Overall, these pretraining
experiments provide insights into how architectural inductive biases interact
with dataset characteristics to influence translation performance.

####  4.5.4 Joint Training

The motivation behind conducting this experiment was to explore whether the
simplest approach of padding channels could enable the model to learn shared
features from both datasets which is multi-language, multi layout and multi
subject. It is important to note that the GWilliams dataset comprises only 208
channels, thus we applied padding on the channel dimensions to match the
Schoffelen’s channel dimension of 273. We are concerned that this approach
might introduce significant interference, rendering both datasets untrainable.
However, our experiment revealed that our worries were unfounded. Despite the
decrease in model performance caused by the simple channel padding method, it
is indeed effective in integrating the two datasets.

Table 8: Performance of joint training on two datasets.

Dataset | BLEU-N (%) | ROUGE-1 (%)
---|---|---
N=1 | N=2 | N=3 | N=4 | F | P | R
GWilliams | 55.13 | 52.91 | 51.28 | 49.89 | 57.33 | 64.22 | 57.17
Schoffelen | 45.12 | 39.00 | 35.48 | 33.03 | 45.15 | 46.47 | 45.49

We got 55.13 and 45.12 BLEU-1 scores on GWilliams and Schoffelen respectively,
which is 5-7 points lower compared to training solely from scratch. The
results demonstrate that, although the performance was slightly affected, the
simple data padding method was able to effectively fuse the two datasets.

However, we acknowledge that the current method we employed is rudimentary,
and future work will involve exploring more advanced techniques for dataset
integration. This includes incorporating physical information from the sensors
and employing more sophisticated approaches to merge the datasets. By doing
so, we can further enhance the performance and effectiveness of joint training
in multi-language and multi-layout neural transcribing systems.

####  4.5.5 Scaling Law

To investigate the presence of a scaling law phenomenon in our model, where
larger model sizes correspond to improved performance, we conduct an
experiment using the Whisper series models at various sizes: tiny, base,
small, medium, and large. We limit the GPU resources for training to a single
H800 card with 80 GiB of memory, which results in batch size values of 128,
64, 16, 8, and 4, from smallest to largest model.

Figure 1: Performance on different model sizes. Black numbers represent the
BLEU-1 score, green numbers are effective epochs after which the evaluation
loss does not descend. Note that each experiment runs 120 epochs.

This experiment utilizes the Schoffelen dataset and sets a maximum training
epoch of 120. The best-performing checkpoint is selected for each model size.
The results, depicted in Fig. 1, indicate that when faced with limited data,
increasing the model size leads to a consistent improvement in model
performance until a certain size is reached.

This finding suggests that scaling up the model size can be a valuable
strategy for enhancing the capabilities and effectiveness of the model in this
particular task. However, it is important to consider that larger models may
suffer from degradation due to the limited availability of data. Therefore, it
is crucial to carefully select an appropriate model size.

Additionally, we observe that the number of epochs required to achieve
convergence decreases by approximately half as the model size increases. This
observation implies that, even under GPU constraints, increasing the size of
the model significantly improves performance without imposing a prohibitive
training time penalty. Larger models are able to reach convergence more
quickly, despite their increased capacity, resulting in training durations
comparable to smaller models. By carefully balancing the bias-variance
tradeoff for each specific task and dataset, researchers can fine-tune the
model to generalize optimally within the limitations of the available data,
thereby maximizing decoding quality given the computational resources at hand.

####  4.5.6 Data augmentation

There are a lot of research works that used data augmentation in EEG/MEG
signal, Bai et al. (Bai et al., 2023) used time-wise masking. Now, we explore
three kinds of augmentation, which are masking, noise injecting, and shifting,
and compare them with baseline using Whisper-Base in Sec. 4.5.5 which does not
have data augmentation. We set three types of masking, which are time-wise,
channel-wise, and block-wise. The mask unit is [channel,40],[1,time],[1,40],
for the three masking methods respectively. We set two signal-to-noise ratios
(SNR), 0dB and 15dB, of Gaussian noise with a probability of 0.5 and 1, which
compose 4 settings. We also set 0.5 and 1 probability of data shifting to see
if moving the neural data uniformly in the range of max input duration, along
with the timestamps can help the model reduce over-fitting. Fig. 2 shows that
block masking can improve performance a lot when the mask ratio is relatively
low. However, when the masking ratio is bigger, the performance seems to
decrease quickly. It can also be noticed that time masking is much worse than
the other two mask types to our surprise. It may indicate that the redundancy
of the MEG signal is much lower at this sampling rate compared to images which
usually use a mask ratio of 75% (He et al., 2022), when the mask ratio is too
big to lose essential information, it can be hard for the model to recover
information back correctly. Shifting gains no improvements under our
experimental setting, oppositely it degrades model performance significantly.
To our surprise, adding strong noise which is 0dB can enhance the model’s
performance greatly. But for weak noise of SNR 15dB, there is no benefit in
our experimental setting. It reminds us that the MEG signal is very noisy,
thus using strong noise may let the model better eliminate noise. Therefore,
we recommend utilizing block masking with a small mask ratio and low signal-
to-noise ratio (SNR) noise for data augmentation as they have demonstrated
their effectiveness in enhancing the performance of the models..

Figure 2: Performance changes with different data augmentations. The
probability means the likelihood of adding augmentation on each data segment.

####  4.5.7 Modification

In line with (Défossez et al., 2023), we modified the whisper encoder by
replacing the first convolution layer with two separate convolution layers to
increase dimensionality and reduce temporal width. Comparing our base model
with and without the point-wise convolution layer (Table 9), we observed a
significant drop in performance without the layer. This highlights the crucial
role of the point-wise convolution layer in achieving good results by
capturing important data features and representations.

Table 9: Performance of dropping the point-wise convolution layer. - indicates
deletion of module.

Model | BLEU-N (%) | ROUGE-1 (%)
---|---|---
N=1 | N=2 | N=3 | N=4 | F | P | R
baseline | 41.65 | 36.47 | 33.51 | 31.39 | 44.86 | 46.71 | 44.74
\- conv | 34.81 | 28.11 | 24.42 | 22.01 | 37.28 | 39.91 | 37.17

####  4.5.8 Data Ratio

We explore the ratio and data used in training, locking validation set and
test set to ensure fairness. In Fig. 3, we can see that performance rises as
data is adding as expected. Higher performance may come if more data is
provided since the curve seems not have reached plateau yet.

Figure 3: Fine-tune layers and data ratio.

####  4.5.9 Layers in fine-tuning

Also in Fig. 3, it can be observed that the more layers we fine-tune, the
better the model will be. Note, that it is only on this base model, and may
not be applicable to all.

##  5 Limitations

While NeuSpeech shows promising results on benchmark MEG-to-text tasks, there
remain opportunities to strengthen the scientific validity and real-world
applicability of the framework. Future work could focus on incorporating less
structured brain signals to better validate translatability of direct human
thought, including physical sensor position data for improved accuracy, and
exploring augmented datasets to address data scarcity issues that currently
limit the model’s ability to generalize and generate completely correct novel
sentences. Addressing these limitations through collecting additional training
data and enriching the inputs could help advance this emerging field by
propelling NeuSpeech towards more precisely decoding thought from brain
signals in varied real-world scenarios.

##  6 Conclusion

This paper explores MEG-to-text using existing speech recognition training,
introducing the concept of decoding neural signal as speech. This approach
brings high performance on two datasets on standard metrics, and stretches to
decode MEG signals without using teacher forcing. We also demonstrate the
feasibility of pre-training or jointly training large neural speech models
across layouts and languages. We have also explored the effects of
augmentation and so on. Despite these advances, the results of hold-out
scenarios are relatively low due to the limited number of sentences and the
noisy nature of the MEG signal. In our future work, we will solve the problem
of limited training sets and also find more mechanisms to fuse multi-
modalities, multi-layers with physical information embeddings.

##  7 Impact Statements

While promising for advancing AI, decoding thoughts from neurodata raises
serious privacy, ethical and regulatory concerns if misapplied. We conduct
this research openly and transparently, stressing that voluntary informed
consent is vital for any future applications. Continued responsible
development with strong safeguards is needed to ensure this capability is
never abused or used to infringe on individuals’ psychological autonomy,
privacy or security without permission. As an emerging field at the interface
of neuroscience and technology, clear guidelines and oversight also require
discussion to adequately address the risks as this work progresses.

## References

  * Anumanchipalli et al. (2019) Anumanchipalli, G. K., Chartier, J., and Chang, E. F.  Speech synthesis from neural decoding of spoken sentences.  _Nature_ , 568(7753):493–498, April 2019.  ISSN 1476-4687.  doi: 10.1038/s41586-019-1119-1.  URL http://dx.doi.org/10.1038/s41586-019-1119-1.
  * Baevski et al. (2020) Baevski, A., Zhou, Y., Mohamed, A., and Auli, M.  wav2vec 2.0: A framework for self-supervised learning of speech representations.  _Advances in Neural Information Processing Systems_ , 33:12449–12460, 2020.
  * Bai et al. (2023) Bai, Y., Wang, X., Cao, Y.-p., Ge, Y., Yuan, C., and Shan, Y.  Dreamdiffusion: Generating high-quality images from brain eeg signals.  _arXiv preprint arXiv:2306.16934_ , 2023.
  * Csaky et al. (2023) Csaky, R., van Es, M. W., Jones, O. P., and Woolrich, M.  Interpretable many-class decoding for meg.  _NeuroImage_ , 282:120396, 2023.
  * Dash et al. (2020) Dash, D., Ferrari, P., and Wang, J.  Decoding imagined and spoken phrases from non-invasive neural (meg) signals.  _Frontiers in neuroscience_ , 14:290, 2020.
  * Duan et al. (2023) Duan, Y., Zhou, C., Wang, Z., Wang, Y.-K., and teng Lin, C.  Dewave: Discrete encoding of eeg waves for eeg to text translation.  In _Thirty-seventh Conference on Neural Information Processing Systems_ , 2023.  URL https://openreview.net/forum?id=WaLI8slhLw.
  * Défossez et al. (2023) Défossez, A., Caucheteux, C., Rapin, J., Kabeli, O., and King, J.-R.  Decoding speech perception from non-invasive brain recordings.  _Nature Machine Intelligence_ , 5(10):1097–1107, October 2023.  ISSN 2522-5839.  doi: 10.1038/s42256-023-00714-5.  URL http://dx.doi.org/10.1038/s42256-023-00714-5.
  * Feng et al. (2023) Feng, C., Cao, L., Wu, D., Zhang, E., Wang, T., Jiang, X., Ding, H., Zhou, C., Chen, J., Wu, H., et al.  A high-performance brain-to-sentence decoder for logosyllabic language.  2023\.
  * Ghazaryan et al. (2023) Ghazaryan, G., van Vliet, M., Saranpää, A., Lammi, L., Lindh-Knuutila, T., Hultén, A., Kivisaari, S., and Salmelin, R.  Trials and tribulations when attempting to decode semantic representations from meg responses to written text.  _Language, Cognition and Neuroscience_ , pp. 1–12, 2023.
  * Gwilliams et al. (2023) Gwilliams, L., Flick, G., Marantz, A., Pylkkänen, L., Poeppel, D., and King, J.-R.  Introducing meg-masc a high-quality magneto-encephalography dataset for evaluating natural speech processing.  _Scientific Data_ , 10(1), December 2023.  ISSN 2052-4463.  doi: 10.1038/s41597-023-02752-5.  URL http://dx.doi.org/10.1038/s41597-023-02752-5.
  * He et al. (2022) He, K., Chen, X., Xie, S., Li, Y., Dollár, P., and Girshick, R.  Masked autoencoders are scalable vision learners.  In _Proceedings of the IEEE/CVF conference on computer vision and pattern recognition_ , pp. 16000–16009, 2022.
  * Lewis et al. (2020) Lewis, M., Liu, Y., Goyal, N., Ghazvininejad, M., Mohamed, A., Levy, O., Stoyanov, V., and Zettlemoyer, L.  BART: Denoising sequence-to-sequence pre-training for natural language generation, translation, and comprehension.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pp. 7871–7880, 2020.  doi: 10.18653/v1/2020.acl-main.703.  URL https://www.aclweb.org/anthology/2020.acl-main.703.
  * Lin (2004) Lin, C.-Y.  Rouge: A package for automatic evaluation of summaries.  In _Text summarization branches out_ , pp. 74–81, 2004.
  * Liu et al. (2023) Liu, Y., Zhao, Z., Xu, M., Yu, H., Zhu, Y., Zhang, J., Bu, L., Zhang, X., Lu, J., Li, Y., Ming, D., and Wu, J.  Decoding and synthesizing tonal language speech from brain activity.  _Science Advances_ , 9(23), June 2023.  ISSN 2375-2548.  doi: 10.1126/sciadv.adh0478.  URL http://dx.doi.org/10.1126/sciadv.adh0478.
  * Metzger et al. (2022) Metzger, S. L., Liu, J. R., Moses, D. A., Dougherty, M. E., Seaton, M. P., Littlejohn, K. T., Chartier, J., Anumanchipalli, G. K., Tu-Chan, A., Ganguly, K., and Chang, E. F.  Generalizable spelling using a speech neuroprosthesis in an individual with severe limb and vocal paralysis.  _Nature Communications_ , 13(1), November 2022.  ISSN 2041-1723.  doi: 10.1038/s41467-022-33611-3.  URL http://dx.doi.org/10.1038/s41467-022-33611-3.
  * Metzger et al. (2023) Metzger, S. L., Littlejohn, K. T., Silva, A. B., Moses, D. A., Seaton, M. P., Wang, R., Dougherty, M. E., Liu, J. R., Wu, P., Berger, M. A., Zhuravleva, I., Tu-Chan, A., Ganguly, K., Anumanchipalli, G. K., and Chang, E. F.  A high-performance neuroprosthesis for speech decoding and avatar control.  _Nature_ , 620(7976):1037–1046, August 2023.  ISSN 1476-4687.  doi: 10.1038/s41586-023-06443-4.  URL http://dx.doi.org/10.1038/s41586-023-06443-4.
  * Papineni et al. (2002) Papineni, K., Roukos, S., Ward, T., and Zhu, W.-J.  Bleu: a method for automatic evaluation of machine translation.  In _Proceedings of the 40th annual meeting of the Association for Computational Linguistics_ , pp. 311–318, 2002.
  * Radford et al. (2023) Radford, A., Kim, J. W., Xu, T., Brockman, G., McLeavey, C., and Sutskever, I.  Robust speech recognition via large-scale weak supervision.  In _International Conference on Machine Learning_ , pp. 28492–28518. PMLR, 2023.
  * Schoffelen et al. (2019) Schoffelen, J.-M., Oostenveld, R., Lam, N. H. L., Uddén, J., Hultén, A., and Hagoort, P.  A 204-subject multimodal neuroimaging dataset to study language processing.  _Scientific Data_ , 6(1), April 2019.  ISSN 2052-4463.  doi: 10.1038/s41597-019-0020-y.  URL http://dx.doi.org/10.1038/s41597-019-0020-y.
  * Tang et al. (2022) Tang, J., LeBel, A., Jain, S., and Huth, A. G.  Semantic reconstruction of continuous language from non-invasive brain recordings.  September 2022.  doi: 10.1101/2022.09.29.509744.  URL http://dx.doi.org/10.1101/2022.09.29.509744.
  * Wang et al. (2020) Wang, R., Chen, X., Khalilian-Gourtani, A., Chen, Z., Yu, L., Flinker, A., and Wang, Y.  Stimulus speech decoding from human cortex with generative adversarial network transfer learning.  In _2020 IEEE 17th International Symposium on Biomedical Imaging (ISBI)_. IEEE, April 2020.  doi: 10.1109/isbi45749.2020.9098589.  URL http://dx.doi.org/10.1109/isbi45749.2020.9098589.
  * Wang & Ji (2022) Wang, Z. and Ji, H.  Open vocabulary electroencephalography-to-text decoding and zero-shot sentiment classification.  In _Proceedings of the AAAI Conference on Artificial Intelligence_ , volume 36, pp. 5350–5358, 2022.
  * Willett et al. (2021) Willett, F. R., Avansino, D. T., Hochberg, L. R., Henderson, J. M., and Shenoy, K. V.  High-performance brain-to-text communication via handwriting.  _Nature_ , 593(7858):249–254, May 2021.  ISSN 1476-4687.  doi: 10.1038/s41586-021-03506-2.  URL http://dx.doi.org/10.1038/s41586-021-03506-2.
  * Willett et al. (2023) Willett, F. R., Kunz, E. M., Fan, C., Avansino, D. T., Wilson, G. H., Choi, E. Y., Kamdar, F., Hochberg, L. R., Druckmann, S., Shenoy, K. V., and Henderson, J. M.  A high-performance speech neuroprosthesis.  January 2023.  doi: 10.1101/2023.01.21.524489.  URL http://dx.doi.org/10.1101/2023.01.21.524489.
  * Xi et al. (2023) Xi, N., Zhao, S., Wang, H., Liu, C., Qin, B., and Liu, T.  Unicorn: Unified cognitive signal reconstruction bridging cognitive signals and human language.  _arXiv preprint arXiv:2307.05355_ , 2023.
  * Zhang et al. (2023) Zhang, Q., Chen, M., Bukharin, A., He, P., Cheng, Y., Chen, W., and Zhao, T.  Adaptive budget allocation for parameter-efficient fine-tuning.  _arXiv preprint arXiv:2303.10512_ , 2023.


# Mpemba effect in anisotropically driven inelastic Maxwell gases

Apurba Biswas1,2, V. V. Prasad3, and R. Rajesh1,2 1 The Institute of
Mathematical Sciences, C.I.T. Campus, Taramani, Chennai 600113, India  2 Homi
Bhabha National Institute, Training School Complex, Anushakti Nagar, Mumbai
400094, India  3 Department of Physics, Cochin University of Science and
Technology, Kochi 682022, India<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Through an exact analysis, we show the existence of Mpemba effect in an
anisotropically driven inelastic Maxwell gas, a simplified model for granular
gases, in two dimensions. Mpemba effect refers to the couterintuitive
phenomenon of a hotter system relaxing to the steady state faster than a
cooler system, when both are quenched to the same lower temperature. The
Mpemba effect has been illustrated in earlier studies on isotropically driven
granular gases, but its existence requires non-stationary initial states,
limiting experimental realisation. In this paper, we demonstrate the existence
of the Mpemba effect in anisotropically driven granular gases even when the
initial states are non-equilibrium steady states. The precise conditions for
the Mpemba effect, its inverse, and the stronger version, where the hotter
system cools exponentially faster are derived.

##  1 Introduction

In recent times, there has been considerable interest in the Mpemba effect, a
counterintuitive phenomenon wherein a hot system equilibrates faster than a
cooler system when quenched to a low temperature. It was initially predicted
for water [1, 2]. Though many reasons have been attributed to the cause of the
Mpemba effect in water including convection [3], evaporation [4], dissolved
gases [5], supercooling [6], hydrogen bonding [7, 8, 9] and non-equipartition
of energy [10], the precise cause is still debated. One such study has even
cast doubts about the existence of the Mpemba effect in water [11]. Though the
Mpemba effect, as described in the case of water, involves a phase transition
where the final phase is ice and the initial phase water or steam, similar
Mpemba effect has been observed in other physical systems that does not
involve a phase transition. The other physical systems where this effect has
been demonstrated experimentally includes clathrate hydrates [12], magnetic
alloys [13], polylactides [14] and more recently in colloidal systems [15,
16].

Analysis of model-based systems also shows the existence of the Mpemba effect
in spin systems [17, 18, 19, 20], three state Markovian systems [21], spin
glasses [22], molecular gases in contact with a thermal reservoir [23, 24, 25]
and granular systems [26, 27, 28, 29, 30]. For the case of spin systems [17,
18, 19, 20] and three state Markovian systems [21], the initial probability
distributions describing the hot and the cold systems correspond to their
equilibrium (Boltzmann) distribution and then they are evolved to the
equilibrium of the final cold temperature following the Markovian dynamics.
The exact condition for the existence of the Mpemba effect is derived by
analysing the distance between the probability distributions during the
relaxation process. Moreover, such systems also show the existence of the
inverse Mpemba effect [21] where an initially colder system can heat up faster
than an initially warmer system, the strong Mpemba effect [18] where certain
initial states lead to an exponentially faster cooling and also exhibit
optimal heating protocols [17] in which precooling leads to faster heating.
For the case of spin glasses, two systems are prepared which are in contact
with different temperature thermal baths. The time evolution of their energy
density (instantaneous energy per spin) is analysed to demonstrate the Mpemba
effect when both the systems are quenched using a cold temperature thermal
bath. For the systems of molecular gases (elastic collisions) in contact with
a thermal bath [23, 24], the Mpemba effect is analysed using the mean kinetic
energy of the constituent molecules. For the case of molecular gas of single
species, in contact with a thermal bath [23], the Mpemba effect is due to the
coupling of mean kinetic energy with the excess kurtosis of the velocity
distribution function in the presence of non-linear viscous drag. On the other
hand, for binary mixture of molecular gases [24, 25] in contact with a
background fluid, the Mpemba effect is due to the coupling of mean kinetic
energies of the individual components of the binary gas. In both cases for
molecular gases, the Mpemba effect is illustrated only for non-stationary
initial states.

In this paper, we focus on granular gases, a dilute composition of driven
inelastic particles. Granular gases are of special interest as it is one
possible area where a strong interplay between experiment and theoretical
analysis for an interacting particle system is possible. At the same time, it
is also an example of a system that is far from equilibrium. The Mpemba effect
has been demonstrated in driven granular gases in few different contexts. To
study the Mpemba effect in granular systems, two systems are prepared at two
different granular temperatures (mean kinetic energy of particles). On
quenching to a lower temperature (by changing the driving parameters), the
Mpemba effect is said to be present if the temperatures of the two systems
cross each other while relaxing to the final stationary state. For a system of
smooth monodispersed particles [26], the Mpemba effect is achieved by the
coupling of mean kinetic energy with the excess kurtosis of the velocity
distribution function. An exact analysis was possible for the case of an
inelastic Maxwell gas, wherein the rate of collision is simplified to be
independent of relative velocity, and it was shown that there has to be non-
trivial correlations among the initial velocities of the particles to achieve
the Mpemba effect [29]. The Mpemba effect was also demonstrated for a system
of rough granular gas [27], granular gas of viscoelastic particles [28] and in
inertial suspensions of granular particles [31]. In all these analysis, the
initial states of the systems are non-stationary for the Mpemba effect to be
achieved. This is a drawback, as achieving special non-stationary states in
experiments is much harder than attractive stationary states.

To achieve the Mpemba effect in a granular system with stationary initial
conditions, a couple of systems have been put forward. Through an exact
analysis of a driven binary granular Maxwell gases [29], it was shown that the
coupling between the mean kinetic energies of the two components of the binary
gas leads to the Mpemba effect, the inverse Mpemba effect and the strong
Mpemba effect starting from steady state initial conditions. Here, a mechanism
of driving the two types of particles differently is required, which may be
difficult to achieve in practice. For a monodispersed gas in two dimensions,
it was recently shown that it is possible to achieve the Mpemba effect, its
inverse and the strong counterpart with initial stationary states provided the
driving is anisotropic (different in the two directions) [30]. This was
established based on an analysis of the Enskog-Boltzmann equation for driven
granular gases with the simplifying assumption that the velocity distribution
is a gaussian. By linearising the theory about the stationary states, it is
shown that the Mpemba effect can be achieved by simply tuning the driving
strengths, thus making it an effective system for experimental realisation of
the effect. Results from event-driven simulations are consistent with the
results from the linearised theory [30].

In this paper, we do an exact analysis of the system of monodispersed
inelastic gas with anisotropic driving based on the inelastic Maxwell model in
two dimensions. Compared to the system studied in Ref. [30] where the rate of
collision is proportional to the relative velocity, in the Maxwell gas, the
rate of collision is independent of the relative velocity. While this makes
the Maxwell gas more unrealistic, it renders it more amenable to exact
analysis, at the same time retaining the qualitative features. This
advantageous feature has been exploited in obtaining more rigorous results in
both freely cooling granular gas [32, 33, 34, 35] as well as in the velocity
distributions of driven granular gases [36, 37, 38, 39, 40, 41, 42]. The
equations for the time evolution of the relevant two point velocity
correlations for the Maxwell gas form a closed set of equations [43]. We
analyse these equations to determine the condition and the parameter regime
for the existence of the Mpemba effect. With our exact analysis of the
anisotropically driven Maxwell gas, we are able to put the results of Ref.
[30], which depended on many simplifying assumptions, on a more sound footing.
We show that the Mpemba and the inverse Mpemba effects exist for steady state
initial conditions which can be prepared by tuning the physical parameters
defining the system. In this analysis, we also demonstrate the existence of
the strong Mpemba effect where for certain specific initial steady states, the
equilibration rate is exponentially faster compared to any other initial
steady states.

The remainder of the paper is organised as follows. Section 2 contains the
definition of the model. In Sec. 3, we show that the time evolution of two
point velocity-velocity correlations do not depend on higher order
correlations and form a closed set of equations. This allows for an exact
calculation of the steady state mean kinetic energies. In Sec. 4, we define
the Mpemba effect and demonstrate its existence for an anisotropically driven
granular gas. In Sec. 5, we discuss the case where driving is limited to only
one direction. Section 6 contains the summary of results and discussion of
their various implications.

##  2  The Model

Consider a monodispersed granular gas composed of $N$ identical particles. We
label the particles by $i=1,\ldots,N$ and denote their two dimensional
velocities by $\emph{{v}}_{\emph{{i}}}=(v_{ix},v_{iy})$. These velocities
evolve in time through momentum conserving inelastic binary collisions and
external driving. A pair of particles $i$ and $j$ collide at a rate
$\lambda_{c}/N$. The factor $1/N$ in the collision rates ensures that the
total rate of collisions between $N[N-1]/2$ pairs of particles are
proportional to the system size $N$. The new velocities
$\emph{{v}}^{\prime}_{\emph{{i}}}$ and $\emph{{v}}^{\prime}_{\emph{{j}}}$ are
given by

$\displaystyle\begin{split}\boldsymbol{v}^{\prime}_{i}=\boldsymbol{v}_{i}-\alpha[(\boldsymbol{v}_{i}-\boldsymbol{v}_{j}).\boldsymbol{\hat{e}}]\boldsymbol{\hat{e}},\\\
\boldsymbol{v}^{\prime}_{j}=\boldsymbol{v}_{j}+\alpha[(\boldsymbol{v}_{i}-\boldsymbol{v}_{j}).\boldsymbol{\hat{e}}]\boldsymbol{\hat{e}},\end{split}$
(1)

where

$\displaystyle\alpha=\frac{1+r}{2},$ (2)

$r$ being the co-efficient of restitution, and $\boldsymbol{\hat{e}}$ is the
unit vector along the line joining the centres of the particles at contact. We
assume that $\boldsymbol{\hat{e}}$ takes a value uniformly from $[0,2\pi)$ for
each collision. In addition to collisions, the system evolves through external
driving. Each particle is driven at a rate $\lambda_{d}$. During a driving
event, the new velocity $\emph{{v}}^{\prime}_{\emph{{i}}}$ is given by

$\displaystyle v^{\prime}_{ix}=-r_{wx}v_{ix}+\eta_{x},\quad-1<r_{wx}\leq 1,$
$\displaystyle v^{\prime}_{iy}=-r_{wy}v_{iy}+\eta_{y},\quad-1<r_{wy}\leq 1,$
(3)

where $r_{wx}$ and $r_{wy}$ are parameters of driving and $\boldsymbol{\eta}$
is a noise taken from a fixed distribution $\Phi(\boldsymbol{\eta})$. We
denote the second moment of the noise distribution by $\sigma^{2}_{x}$ and
$\sigma^{2}_{y}$:

$\sigma^{2}_{k}=\int^{\infty}_{-\infty}d\eta_{k}\eta^{2}_{k}\Phi(\boldsymbol{\eta}),\quad
k=x,y.$ (4)

Note that $\sigma^{2}_{x}\neq\sigma^{2}_{y}$ or $r_{wx}\neq r_{wy}$
corresponds to anisotropic driving and will introduce an anisotropy in the
resultant velocity distribution of the particles. Such a driving scheme [Eq.
(3)] leads the system to a steady state and has been used extensively in
earlier studies [39, 40, 41]. The physical motivations for the form of driving
may be found in Refs. [44, 43], where positive $r_{w}$’s can be identified as
the coefficient of restitution of collisions between particle and a vibrating
wall.

The model has two simplifying features. The spatial degrees of freedom have
been neglected. This corresponds to the well-mixed limit where the spatial
correlations between particles are ignored. In addition, the collision rates
are independent of the relative velocity of the colliding particles. This
corresponds to the so called Maxwell limit.

Let $P(\boldsymbol{v},t)$ denote the probability that a randomly chosen
particle has velocity $\boldsymbol{v}$ at time $t$. Its time evolution is
given by

$\displaystyle\frac{d}{dt}P(\boldsymbol{v},t)$
$\displaystyle=\lambda_{c}\int\int\int
d\boldsymbol{\hat{e}}d\boldsymbol{v}_{1}d\boldsymbol{v}_{2}P(\boldsymbol{v}_{1},t)P(\boldsymbol{v}_{2},t)\delta(\boldsymbol{v}_{1}-\alpha[(\boldsymbol{v}_{1}-\boldsymbol{v}_{2}).\boldsymbol{\hat{e}}]\boldsymbol{\hat{e}}-\boldsymbol{v})$
$\displaystyle+\lambda_{d}\int\int
d\boldsymbol{\eta}d\boldsymbol{v}_{1}\Phi(\boldsymbol{\eta})P(\boldsymbol{v}_{1},t)\delta[-r_{w}\boldsymbol{v}_{1}+\boldsymbol{\eta}-\boldsymbol{v}]-\lambda_{c}P(\boldsymbol{v},t)-\lambda_{d}P(\boldsymbol{v},t),$
(5)

where the first and third terms on the right hand side describe the gain and
loss terms due to collisions while the second and fourth terms describe the
gain and loss terms due to driving.

##  3 Two point correlations

We are interested in the time evolution of the following two-point correlation
functions:

$\displaystyle E_{x}(t)=\frac{1}{N}\sum^{N}_{i=1}\langle
v^{2}_{ix}(t)\rangle,$ $\displaystyle
C_{x}(t)=\frac{1}{N(N-1)}\sum^{N}_{i=1}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}\langle v_{ix}(t)v_{jx}(t)\rangle,$ $\displaystyle
E_{y}(t)=\frac{1}{N}\sum^{N}_{i=1}\langle v^{2}_{iy}(t)\rangle,$
$\displaystyle
C_{y}(t)=\frac{1}{N(N-1)}\sum^{N}_{i=1}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}\langle v_{iy}(t)v_{jy}(t)\rangle,$ $\displaystyle
E_{xy}(t)=\frac{1}{N}\sum^{N}_{i=1}\langle v_{ix}(t)v_{iy}(t)\rangle,$
$\displaystyle
C_{xy}(t)=\frac{1}{N(N-1)}\sum^{N}_{i=1}\sum_{\begin{subarray}{c}j=1\\\ j\neq
i\end{subarray}}^{N}\langle v_{ix}(t)v_{jy}(t)\rangle.$ (6)

$E_{x}(t)$ and $E_{y}(t)$ denote the mean kinetic energies of the particles
along $x$ and $y$ directions respectively. $E_{xy}(t)$ denote the correlations
between $v_{x}$ and $v_{y}$ of the same particle whereas $C_{x}(t)$,
$C_{y}(t)$ and $C_{xy}(t)$ denote the velocity-velocity correlations between
pairs of particles. The time evolution for these correlation functions can be
obtained starting from Eq. (5) [39, 40, 41, 43, 29]. These may be written
compactly in a matrix form as

$\displaystyle\frac{d\boldsymbol{\tilde{\Sigma}}(t)}{dt}=\boldsymbol{\tilde{R}}\boldsymbol{\tilde{\Sigma}}(t)+\boldsymbol{\tilde{D}},$
(7)

where the column vectors $\boldsymbol{\tilde{\Sigma}}(t)$ and
$\boldsymbol{\tilde{D}}$ are given by:

$\displaystyle\boldsymbol{\tilde{\Sigma}}(t)=[E_{x}(t),E_{y}(t),E_{xy}(t),C_{x}(t),C_{y}(t),C_{xy}(t)]^{T},$
(8)
$\displaystyle\boldsymbol{\tilde{D}}=[\lambda_{d}\sigma^{2}_{x},\lambda_{d}\sigma^{2}_{y},0,0,0,0]^{T}.$
(9)

While the matrix $\boldsymbol{\tilde{R}}$ can be written for any $N$, in the
thermodynamic limit $N\rightarrow\infty$, it simplifies to

$\boldsymbol{\tilde{R}}\\!=\\!\left[\begin{array}[]{cccccc}A_{1}+A^{xx}_{4}&A_{2}&0&-A_{1}&-A_{2}&0\\\
A_{2}&A_{1}+A^{yy}_{4}&0&-A_{2}&-A_{1}&0\\\ 0&0&-A_{3}+A^{xy}_{4}&0&0&A_{3}\\\
0&0&0&2A^{x}_{5}&0&0\\\ 0&0&0&0&2A^{y}_{5}&0\\\
0&0&0&0&0&A^{x}_{5}+A^{y}_{5}\end{array}\right].$ (10)

The constants $A_{1},A_{2},A_{3},A^{ij}_{4},A^{i}_{5}$ are given by:

$\displaystyle\begin{split}&A_{1}=\frac{3}{4}\lambda_{c}\alpha^{4}-\lambda_{c}\alpha,\quad
A_{2}=\frac{\lambda_{c}\alpha^{4}}{4},\\\
&A_{3}=\lambda_{c}\alpha(1-\frac{\alpha}{2}),\quad
A^{ij}_{4}=-\lambda_{d}(1-r_{wi}r_{wj}),\\\
&A^{i}_{5}=-\lambda_{d}(1+r_{wi}),\quad\textnormal{where}~{}i,j\in(x,y).\end{split}$
(11)

In the steady state, the left-hand side of Eq. (7) equals zero. After taking
the thermodynamic limit ($N\rightarrow\infty$), we obtain the steady state
values of the different two point correlation functions as

$\displaystyle
E_{x}=\frac{\lambda_{d}\big{[}\big{(}4\lambda_{d}(1-r^{2}_{wy})+\lambda_{c}\alpha(4-3\alpha)\big{)}\sigma^{2}_{x}+\alpha^{2}\lambda_{c}\sigma^{2}_{y}\big{]}}{\mathcal{F}},$
(12) $\displaystyle
E_{y}=\frac{\lambda_{d}\big{[}\big{(}4\lambda_{d}(1-r^{2}_{wx})+\lambda_{c}\alpha(4-3\alpha)\big{)}\sigma^{2}_{y}+\alpha^{2}\lambda_{c}\sigma^{2}_{x}\big{]}}{\mathcal{F}},$
(13) $\displaystyle E_{xy}=C_{x}=C_{y}=C_{xy}=0,$ (14)

where

$\displaystyle\mathcal{F}=4\lambda_{d}^{2}(1-r^{2}_{wx})(1-r^{2}_{wy})+\alpha\lambda_{c}\lambda_{d}(4-3\alpha)(2-r^{2}_{wx}-r^{2}_{wy})+2\alpha^{2}\lambda^{2}_{c}(2-3\alpha+\alpha^{2}).$
(15)

From the structure of $\boldsymbol{\tilde{R}}$ [see Eq. (10)], it is evident
that the time evolution of velocity-velocity correlations only depend
(linearly) on other velocity-velocity correlations and do not depend on the
mean kinetic energies. Thus, if in the initial state, these correlations are
zero, then they remain zero for all times. Since we will be considering only
initial states that are stationary, the velocity-velocity correlations are
initially zero [see Eq. (14)] and will continue to remain zero for all times.

We therefore set these velocity correlations to zero and write the time
evolution for only the non-zero quantities, $E_{x}$ and $E_{y}$ as:

$\frac{d\boldsymbol{\Sigma}(t)}{dt}=~{}\boldsymbol{R}\boldsymbol{\Sigma}(t)+\boldsymbol{S},$
(16)

where

$\displaystyle\boldsymbol{\Sigma}(t)=\begin{bmatrix}E_{x}(t),E_{y}(t)\end{bmatrix}^{T},$
(17)
$\displaystyle\boldsymbol{S}=\begin{bmatrix}\lambda_{d}\sigma^{2}_{x},\lambda_{d}\sigma^{2}_{y}\end{bmatrix}^{T},$
(18)

and $\boldsymbol{R}$ is a $2\times 2$ matrix, whose entries are given by

$\displaystyle\begin{split}R_{11}=&\frac{3}{4}\lambda_{c}\alpha^{2}-\lambda_{c}\alpha-\lambda_{d}(1-r^{2}_{wx}),\quad
R_{12}=\frac{\lambda_{c}}{4}\alpha^{2},\\\
R_{22}=&\frac{3}{4}\lambda_{c}\alpha^{2}-\lambda_{c}\alpha-\lambda_{d}(1-r^{2}_{wy}),\quad
R_{21}=\frac{\lambda_{c}}{4}\alpha^{2}.\\\ \end{split}$ (19)

It is convenient to work in a different set of variables than $E_{x}(t)$ and
$E_{y}(t)$. We introduce the total energy, $E_{tot}$, and the difference in
energies, $E_{dif}$, as:

$\displaystyle E_{tot}=E_{x}+E_{y},$ (20) $\displaystyle E_{dif}=E_{x}-E_{y}.$
(21)

Note that since the driving is anisotropic, $E_{dif}\neq 0$ in general.

The time evolution equations for $E_{tot}$ and $E_{dif}$ can be expressed,
starting from Eq. (16), as

$\frac{d\boldsymbol{E}(t)}{dt}=-~{}\boldsymbol{\chi}\boldsymbol{E}(t)+\boldsymbol{D},$
(22)

where

$\displaystyle\boldsymbol{E}(t)=\begin{bmatrix}E_{tot}(t),E_{dif}(t)\end{bmatrix}^{T},$
(23)
$\displaystyle\boldsymbol{D}=\begin{bmatrix}\lambda_{d}(\sigma^{2}_{x}+\sigma^{2}_{y}),\lambda_{d}(\sigma^{2}_{x}-\sigma^{2}_{y})\end{bmatrix}^{T},$
(24)

and $\boldsymbol{\chi}$ is a $2\times 2$ matrix with the components of the
matrix given by

$\displaystyle\begin{split}\chi_{11}&=\frac{2\lambda_{c}\alpha(1-\alpha)+\lambda_{d}(2-r^{2}_{wx}-r^{2}_{wy})}{2},\quad\chi_{12}=\frac{\lambda_{d}(r^{2}_{wy}-r^{2}_{wx})}{2},\\\
\chi_{22}&=\frac{\lambda_{c}\alpha(2-\alpha)+\lambda_{d}(2-r^{2}_{wx}-r^{2}_{wy})}{2},\quad\chi_{21}=\frac{\lambda_{d}(r^{2}_{wy}-r^{2}_{wx})}{2}.\end{split}$
(25)

Equation (22) can be solved exactly by linear decomposition using the
eigenvalues $\lambda_{\pm}$ of $\boldsymbol{\chi}$:

$\lambda_{\pm}=\frac{1}{4}\Big{[}2\lambda_{d}(2-r^{2}_{wx}-r^{2}_{wy})+\alpha\lambda_{c}(4-3\alpha)\pm\sqrt{4\lambda^{2}_{d}(r^{2}_{wy}-r^{2}_{wx})^{2}+\alpha^{4}\lambda^{2}_{c}}\Big{]}.$
(26)

It is straightforward to show that $\lambda_{\pm}>0$ with
$\lambda_{+}>\lambda_{-}$. The solution for $E_{tot}(t)$ and $E_{dif}(t)$ is

$\displaystyle\begin{split}E_{tot}(t)-\langle
E_{tot}\rangle&=K_{+}e^{-\lambda_{+}t}+K_{-}e^{-\lambda_{-}t},\\\
E_{dif}(t)-\langle
E_{dif}\rangle&=L_{+}e^{-\lambda_{+}t}+L_{-}e^{-\lambda_{-}t},\end{split}$
(27)

where $\langle E_{tot}\rangle$ and $\langle E_{dif}\rangle$ are steady state
values of $E_{tot}(t)$ and $E_{dif}(t)$ respectively. The coefficients
$K_{+},K_{-},L_{+}$ and $L_{-}$ along with $\langle E_{tot}\rangle$ and
$\langle E_{dif}\rangle$ are given in Eq. (42). These coefficients depend only
on the system parameters and initial conditions. Equation (27) gives the full
time dependent solution for the energies, and we utilise them to demonstrate
the Mpemba effect.

##  4 The Mpemba effect in an anisotropically driven gas

In this section, we show and determine the conditions for the existence of the
Mpemba effect in the anisotropically driven monodispersed Maxwell gas based on
the analysis of $E_{tot}(t)$ and $E_{dif}(t)$ [see Eq. (27)]. The Mpemba
effect in granular systems has been defined as follows [26, 27, 29, 28, 30].
Consider two systems $P$ and $Q$ which have identical parameters except for
the pair of driving strengths, $\sigma^{2}_{x}$ and $\sigma^{2}_{y}$. We will
choose $E_{tot}$ of $P$ to be higher. Note that the systems $P$ and $Q$ are
initially in steady states. We denote their steady state values for the
energies by $[E^{P}_{tot}(0),E^{P}_{dif}(0)]$ and
$[E^{Q}_{tot}(0),E^{Q}_{dif}(0)]$ respectively. Both the systems are then
quenched to a common steady state having lower energy compared to the initial
steady state energies of $P$ and $Q$. This is achieved by changing the driving
strengths of $P$ and $Q$ to the common driving strengths, $\sigma^{2}_{x}$ and
$\sigma^{2}_{y}$ of the final steady state, keeping all the other parameters
of both the systems constant.

We say that the Mpemba effect is present if the two trajectories
$E^{P}_{tot}(t)$ and $E^{Q}_{tot}(t)$ cross each other at some finite time
$t=\tau$ at which

$E^{P}_{tot}(\tau)=E^{Q}_{tot}(\tau).$ (28)

To obtain the value of $\tau$, we equate the energies for $P$ and $Q$ from Eq.
(27) to obtain

$K^{P}_{+}e^{-\lambda_{+}\tau}+K^{P}_{-}e^{-\lambda_{-}\tau}=K^{Q}_{+}e^{-\lambda_{+}\tau}+K^{Q}_{-}e^{-\lambda_{-}\tau},$
(29)

whose solution is

$\tau=\frac{1}{\lambda_{+}-\lambda_{-}}\ln[\frac{K^{P}_{+}-K^{Q}_{+}}{K^{Q}_{-}-K^{P}_{-}}\Big{]}.$
(30)

In terms of the parameters of the initial steady states, $\tau$ reduces to

$\tau=\frac{1}{\lambda_{+}-\lambda_{-}}\ln[\frac{\chi_{12}\Delta
E_{dif}-(\lambda_{-}-\chi_{11})\Delta E_{tot}}{\chi_{12}\Delta
E_{dif}-(\lambda_{+}-\chi_{11})\Delta E_{tot}}\Big{]},$ (31)

where

$\displaystyle\begin{split}&\Delta E_{tot}=E^{P}_{tot}(0)-E^{Q}_{tot}(0),\\\
&\Delta E_{dif}=E^{P}_{dif}(0)-E^{Q}_{dif}(0).\end{split}$ (32)

For the Mpemba effect to be present, we require that $\tau>0$. Since
$\lambda_{+}>\lambda_{-}$, the argument of logarithm in Eq. (31) should be
greater than one. Simplifying, we obtain the criterion for the crossing of the
two trajectories to be

$\frac{\Delta E_{tot}}{\Delta
E_{dif}}<\frac{2\lambda_{d}(r^{2}_{wy}-r^{2}_{wx})}{\lambda_{c}\alpha^{2}+\sqrt{4\lambda^{2}_{d}(r^{2}_{wy}-r^{2}_{wx})^{2}+\alpha^{4}\lambda^{2}_{c}}}.$
(33)

The right hand side of Eq. (33) depends only on the intrinsic parameters of
the system and it is always less than one (since $\alpha,\lambda_{c}>0$). On
the other hand, the ratio $\Delta E_{tot}/\Delta E_{dif}$, depends on the
initial steady state energies of $P$ and $Q$ [see Eq. (32)]. In the stationary
state, the ratio $\Delta E_{tot}/\Delta E_{dif}$ is given by

$\displaystyle\frac{\Delta E_{tot}}{\Delta
E_{dif}}=\frac{\big{[}\lambda_{d}(1-r^{2}_{wy})+\alpha\lambda_{c}(2-\alpha)\big{]}\Delta\sigma^{2}_{x}+\big{[}2\lambda_{d}(1-r^{2}_{wx})+\alpha\lambda_{c}(2-\alpha)\big{]}\Delta\sigma^{2}_{y}}{2\big{[}\lambda_{d}(1-r^{2}_{wy})+\alpha\lambda_{c}(1-\alpha)\big{]}\Delta\sigma^{2}_{x}-\big{[}\lambda_{d}(1-r^{2}_{wx})+\alpha\lambda_{c}(1-\alpha)\big{]}\Delta\sigma^{2}_{y}},$
(34)

where,

$\displaystyle\Delta\sigma^{2}_{i}=(\sigma^{P}_{i})^{2}-(\sigma^{Q}_{i})^{2},\quad
i\in(x,y).$ (35)

Equation (34) shows that the ratio $\Delta E_{tot}/\Delta E_{dif}$ depends on
the intrinsic parameters of the system as well as the driving strengths,
$\sigma^{2}_{x}$ and $\sigma^{2}_{y}$. As a result, the driving strengths can
be appropriately tuned, keeping all the other intrinsic parameters identical
for both $P$ and $Q$, to prepare the initial conditions that satisfy Eq. (33).
In Fig. 1(a), we consider such a situation where Eq. (33) is satisfied. Here,
the systems $P$ and $Q$ have identical intrinsic parameters but the pair of
driving strengths, $\sigma^{2}_{x}$ and $\sigma^{2}_{y}$, are different for
the two systems. The trajectories cross at the point as predicted by Eq. (31).
It is clear that though $P$ has larger initial energy than $Q$, it relaxes
faster compared to the latter.

Figure 1: (a) The time evolution of the total energy, $E_{tot}(t)$ for
anisotropically driven systems $P$ and $Q$ of a two dimensional inelastic
Maxwell gas, driven along both the directions of the plane, with initial
conditions $E^{P}_{tot}(0)$=20.27, $E^{Q}_{tot}(0)$=17.32,
$E^{P}_{dif}(0)$=-7.93 and $E^{Q}_{dif}(0)$=6.26 such that
$E^{P}_{tot}(0)>E^{Q}_{tot}(0)$, which satisfies the condition for the Mpemba
effect as described in Eq. (33). The other parameters decribing the systems
are chosen to be $r$=0.3, $r_{wx}=0.88$ and $r_{wy}=0.39$. $P$ relaxes to the
steady state faster than $Q$, though its initial energy is larger. The time at
which the trajectories cross each other is $\tau=0.73$ as given by Eq. (31).
(b) $\Delta E_{tot}/\Delta E_{dif}$–$r$ phase diagram showing regions where
the Mpemba effect is observed and $r$ is the coefficient of restitution. The
line corresponds to a particular choice of the parameters $r_{wx}=0.2$,
$r_{wy}=0.8$, $\lambda_{c}=1.0$ and $\lambda_{d}=1.0$. The region below the
line given by Eq. (33) denotes the set of steady state initial conditions that
show the Mpemba effect whereas the region on the other side of the line
corresponds to initial states that do not show the Mpemba effect.

Figure 1(b) illustrates the phase space (initial conditions), based on Eq.
(33), where the Mpemba effect is observable. In the figure, the line denotes
the variation of right hand side of Eq. (33) with $r$, keeping all the other
intrinsic parameters of the system as constant. The figure corresponds to a
particular choice of the parameters $r_{wx}$, $r_{wy}$, $\lambda_{c}$ and
$\lambda_{d}$. If the ratio $\Delta E_{tot}/\Delta E_{dif}$ which depends on
the initial conditions of $P$ and $Q$, falls in the region below (above) the
line in the phase diagram [see Fig. 1(b)], then the system exhibits (does not
exhibit) the Mpemba effect.

For steady state initial conditions, the ratio $\Delta E_{tot}/\Delta E_{dif}$
is given by Eq. (34). As the ratio $\Delta E_{tot}/\Delta E_{dif}$ is a
function of the driving strengths, $\sigma^{2}_{x}$ and $\sigma^{2}_{y}$ [see
Eq. (34)], they can be appropriately tuned, independently for the systems $P$
and $Q$ as well as along the $x$ and $y$ directions, to access the entire
region of phase space where the Mpemba effect is observable.

Note that one can introduce anisotropy in the mean kinetic energies by simply
considering the case $\sigma^{2}_{x}\neq\sigma^{2}_{y}$, and keeping
$r_{wx}=r_{wy}$ [see Eqs. (12) and (13)]. But in that case, the condition for
the Mpemba effect reduces to $\Delta E_{tot}<0$ [see Eq. (33)] which is not
possible to realise as we have assumed $\Delta
E_{tot}=E^{P}_{tot}(0)-E^{Q}_{tot}(0)>0$. Therefore, to demonstrate the Mpemba
effect, we restrict ourselves to the case $r_{wx}\neq r_{wy}$.

###  4.1 The inverse Mpemba effect

Consider now the case where a system is heated instead of being cooled unlike
the direct Mpemba effect. Now if an initially colder system heats up faster
than a system at an intermediate one then it is called the inverse Mpemba
effect. We follow the same analysis as in the direct Mpemba effect. The
condition for the inverse Mpemba effect is same as that for the direct Mpemba
effect as given in Eq. (33). We prepare two systems $P$ and $Q$ such that $P$
has a higher initial total energy than $Q$ and also satisfy the condition for
the inverse Mpemba effect [Eq. (33)]. We then quench both the systems to a
common steady state having higher total energy compared to the initial total
energies of both $P$ and $Q$. The cross-over time $\tau$ at which the
trajectories $E^{P}_{tot}(t)$ and $E^{Q}_{tot}(t)$ cross is given by Eq. (31).
An example is illustrated in Fig. 2.

Figure 2: The time evolution of the total energy, $E_{tot}(t)$ for
anisotropically driven systems $P$ and $Q$ of a two dimensional inelastic
Maxwell gas, driven along both the directions of the plane, with initial
conditions $E^{P}_{tot}(0)$=20.27, $E^{Q}_{tot}(0)$=17.32,
$E^{P}_{dif}(0)$=-7.93 and $E^{Q}_{dif}(0)$=6.26 such that
$E^{P}_{tot}(0)>E^{Q}_{tot}(0)$, which satisfies the condition for the inverse
Mpemba effect as described in Eq. (33). The other parameters decribing the
systems are chosen to be $r$=0.3, $r_{wx}=0.88$ and $r_{wy}=0.39$. $P$ relaxes
to the steady state slower than $Q$, though its initial energy is larger. The
time at which the trajectories cross each other is $\tau=0.73$ as given by Eq.
(31).

The phase space of the initial steady states that satisfy the condition for
the inverse Mpemba effect turns out to be the same as that for the direct
Mpemba effect and is given by Eq. (34). Thus, Fig. 1(b) also illustrates the
valid initial steady states given by Eq. (34) that satisfy the condition [Eq.
(33)] where the inverse Mpemba effect is observable.

###  4.2 The strong Mpemba effect

It can be shown that there exists certain initial conditions such that the
system at higher temperature relaxes to a final steady state exponentially
faster compared to other initial conditions. This phenomenon is called the
strong Mpemba effect. The effect may be realised when the coefficient
($K_{-}$) associated with the slower relaxation rate in the time evolution of
total kinetic energy, $E_{tot}(t)$ [see Eq. (27)] vanishes.

Setting the coefficient $K_{-}$ [given by Eq. (42)] to zero, we obtain

$\displaystyle
E_{tot}(0)=\frac{2\lambda_{d}(r^{2}_{wy}-r^{2}_{wx})}{\lambda_{c}\alpha^{2}+\sqrt{4\lambda^{2}_{d}(r^{2}_{wy}-r^{2}_{wx})^{2}+\alpha^{4}\lambda^{2}_{c}}}E_{dif}(0)+c,$
(36)

where

$\displaystyle
c=\frac{2\lambda_{d}\left[2(\sigma^{2}_{x}+\sigma^{2}_{y})+\lambda_{d}(\sigma^{2}_{x}-\sigma^{2}_{y})(r^{2}_{wx}-r^{2}_{wy})\right]}{2\lambda_{d}(2-r^{2}_{wx}-r^{2}_{wy})+\alpha\lambda_{c}(4-3\alpha)-\sqrt{4\lambda^{2}_{d}(r^{2}_{wy}-r^{2}_{wx})^{2}+\alpha^{4}\lambda^{2}_{c}}}.$
(37)

The solution of Eq. (36) in terms of $E_{tot}(0)$ and $E_{dif}(0)$ yields the
set of initial states whose relaxation is exponentially faster than the set of
generic states. Among these initial states one would like to determine the
ones which are steady states. The steady state ratio of
$E_{tot}(0)/E_{dif}(0)$ [see Eq. (42)] for a system is given by

$\displaystyle\frac{E_{tot}(0)}{E_{dif}(0)}=\frac{\begin{aligned}
&\big{[}\lambda_{d}(1-r^{2}_{wy})+\alpha\lambda_{c}(2-\alpha)\big{]}\sigma^{2}_{x}+\big{[}2\lambda_{d}(1-r^{2}_{wx})+\alpha\lambda_{c}(2-\alpha)\big{]}\sigma^{2}_{y}\end{aligned}}{\begin{aligned}
&2\big{[}\lambda_{d}(1-r^{2}_{wy})+\alpha\lambda_{c}(1-\alpha)\big{]}\sigma^{2}_{x}-\big{[}\lambda_{d}(1-r^{2}_{wx})+\alpha\lambda_{c}(1-\alpha)\big{]}\sigma^{2}_{y}\end{aligned}},$
(38)

and is a function of the driving strengths, $\sigma^{2}_{x}$ and
$\sigma^{2}_{y}$, as all other parameters are kept constant. One observes that
the valid steady states with initial energies, $E_{tot}(0)$ and $E_{dif}(0)$
that satisfy the condition for the strong Mpemba effect [see Eq. (36)] can be
obtained by appropriately tuning the driving strengths.

Thus, for a system of monodispersed Maxwell gas in two dimensions, there
exists steady state initial conditions that satisfy the condition given by Eq.
(36) and hence approach the final steady state exponentially faster compared
to any other similar system whose initial energies lie slightly below or above
the line. An example of the strong Mpemba effect is shown in Fig. 3.

Figure 3: The time evolution of the total energy, $E_{tot}(t)$ for
anisotropically driven systems $P$ and $Q$ of a two dimensional inelastic
Maxwell gas, driven along both the directions of the plane, with
$E^{P}_{tot}(0)=95.89$, $E^{Q}_{tot}(0)=61.57$, $E^{P}_{diff}(0)=-59.82$ and
$E^{Q}_{diff}(0)=-4.26$ such that $E^{P}_{tot}(0)>E^{Q}_{tot}(0)$. These
initial values satisfy both the conditions for the Mpemba effect as described
in Eq. (33) as well as those for the strong Mpemba effect (for system $P$) as
described in Eq. (36). The other parameters defining the system are chosen to
be $r=0.2$, $r_{wx}=0.88$ and $r_{wy}=0.49$. $P$ equilibrates to the final
state at an exponentially faster rate compared to $Q$ and the time at which
the trajectories cross each other is $\tau=4.14$ as given by Eq. (31).

##  5 Special case when the driving is only in one direction

In Sec. 4, we discussed the possibility of the Mpemba effect in the case of
anisotropically driven monodispersed Maxwell gas where the particles are
driven along both the directions. We now consider a similar system but the
driving is restricted to one direction. We follow the same analysis as in Sec.
4. Here, for the case when particles are driven only along $x$-direction (say)
with driving strengths, $\sigma^{2}_{x}\neq 0$ and $\sigma^{2}_{y}=0$, the
time evolution of mean kinetic energies $E_{x}$ and $E_{y}$ is

$\displaystyle\begin{split}&\frac{dE_{x}(t)}{dt}=E_{x}\big{[}\lambda_{c}\alpha(\frac{3}{4}\alpha-1)-\lambda_{d}(1-r^{2}_{wx})\big{]}+E_{y}\big{[}\frac{\lambda_{c}}{4}\alpha^{2}\big{]}+\lambda_{d}\sigma^{2}_{x},\\\
&\frac{dE_{y}(t)}{dt}=E_{x}\big{[}\frac{\lambda_{c}}{4}\alpha^{2}\big{]}+E_{y}\big{[}\lambda_{c}\alpha(\frac{3}{4}\alpha-1)\big{]}.\end{split}$
(39)

The time evolution for the quantities $E_{tot}$ and $E_{dif}$ are given by Eq.
(22) but now the column matrix $\boldsymbol{D}$ takes the form

$\displaystyle\boldsymbol{D}$
$\displaystyle=\begin{bmatrix}\lambda_{d}\sigma^{2}_{x},\lambda_{d}\sigma^{2}_{x}\end{bmatrix}^{T},$
(40)

The solutions for $E_{tot}(t)$ and $E_{dif}(t)$ are obtained in the similar
way as in Eq. (27) with the coefficients $K_{+},K_{-},L_{+}$ and $L_{-}$ along
with the steady state energies $\langle E_{tot}\rangle$ and $\langle
E_{dif}\rangle$ given in Eq. (43).

We now consider two systems labeled as $P$ and $Q$ with different initial
conditions $[E^{P}_{tot}(0),E^{P}_{dif}(0)]$ and
$[E^{Q}_{tot}(0),E^{Q}_{dif}(0)]$ where $E^{P}_{tot}(0)>E^{Q}_{tot}(0)$. Both
the systems are quenched to a common steady state whose total energy is
smaller than the initial total energies of $P$ and $Q$. This is achieved by
changing the driving strengths of $P$ and $Q$ to the common driving strengths,
$\sigma^{2}_{x}\neq 0$ and $\sigma^{2}_{y}=0$ of the final steady state,
keeping all the other parameters constant for both the systems.

The condition for the Mpemba effect to be present is the same as that derived
for the more general case [see Eq. (33)]. In Fig. 4(a), we consider such a
situation where Eq. (33) is satisfied and hence the systems $P$ and $Q$ show
the Mpemba effect. The trajectories cross at the point as predicted by Eq.
(31).

In Fig. 4(b), we identify the region of phase space (initial condition) where
the Mpemba effect is observable, based on Eq. (33). In the figure, the line
denotes the variation of right hand side of Eq. (33) with $r$, keeping all the
other intrinsic parameters of the system as constant. The figure corresponds
to a particular choice of the parameters $r_{wx}$, $\lambda_{c}$ and
$\lambda_{d}$. The region below the line in the phase diagram corresponds to
the initial conditions $\Delta E_{tot}/\Delta E_{dif}$ [see Fig. 4(b)] that
show the Mpemba effect [Eq. (33)] whereas the other region does not show the
effect.

Here, we consider that the systems $P$ and $Q$ have identical intrinsic
parameters once the quench is done to the common steady state. However, these
intrinsic parameters that characterise the initial conditions of $P$ and $Q$
or equivalently $\Delta E_{tot}/\Delta E_{dif}$, could be different. As a
result, one can tune these intrinsic parameters differently for $P$ and $Q$ to
obtain initial steady states that satisfy the condition given by Eq. (33) and
hence show the Mpemba effect.

However, when the intrinsic parameters other than driving strength is kept the
same (both before and after the quench), the ratio $\Delta E_{tot}/\Delta
E_{dif}$ for initial steady states has a simple form:

$\frac{\Delta E_{tot}}{\Delta E_{dif}}=\frac{(2-\alpha)}{2(1-\alpha)}\geq
1.5.$ (41)

Note that $\alpha\in[1/2,1]$ and hence the ratio in Eq. (41) is always larger
than or equal to 1.5. However, we know from Eq. (33) that for the Mpemba
effect to be present, $\Delta E_{tot}/\Delta E_{dif}<1$. Thus, Eq. (41) does
not satisfy the required condition for the existence of the Mpemba effect. We
conclude that for initial states that correspond to steady states where $P$
and $Q$ have identical intrinsic parameters except for the driving strength,
the Mpemba effect is not possible when the driving is restricted to one
direction.

Figure 4: (a) The time evolution of the total energy, $E_{tot}(t)$ for
anisotropically driven systems $P$ and $Q$ of a two dimensional inelastic
Maxwell gas driven along a single direction with initial conditions
$E^{P}_{tot}(0)$=28, $E^{Q}_{tot}(0)$=22, $E^{P}_{dif}(0)$=26 and
$E^{Q}_{dif}(0)$=5 such that $E^{P}_{tot}(0)>E^{Q}_{tot}(0)$, which satisfies
the condition for the Mpemba effect as described in Eq. (33). The other
parameters decribing the systems are chosen to be $r$=0.5 and $r_{wx}=0.6$.
$P$ relaxes to the steady state faster than $Q$, though its initial energy is
larger. The time at which the trajectories cross each other is $\tau=1.07$ as
given by Eq. (31). (b) $\Delta E_{tot}/\Delta E_{dif}$–$r$ phase diagram
showing regions where the Mpemba effect is observed and $r$ is the coefficient
of restitution. The line corresponds to a particular choice of the parameters
$r_{wx}=0.5$, $\lambda_{c}=1.0$ and $\lambda_{d}=1.0$. The region below the
line given by Eq. (33) shows the Mpemba effect whereas the region on the other
side of the line does not show the Mpemba effect.

##  6 Summary and discussion

In this paper, we have shown an exact analysis for the existence of the Mpemba
effect, the inverse Mpemba effect and the strong Mpemba effect in an
anisotropically driven inelastic Maxwell gas in two dimensions. The Maxwell
model for granular gases is a simplified model where the rate of collision
between the granular particles is independent of their relative velocities. In
addition, we assumed the well-mixed limit such that the spatial correlations
were ignored. The model allows for an exact solution as the two-point velocity
correlations form a coupled set of linear equations.

We show that anisotropic driving leads to the existence of the Mpemba effect
starting from steady state initial conditions unlike the case of isotropic
driving which required the initial conditions to be non-stationary. We
considered two different cases of anisotropic driving in two dimensions: when
particles are driven along one direction only and the other case where
particles are driven along both the directions. In both the cases, we show
that the Mpemba effect can exist for initial conditions which are valid steady
states characterised by the parameters of the system. We also demonstrated the
existence of the inverse Mpemba effect where a system is heated instead of
being cooled. Here, an initially colder system equilibrates to a final high
temperature state faster than an initially warmer system. We also derived the
condition for the existence of the strong Mpemba effect where for certain
initial states, a system equilibrates at an exponentially faster rate compared
to any other states.

Our results place the results obtained by us in an earlier work [30] on an
anisotropically driven granular gas with a more realistic velocity dependent
collision rate, on a more rigorous footing. One difference between results for
the two models is that in the Maxwell gas, the results depend on whether the
anisotropy is applied through different driving strengths
($\sigma^{2}_{x}\neq\sigma^{2}_{y}$) or through different driving parameters
($r_{wx}\neq r_{wy}$). For the Maxwell gas, if all other parameters are kept
the same, then the latter condition is necessary. This difference is due to
the lack of non-linear coupling between $E_{tot}$ and $E_{dif}$ in the Maxwell
model.

##  Appendix A Coefficients for the time evolutions of $E_{tot}$ and $E_{dif}$

In this Appendix, we solve for the time evolutions of $E_{tot}$ and $E_{dif}$
for the anisotropically driven inelastic Maxwell gas in two dimensions. We
consider two cases of anisotropic driving: when the external driving is
applied along both the directions with different driving strengths and another
case where the driving is along one direction only as described in A.1 and A.2
respectively.

###  A.1 Driven along both directions of the plane

We consider an inelastic Maxwell gas in two dimensions where the particles are
driven along both directions of the plane at rate $\lambda_{d}$ and with
different driving strengths, $\sigma^{2}_{x}$ and $\sigma^{2}_{y}$
respectively. The time evolutions of $E_{tot}$ and $E_{dif}$ are given as in
Eq. (27) where the coefficients $K_{+},K_{-},L_{+}$ and $L_{-}$ along with the
steady state energies $\langle E_{tot}\rangle$ and $\langle E_{dif}\rangle$
are given by

$\displaystyle K_{+}$
$\displaystyle=\frac{1}{\gamma}\Big{[}-(\lambda_{-}-\chi_{11})E_{tot}(0)+\chi_{12}E_{dif}(0)-\frac{\lambda_{d}}{\lambda_{+}}\big{[}\big{(}\chi_{12}-(\lambda_{-}-\chi_{11})\big{)}\sigma^{2}_{x}-\big{(}\chi_{12}+(\lambda_{-}-\chi_{11})\big{)}\sigma^{2}_{y}\big{]}\Big{]},$
$\displaystyle K_{-}$
$\displaystyle=\frac{1}{\gamma}\Big{[}(\lambda_{+}-\chi_{11})E_{tot}(0)-\chi_{12}E_{dif}(0)+\frac{\lambda_{d}}{\lambda_{-}}\big{[}\big{(}\chi_{12}-(\lambda_{+}-\chi_{11})\big{)}\sigma^{2}_{x}-\big{(}\chi_{12}+(\lambda_{+}-\chi_{11})\big{)}\sigma^{2}_{y}\big{]}\Big{]},$
$\displaystyle\langle E_{tot}\rangle$
$\displaystyle=\frac{\lambda_{d}}{\gamma}\Big{[}\frac{\big{(}\chi_{12}-(\lambda_{-}-\chi_{11})\big{)}\sigma^{2}_{x}-\big{(}\chi_{12}+(\lambda_{-}-\chi_{11})\big{)}\sigma^{2}_{y}}{\lambda_{+}}-\frac{\big{(}\chi_{12}-(\lambda_{+}-\chi_{11})\big{)}\sigma^{2}_{x}-\big{(}\chi_{12}+(\lambda_{+}-\chi_{11})\big{)}\sigma^{2}_{y}}{\lambda_{-}}\Big{]},$
$\displaystyle L_{+}$
$\displaystyle=\frac{1}{\gamma}\Big{[}-\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}}E_{tot}(0)+(\lambda_{+}-\chi_{11})E_{dif}(0)-\frac{\lambda_{d}}{\lambda_{+}\chi_{12}}\big{[}(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})(\sigma^{2}_{x}-\sigma^{2}_{y})\big{]}\Big{]},$
$\displaystyle L_{-}$
$\displaystyle=\frac{1}{\gamma}\Big{[}\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}}E_{tot}(0)-(\lambda_{+}-\chi_{11})E_{dif}(0)+\frac{\lambda_{d}}{\lambda_{-}\chi_{12}}\big{[}(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})(\sigma^{2}_{x}-\sigma^{2}_{y})\big{]}\Big{]},$
$\displaystyle\langle E_{dif}\rangle$
$\displaystyle=\frac{\lambda_{d}}{\chi_{12}\gamma}\Big{[}(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})(\sigma^{2}_{x}-\sigma^{2}_{y})\big{(}\frac{1}{\lambda_{+}}-\frac{1}{\lambda_{-}}\big{)}\Big{]},$
$\displaystyle\gamma$ $\displaystyle=\lambda_{+}-\lambda_{-}.$ (42)

###  A.2 Driven along a single direction of the plane

Here, we consider an inelastic Maxwell gas in two dimensions where the
particles are driven along a single direction (say along $x$ direction) at a
rate $\lambda_{d}$ and with driving strengths, $\sigma^{2}_{x}\neq 0$ and
$\sigma^{2}_{y}=0$. The time evolutions of $E_{tot}$ and $E_{dif}$ are given
as in Eq. (27) where the coefficients $K_{+},K_{-},L_{+}$ and $L_{-}$ along
with the steady state energies $\langle E_{tot}\rangle$ and $\langle
E_{dif}\rangle$ are now given by

$\displaystyle\begin{split}K_{+}&=\frac{1}{\gamma}\Big{[}(-\lambda_{-}+\chi_{11})E_{tot}(0)+\chi_{12}E_{dif}(0)-\frac{\chi_{12}-\lambda_{-}+\chi_{11})}{\lambda_{+}}\lambda_{d}\sigma^{2}_{x}\Big{]},\\\
K_{-}&=\frac{1}{\gamma}\Big{[}(\lambda_{+}-\chi_{11})E_{tot}(0)-\chi_{12}E_{dif}(0)+\frac{\chi_{12}-\lambda_{+}+\chi_{11})}{\lambda_{-}}\lambda_{d}\sigma^{2}_{x}\Big{]},\\\
\langle
E_{tot}\rangle&=\frac{1}{\gamma}\Big{[}\frac{\chi_{12}-(\lambda_{-}-\chi_{11})}{\lambda_{+}}-\frac{\chi_{12}-(\lambda_{+}-\chi_{11})}{\lambda_{-}}\Big{]}\lambda_{d}\sigma^{2}_{x},\\\
\langle
E_{dif}\rangle&=\frac{1}{\gamma}\Big{[}\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}\lambda_{+}}-\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}\lambda_{-}}\Big{]}\lambda_{d}\sigma^{2}_{x},\\\
L_{+}&=\frac{1}{\gamma}\Big{[}-\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}}E_{tot}(0)+(\lambda_{+}-\chi_{11})E_{dif}(0)-\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}\lambda_{+}}\lambda_{d}\sigma^{2}_{x}\Big{]},\\\
L_{-}&=\frac{1}{\gamma}\Big{[}\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}}E_{tot}(0)-(\lambda_{-}-\chi_{11})E_{dif}(0)+\frac{(\lambda_{+}-\chi_{11})(\lambda_{-}-\chi_{11})}{\chi_{12}\lambda_{-}}\lambda_{d}\sigma^{2}_{x}\Big{]},\\\
\gamma&=\lambda_{+}-\lambda_{-}.\end{split}$ (43)

## References

  * [1] Lee H Meteorologica, Loeb classical library: Greek authors p.87 (Harvard University Press, Cambridge, MA, 1952)
  * [2] Mpemba E B and Osborne D G 1969 Phys. Educat. 4 172–175
  * [3] Vynnycky M and Kimura S 2015 Int. J. Heat Mass Transf. 80 243 – 255
  * [4] Mirabedin S M and Farhadi F 2017 Int. J. Refrig. 73 219 – 225
  * [5] Katz J I 2009 Am. J. Phys. 77 27–29
  * [6] Auerbach D 1995 Am. J. Phys. 63 882–885
  * [7] Zhang X, Huang Y, Ma Z, Zhou Y, Zhou J, Zheng W, Jiang Q and Sun C Q 2014 Phys. Chem. Chem. Phys. 16(42) 22995–23002
  * [8] Tao Y, Zou W, Jia J, Li W and Cremer D 2017 J. Chem. Theory Comput. 13 55–76
  * [9] Jin J and Goddard III W A 2015 J. Phys. Chem. C 119 2622–2629
  * [10] Gijón A, Lasanta A and Hernández E 2019 Phys. Rev. E 100 032103
  * [11] Burridge H C and Linden P F 2016 Sci. Rep. 6 37665
  * [12] Ahn Y H, Kang H, Koh D Y and Lee H 2016 Korean J. Chem. Eng. 33 1903–1907
  * [13] Chaddah P, Dash S, Kumar K and Banerjee A 2010 arXiv preprint arXiv:1011.3598
  * [14] Hu C, Li J, Huang S, Li H, Luo C, Chen J, Jiang S and An L 2018 Cryst. Growth Des. 18 5757–5762
  * [15] Kumar A and Bechhoefer J 2020 Nature 584 64–68
  * [16] Kumar A, Chetrite R and Bechhoefer J 2021 arXiv preprint arXiv:2104.12899
  * [17] Gal A and Raz O 2020 Phys. Rev. Lett. 124(6) 060602
  * [18] Klich I, Raz O, Hirschberg O and Vucelja M 2019 Phys. Rev. X 9(2) 021060
  * [19] Klich I and Vucelja M 2018 arXiv preprint arXiv:1812.11962
  * [20] Das S K and Vadakkayil N 2021 Phys. Chem. Chem. Phys.
  * [21] Lu Z and Raz O 2017 Proc. Natl. Acad. Sci. USA 114 5083–5088
  * [22] Baity-Jesi M, Calore E, Cruz A, Fernandez L A, Gil-Narvión J M, Gordillo-Guerrero A, Iñiguez D, Lasanta A, Maiorano A, Marinari E et al. 2019 Proc. Natl. Acad. Sci. USA 116 15350–15355
  * [23] Santos A and Prados A 2020 Phys. Fluids 32 072010
  * [24] González R G, Khalil N and Garzó V 2020 arXiv preprint arXiv:2010.14215
  * [25] González R G and Garzó V 2020 arXiv preprint arXiv:2011.13237
  * [26] Lasanta A, Vega Reyes F, Prados A and Santos A 2017 Phys. Rev. Lett. 119(14) 148001
  * [27] Torrente A, López-Castaño M A, Lasanta A, Reyes F V, Prados A and Santos A 2019 Phys. Rev. E 99(6) 060901
  * [28] Mompó E, Castaño M, Torrente A, Reyes F V and Lasanta A 2020 arXiv preprint arXiv:2006.00241
  * [29] Biswas A, Prasad V V, Raz O and Rajesh R 2020 Phys. Rev. E 102(1) 012906
  * [30] Biswas A, Prasad V V and Rajesh R 2021 arXiv preprint arXiv:2104.08730
  * [31] Takada S, Hayakawa H and Santos A 2021 Phys. Rev. E 103(3) 032901
  * [32] Ben-Naim E and Krapivsky P L 2002 Phys. Rev. E 66(1) 011309
  * [33] Baldassarri A, Marconi U M B and Puglisi A 2002 Europhys. Lett. 58 14–20
  * [34] Ernst M H and Brito R 2002 Europhys. Lett. 58 182–187
  * [35] Krapivsky P L and Ben-Naim E 2002 J. Phys. A 35 L147–L152
  * [36] Ernst M H and Brito R 2002 Phys. Rev. E 65(4) 040301
  * [37] Antal T, Droz M and Lipowski A 2002 Phys. Rev. E 66(6) 062301
  * [38] Santos A and Ernst M H 2003 Phys. Rev. E 68(1) 011305
  * [39] Prasad V V, Sabhapandit S and Dhar A 2014 Europhys. Lett. 104 54003
  * [40] Prasad V V, Sabhapandit S and Dhar A 2014 Phys. Rev. E 90(6) 062130
  * [41] Prasad V V, Das D, Sabhapandit S and Rajesh R 2017 Phys. Rev. E 95 032909
  * [42] Biswas A, Prasad V V and Rajesh R 2020 J. Stat. Mech.: Theory Exp. 2020(1) 013202
  * [43] Prasad V V, Das D, Sabhapandit S and Rajesh R 2019 J. Stat. Mech.: Theory Exp. 2019(6) 063201
  * [44] Prasad V V and Rajesh R 2019 J. Stat. Phys. 176 1409–1433