#  MUSER: MUltimodal Stress Detection using
Emotion Recognition as an Auxiliary Task

Yiqun Yao${}^{\textrm{{1}}}$, Michalis Papakostas${}^{\textrm{{1}}}$, Mihai
Burzo${}^{\textrm{{2}}}$, Mohamed Abouelenien${}^{\textrm{{3}}}$, Rada
Mihalcea${}^{\textrm{{1}}}$
${}^{\textrm{1}}$Computer Science and Engineering, University of Michigan
${}^{\textrm{2}}$Mechanical Engineering, University of Michigan
${}^{\textrm{3}}$Computer and Information Science, University of Michigan
<EMAIL_ADDRESS>

###### Abstract

The capability to automatically detect human stress can benefit artificial
intelligent agents involved in affective computing and human-computer
interaction. Stress and emotion are both human affective states, and stress
has proven to have important implications on the regulation and expression of
emotion. Although a series of methods have been established for multimodal
stress detection, limited steps have been taken to explore the underlying
inter-dependence between stress and emotion. In this work, we investigate the
value of emotion recognition as an auxiliary task to improve stress detection.
We propose MUSER – a transformer-based model architecture and a novel multi-
task learning algorithm with speed-based dynamic sampling strategy.
Evaluations on the Multimodal Stressed Emotion (MuSE) dataset show that our
model is effective for stress detection with both internal and external
auxiliary tasks, and achieves state-of-the-art results.

##  1 Introduction

Stress is a feeling of emotional or physical tension, as a response to the
environment when people have difficulty dealing with the conditions Dobson and
Smith (2000); Muthukumar and Nachiappan (2013). Stress detection is a
classification task that predicts whether a certain target is under stress.
The task has drawn research attention for two reasons: first, stress detection
plays an important role in applications related to psychological well-being
Cohen et al. (1991), cognitive behavior therapies Tull et al. (2007), and safe
driving Gao et al. (2014); Chen et al. (2017); second, stress is a known
regulator of human emotion mechanisms Tull et al. (2007), and thus research on
stress detection can potentially benefit the development of emotionally
intelligent agents.

The impact of stress on human behavior can be observed through various
modalities. Previous work has considered both unimodal and multimodal stress
detection using acoustic, video and physiological sensor signals Lane et al.
(2015); Jaques et al. (2016); Aigrain et al. (2016); Alberdi et al. (2016);
Bara et al. (2020). However, text-based stress detection remains vastly
underexplored, with some studies Lin et al. (2014) showing the potential for
further research. In recent years, the surge of advanced natural language
understanding models and structures provides a great opportunity for stress
detection systems, especially using the textual modality. In this work, we
focus on the textual and acoustic modalities. For the model architecture, we
use Transformers Vaswani et al. (2017) as a textual encoder and Multi-Layer
Perceptrons (MLP) as an acoustic encoder.

The majority of existing stress detection methods are based on single-task
learning with the binary stress/non-stress labels. However, stress is not an
isolated affective state, but closely related to the expression and regulation
of human emotions. Physiological studies Wang and Saudino (2011) have
demonstrated that emotion and stress share some neural structures, including
prefrontal cortex Taylor and Stanton (2007), anterior cingulate cortex
Pruessner et al. (2008), and amygdala Adolphs (2003). Acoustic studies
Paulmann et al. (2016) have shown that the pitch and amplitude of human
emotional prosody is different under stress and non-stressed status. Inspired
by these studies, our work aims to exploit the inter-dependence between
emotion and stress. Specifically, we investigate the value of emotion
recognition as an auxiliary task for stress detection.

Multi-task learning Pasunuru and Bansal (2017); Gottumukkala et al. (2020);
Guo et al. (2018a); Gong et al. (2019) has proven to be effective for
transferring knowledge between different tasks. Dynamic sampling strategies,
which aim at adaptively adjusting the ratio of samples from different tasks,
are widely used to balance the training schedule. However, strategies based on
gradients Chen et al. (2018b), uncertainty Kendall et al. (2018) or loss Liu
et al. (2019) cannot leverage the validation performances, while some
performance-based strategies Gottumukkala et al. (2020) are impractical if the
metrics for different tasks are not directly comparable (i.e., with different
scale ranges). To this end, we propose a novel speed-based strategy that is
both effective and efficient in the multi-task learning for stress and
emotion.

Our method is evaluated on the Multimodal Stressed Emotion (MuSE) dataset
Jaiswal et al. (2019, 2020), which includes both stress and emotion labels,
making it the ideal benchmark for an in-depth analysis of their inter-
dependence. To test the generalization ability of our method, we also use an
external emotion dataset for the auxiliary task. Multimodal emotion
recognition is a well-studied field with many existing datasets Busso et al.
(2008, 2016); Chen et al. (2018a); Barros et al. (2018); Zadeh et al. (2018).
We choose the OMG-Emotion dataset Barros et al. (2018) as the external
auxiliary task because it is representative and challenging, with numerical
emotion scores instead of categorical labels.

Our paper makes four main contributions. First, we show the inter-dependence
between stress and emotion via quantitative analyses on linguistic and
acoustic features, and propose to use emotion recognition as an auxiliary task
for stress detection. Second, we establish a stress detection model with a
transformer structure, as well as a novel speed-based dynamic sampling
strategy for multi-task learning. We name our framework the MUltimodal Stress
Detector with Emotion Recognition (MUSER). Third, we achieve state-of-the-art
results on the MuSE dataset via multi-task training with stress and emotion
labels. We also achieve competitive results when we use the OMG-Emotion Barros
et al. (2018) dataset as an external auxiliary task. Finally, experimental
results show that our speed-based dynamic sampling significantly outperforms
other widely-used methods.

##  2 Related Work

###  2.1 Unimodal Stress Detection

Stress detection based on textual modality has been studied by Lin et al.
(2014) and Jaiswal et al. (2020), using the Linguistic Inquiry and Word Count
(LIWC) lexicon Pennebaker et al. (2001) to extract features that are
indicative of human emotion. Acoustic features Lane et al. (2015); Paulmann et
al. (2016); Horvath (1982); Lech and He (2014) have also been used for
unimodal stress detection in both physiological and computational studies.

A drawback of the unimodal approaches is that they only have access to partial
information about the expression of stress, while multiple modalities can
potentially be informative at the same time Aigrain et al. (2016). As
demonstrated by previous work on human sentiment and emotion prediction Zadeh
et al. (2016, 2018); Yao et al. (2020), multimodal features usually results in
better performances.

###  2.2 Multimodal Stress Detection

Commonly-used modalities for stress detection include video, audio, text and
physiological signals such as thermal maps from sensors Aigrain et al. (2016);
Alberdi et al. (2016); Lane et al. (2015); Jaques et al. (2016).

Jaiswal et al. Jaiswal et al. (2020) proposed the Multimodal Stressed Emotion
(MuSE) dataset, which includes records from all the commonly-used modalities.
Each video clip is annotated for both stress detection and emotion
recognition. Unimodal and multimodal baselines are provided for each task.
Bara et al. Bara et al. (2020) developed a multimodal deep learning method
that learns modality-independent representations in an unsupervised approach.
However, none of these models leverage the intrinsic connections between
stress and emotion.

Our experiments are conducted on the MuSE dataset using only the textual and
acoustic modalities, to be compatible with most external emotion recognition
tasks. However, our proposed multi-task learning method is model-agnostic and
can be generalized to any structure and any modality combinations.

###  2.3 Emotion Recognition

Widely-used multimodal emotion recognition datasets include SEMAINE McKeown et
al. (2011), IEMOCAP Busso et al. (2008), MOSEI Zadeh et al. (2018) and OMG-
Emotion Barros et al. (2018). Emotion can be annotated either with pre-defined
emotion categories or through two-dimensional scores of activation (arousal)
and valence, according to the self-assessment manikin proposed by Bradley and
Lang (1994). MuSE, in particular, has emotion annotations expressed by
activation and valence scores (1$\sim$9), which is more fine-grained than
categorical definitions (happy, angry, etc.). The OMG-Emotion dataset we use
as external auxiliary task is annotated in the same way with a score range of
0$\sim$1.

###  2.4 Multi-task Learning

Because of the different task natures, balancing the training procedure with
all the tasks is a critical problem for multi-task learning. Loss-balancing
strategies Chen et al. (2018b); Kendall et al. (2018); Liu et al. (2019); Gong
et al. (2019); Guo et al. (2018a); Lample et al. (2017); Yao et al. (2019) are
suitable for situations in which there are multiple training objectives that
can be combined via weighted summation for each data point. In contrast, for
multi-task learning across different datasets, a sampling strategy should be
applied to decide the mixing ratio (how many batches to sample from each task)
in each epoch. To this end, Pasunuru et al. Pasunuru and Bansal (2017) used a
fixed sampling ratio; Guo et al. Guo et al. (2018b) proposed a dynamic
sampling strategy based on reinforcement learning, which depends on the
estimation of Q-values; Gottumukkala et al. Gottumukkala et al. (2020) used a
dynamic sampling procedure based on the gap between multi-task and single-task
results – a performance-based method that requires all the tasks to use the
same set of evaluation metrics. For comparison, our proposed strategy is also
based on how fast the model is learning each task, but does not require the
metrics to be directly comparable.

Table 1: LIWC features that have the top 20 highest regression coefficients in
all three tasks.

Feature | Examples
---|---
bio | eat, blood, pain
health | clinic, flu, pill
relativity | area, stop, exit
body | cheek, hands, spit
ingest | dish, eat, pizza
positive-emo | love, nice, sweet
space | down, in, thin
time | end, until, season
perceptual | observe, hear, feeling

##  3 Expressions of Stress in Data

###  3.1 Dataset

The MuSE dataset Jaiswal et al. (2020) is collected from the multimodal video
recordings of 28 student participants, 9 female and 19 male. Each participant
is invited to a video-recording session before and after the final exam
period; sessions before exams are labeled as stressed, and the remainng ones
are labeled as non-stressed. We use only the records from the monologue sub-
sessions where both acoustic and textual modalities are available. In these
sub-sessions, the participants view five emotion-eliciting questions on the
screen in a sequence, and answer them with monologues. After each monologue,
the participants provide self-assessment scores for activation (calm vs.
excited) and valence (negative vs. positive). The scores range from 1$\sim$9\.
The monologues are segmented into sentences for pre-processing; each sentence
is annotated with the same stress label and emotion scores as the whole
monologue.

We use a train, validation, and test split of 1,853, 200, and 273 sentences,
respectively. Textual features come from the automatic transcripts for the
audio clips of each sentence. Although the sentences come with visual and
thermal features as well, we focus mainly on the textual and acoustic
modalities because this allows us to use almost any external emotion
recognition dataset as our auxiliary task.

Table 2: LIWC features that are among the top 20 highest regression
coefficients unique to stress, activation and valence tasks.

Feature | Examples
---|---
nonfl | er, hmm, umm
affect | happy, cried, abandon
social | mate, talk, child
family | daughter, husband, aunt
past | went, ran, had
money | audit, cash, owe
tentat | maybe, perhaps, guess
feel | feels, touch
sad | crying, grief, sad
negate | no, not, never
anger | hate, kill, annoyed
achieve | earn, hero, win
quant | few, many, much

Table 3: Opensmile eGeMap features that have the top 20 highest regression
coefficients for different tasks.

Task | Feature | Category | Examples
---|---|---|---
All | Amplitude | Energy | F1/F2/F3 mean amplitude
Loudness | Energy | mean loudness
Spectral Flux | Spectral | mean flux in voiced regions
Only in stress | Hammarberg Index | Spectral | hammarberg index in unvoiced regions
Alpha Ratio | Spectral | mean ratio for unvoiced regions
Slope | Spectral | 500-1500Hz in unvoiced regions
Only in activation | HNR | Energy | mean HNR (Harmonics-to-Noise Ratio)
Voice Length | Temporal | mean voiced segment length
Only in valence | Pitch | Frequency | F0 semitone
Formant | Frequency | frequency of formant 3

###  3.2 Characteristics of Stress in Language

In order to analyze the connections between linguistic features that are most
indicative of stress, activation, and valence, we first extract a feature
vector based on the LIWC lexicon Pennebaker et al. (2001). Each dimension of
the vector corresponds to a certain word category and has a value equal to the
count of words observed in that category in the sentence. We then apply
z-normalization on each feature and fit a linear model to predict the
stress/non-stress label, as well as the activation and valence scores. For
each of the three tasks, we pick the features with the top 20 highest absolute
values of the linear classification/regression coefficients, which we assume
to be the key indicators.

Features that appear in top 20 for all three tasks are shown in Table 1. The
features are ranked by the absolute value of their linear coefficients. As
shown, the “positive-emotion” and “perceptual” word classes are critical for
both emotion and stress tasks, which is intuitive because they are a pair of
inter-dependent human affect status. Bio, health, and body words are also on
the list, suggesting that both stress and emotion are closely related to
physiological status and feelings, which is potentially because they share
some neural structures in brain Wang and Saudino (2011). The intersection of
all the three top-indicator sets has nine elements, reflecting a reasonable
overlap.

Table 2 shows the word classes appearing uniquely in the top 20 indicator list
for each task. It is worth noticing that the non-fluent words (er, hmm, etc.)
are the strongest unique indicator of stress, which reflects the differences
in the audio speeches under stressed/non-stressed conditions. We could also
observe that activation is more connected to entities and events, while
valence is more related to personal feelings.

###  3.3 Characteristics of Stress in Speech

For stress indicators in the acoustic modality, we extract 88-dimensional
features using OpenSmile Eyben et al. (2010) with the eGeMaps Eyben et al.
(2015) configuration. We follow Jaiswal et al. (2020) to do speaker-level
z-normalization on each feature, and fit a linear classification/regression
model as we did for the textual features.

Table 3 shows the most indicative acoustic feature classes for all the tasks,
as well as the ones that are unique for each task. Amplitude/loudness is the
strongest indicator class for all tasks, followed by spectral flux, which is a
measure of how quickly the power spectrum of a signal is changing. It also
suggests that stress has a closer relationship with spectral features such as
slope, describing how quickly the spectrum of an audio sound tails off towards
the high frequencies.

The intersection of all three indicator sets has 11 elements, suggesting that
they share many acoustic patterns. For more detailed explanations and examples
of the eGeMaps features please refer to Eyben et al. (2015) and Botelho et al.
(2019).

Regarding the differences in the task nature, as seen in Table 4, the number
of unique indicators for each each and for each modality show that the
activation task is less independent of the stress task than the valence task.
In other words, the activation task has more indicators in common with the
stress task.

Table 4: Number of unique textual and acoustic indicators for stress, activation and valence. Feature | Stress | Activation | Valence
---|---|---|---
LIWC | 4 | 3 | 6
eGeMaps | 7 | 6 | 8
Figure 1: Multimodal fusion architecture for MUSER.

##  4 Method

###  4.1 Auxiliary Tasks

Based on the task inter-dependency demonstrated in Section 3.2 and 3.3, we
propose to use multimodal emotion recognition as an auxiliary task for stress
detection. Since MuSE has both stress and emotion labels, their activation and
valence scores can be used as an internal auxiliary task.

To test the generalization capability of our multi-task learning method, we
choose OMG-Emotion Barros et al. (2018) as an external emotion recognition
dataset for the auxiliary task, which is annotated in the same manner as MuSE
(activation/valence). We download the videos from the training and validation
sets and filter out all the samples where the video link is broken or the
length of automatic transcription is less than 5, resulting in 1,484 videos.
The contents and scenarios in the OMG-Emotion dataset are completely different
from MuSE. We hold out 300 videos as a validation set to enable dynamic
sampling.

Note that stress detection is a binary classification task, while the two
auxiliary emotion tasks have a regressive nature.

###  4.2 Pre-processing

Each utterance in the MuSE dataset is automatically segmented into sentences,
transcribed, and tokenized by a pre-trained BERT tokenizer Devlin et al.
(2019). For the acoustic modality, we use OpenSmile Eyben et al. (2010) with
the eGeMAPS configuration Eyben et al. (2015) to extract 88 utterance-level
statistical features. Following Jaiswal et al. (2020), we perform speaker-
level z-normalization on all acoustic features. For videos in the OMG-Emotion
dataset, we first extract the audio and automatic transcripts, and then do the
same pre-processing as on MuSE.

###  4.3 MUSER: Architecture

We propose MUSER: MUltimodal Stress Detector using Emotion Recognition. The
model structure is based on neural networks. Specifically, we use a
Transformer Vaswani et al. (2017) textual encoder pre-trained with BERT Devlin
et al. (2019), and an MLP-based acoustic encoder to generate representations
on each modality, and fuse them before classification or regression. Our model
architecture is depicted in Figure 1.

Figure 2: Dynamic sampling procedure for MUSER multi-task training.

####  4.3.1 Textual Encoder

For the textual encoder, we use a Transformer neural network pre-trained with
BERT on BookCorpus and English Wikipedia Devlin et al. (2019). Our Transformer
model has 12 layers, 12 attention heads, and 768-dimensional hidden states.
The averaged hidden states on the top level are projected to 256-dimensional
representations by a fully-connected layer.

####  4.3.2 Acoustic Encoder

Our acoustic encoder is a Multi-layer Perceptron network with four hidden
layers and ReLU activation. The input of the acoustic encoder is the OpenSmile
features extracted from the audio speech of each sentence, and the output of
each hidden layer is 256-dimensional.

####  4.3.3 Multimodal Fusion

We fuse the multimodal features by concatenating the top-level 256-dimensional
textual and acoustic representations. For the emotion recognition tasks, the
concatenated representation is fully connected to a single output unit by a
task-specific linear layer with a 0.1 dropout rate. For the stress detection
task, two output units are used to predict the logits for stress and non-
stress labels. A softmax layer is used to compute probabilities and training
loss. Note that in related work Jaiswal et al. (2020); Bara et al. (2020), the
“late fusion” stands for an ensemble method, while MUSER solves the task with
a single model.

###  4.4 MUSER: Multi-task Learning

####  4.4.1 Weight Sharing

We directly share all the trainable parameters (hard sharing) except the task-
specific output layers. For each batch of training samples, only one task is
assigned, and one step of back-propagation is performed according to the task
objective function with the task-specific output layer plugged in.

####  4.4.2 Sampling Strategy

In each epoch of multi-task training, different amounts of training data are
sampled from both the auxiliary task of activation/valence regression and the
main task of stress classification. We explore both uniform sampling and
dynamic sampling strategy to adaptively decide the mixing ratio of the
multiple tasks in each epoch.

Table 5: Comparison with look-alike methods for multi-task learning.

Method | Type | Based on | History
---|---|---|---
Liu et al. (2019) | Loss Balancing | Ratios of training loss | 1
Gottumukkala et al. (2020) | Dynamic Sampling | Validation metrics | 1
MUSER | Dynamic Sampling | Ratios of validation metrics | $n\geq 5$

Uniform Sampling. In our conditions, the number of training samples in the
main task and the auxiliary tasks are approximately on the same scale.
Therefore, an intuitive method is to switch between the tasks with uniform
sampling: for each batch, we first decide which task to train with an equal
chance, and then randomly select 32 (the batch size) samples; the batch is
trained with the corresponding supervision signals (either emotion scores or
stress labels) from the selected task.

Dynamic Sampling. Having an equal number of samples for each task in each
epoch is not the most efficient way for multi-task training because it is not
searching for the most informative task during each epoch. It is more
intuitive that when one task reaches a bottleneck, more samples from the other
tasks should be selected instead.

Motivated by this idea, we propose to dynamically select the task for each
batch according to the model’s speed of learning each task. After each
training epoch, the sampling distribution is updated based on the model’s
current and historical performance on each task on the validation set.
Specifically, for activation and valence tasks, we compute the ratio of the
average rooted mean square error (RMSE) score on the past $n$ epochs to the
RMSE score in the current epoch. The ratios are noted as $r_{a}$ and $r_{v}$,
respectively. For the stress task, we compute the ratio of the accuracy in the
current epoch to the average of the past $n$ epochs, noted as $r_{s}$. The
history length $n$ is picked by hand. The sampling distribution for the next
epoch is then computed as:

$p_{a},p_{v},p_{s}=softmax([r_{a}/\rho,r_{v}/\rho,r_{s}/\rho]),$ (1)

where $\rho$ is the temperature coefficient set to 0.1.

We use the “ratios to history” instead of the validation scores themselves to
compute the distribution because this makes different metrics comparable to
each other, and it is a good estimation of which task is the model currently
learning the fastest. We name this strategy a “speed-based” dynamic sampling.
The sampling procedure is shown in Figure 2, and a comparison to look-alike
multi-task learning methods is included in Table 5.

Table 6: Results: stress detection with single modality.

Models | Accuracy | Precision | Recall | F-score
---|---|---|---|---
MLP+LIWC Jaiswal et al. (2020) | 0.60 | 0.74 | 0.61 | 0.67
MLP+Opensmile Jaiswal et al. (2020) | 0.67 | 0.70 | 0.69 | 0.69
MUSER Textual Encoder | 0.69 | 0.77 | 0.74 | 0.75
MUSER Acoustic Encoder | 0.79 | 0.80 | 0.79 | 0.80

Table 7: Stress detection results with multiple modalities, with or without
pre-training.

Method | Accuracy | Precision | Recall | F-score
---|---|---|---|---
LIWC+Audio Jaiswal et al. (2020) | 0.60 | 0.74 | 0.61 | 0.67
A-modal GRU Bara et al. (2020) | 0.573 | 0.582 | 0.557 | 0.569
MUSER from scratch | 0.821 | 0.834 | 0.828 | 0.831
Activation-100 | 0.842 | 0.866 | 0.832 | 0.849
Valence-100 | 0.823 | 0.830 | 0.839 | 0.834
Activation-valence-stress | 0.819 | 0.841 | 0.813 | 0.827
Valence-activation-stress | 0.828 | 0.854 | 0.817 | 0.835

Table 8: Stress detection results with multi-task learning

Strategy | Accuracy | Precision | Recall | F-score
---|---|---|---|---
Uniform sampling with activation | 0.832 | 0.833 | 0.857 | 0.845
Uniform sampling with valence | 0.823 | 0.837 | 0.828 | 0.832
Uniform sampling with activation & valence | 0.846 | 0.862 | 0.846 | 0.854
Liu et al. (2019) \- modified | 0.842 | 0.856 | 0.846 | 0.851
Speed-based sampling (history=5) | 0.854 | 0.861 | 0.864 | 0.863
Speed-based sampling (history=10) | 0.856 | 0.867 | 0.861 | 0.864

Table 9: Stress detection results with multi-task learning with OMG-Emotions

Strategy | Accuracy | Precision | Recall | F-score
---|---|---|---|---
Uniform sampling with activation & valence | 0.844 | 0.867 | 0.835 | 0.850
Liu et al. (2019) \- modified | 0.836 | 0.820 | 0.886 | 0.852
Speed-based sampling (history=5) | 0.850 | 0.871 | 0.842 | 0.856

##  5 Experiments

###  5.1 Settings

We use an AdamW Loshchilov and Hutter (2018) optimizer with an initial
learning rate of 3e-4 for all our multimodal and multi-task experiments. In
each epoch, we repeatedly sample data with a batch-size of 32 from the main
task or the auxiliary tasks, and apply one-step back-propagation for each
batch, until the total selected number reaches the size of the MuSE training
set. Gradients are clipped to have a maximum norm of 1.0. The history length
$n$ in speed-based dynamic sampling is chosen from {1, 5, 10} according to the
performance on the validation set. We warm up the dynamic sampling by applying
uniform sampling for the first $n$ epochs. The maximum epoch number is
typically set to 1000, while the training process is controlled by early
stopping. For the Transformer textual encoder, we limit the maximum sequence
length to be 128. The evaluation metrics include overall accuracy, as well as
the precision, recall, and f-score for the “stressed” class.

###  5.2 Unimodal Results

For unimodal experiments, we use the textual encoder or the acoustic encoder
independently to compute representations before regression or classification.
For the Transformer textual encoder, we use a learning rate of 2e-5; for the
MLP acoustic model, we use a learning rate of 5e-4. These learning rates are
separately fine-tuned on each unimodal task. Other hyperparameters of the
models are kept the same as the multimodal structure.

Table 6 shows the stress detection results with single modalities. Our
Transformer encoder outperforms the baseline textual model because of its
capability to discover syntactic-level long distance relationships in natural
language and the external linguistic knowledge from the advanced BERT pre-
training; our acoustic model also improves beyond the baseline results,
potentially because we used a more up-to-date version of eGeMaps configuration
and a fine-tuned learning rate.

###  5.3 Multimodal Results

To jointly train with both the textual and acoustic features, we use the
multimodal fusion model introduced in Section 4.3 as a basic architecture.

####  5.3.1 Pre-training

Our MUSER model is trained from scratch to set up a single-task baseline for
multimodal stress detection. Besides, a potential alternative to multi-task
learning is pre-training on the auxiliary tasks and fine-tuning on the main
task. For a complete comparison, we set up several strategies for pre-
training. All the pre-training methods use the internal auxiliary task of
MuSE. The compared methods are as follows:
Activation-100: pre-train for 100 epochs with the activation annotations, then
switch to the main task of stress detection.
Valence-100: pre-train for 100 epochs with the valence annotations, then
switch to the main task of stress detection.
Activation-valence-stress: pre-train for 100 epochs on the activation task,
then 100 epochs on the valence task, and switch to stress detection.
Valence-activation-stress: pre-train for 100 epochs on the valence task, then
100 epochs on the activation task, and switch to stress detection.

The results are presented in Table 7. Among the pre-training and fine-tuning
results, Activation-100 shows the most significant improvement. The second-
best score is the valence-activation-stress order. Thus, we can conclude that
activation is the better auxiliary task under this paradigm. Additionally,
using only one auxiliary task is always better than using two of them; this is
because when the model learns from the second auxiliary task, it “forgets” the
knowledge from the previous task because it lacks a memory mechanism to look
back Hayes et al. (2020).

Pre-training on the emotion recognition tasks using either activation or
valence improves stress detection because the model is equipped with the
capabilities to encode the features and predict emotions before the training
of stress detection task starts.

####  5.3.2 Multi-task Learning on MuSE

For multi-task learning, we compare two sampling strategies: uniform sampling
and our proposed speed-based dynamic sampling. We also implement and modify
the loss-based weight balancing method proposed by Liu et al. (2019) to adjust
the mixing ratios in dynamic sampling instead, and compare it with our
methods. The results using the internal MuSE emotion recognition as an
auxiliary task are shown in Table 8.

Comparing the uniform sampling results with Table 7, we conclude that using
any auxiliary task is better than training from scratch. However, multi-task
training with the activation and valence tasks together is better than using
them separately. This is different from the observations in Table 7 and can be
explained by the differences in the training procedure: in multi-task
learning, the model looks back-and-forth into each task in each epoch, making
it able to memorize the shared knowledge from all the tasks. Additionally,
when the model is optimized for the two emotion tasks at the same time, the
lower-level representation becomes more general and informative because it is
frequently plugged with different task-specific layers.

Comparing the results of using a single auxiliary task of activation vs.
valence, activation leads to better results as compared to valence, which is
in agreement with Table 7. This is further supported by the analyses in Tables
2 and 4: given the lower unique indicator count of the activation task, as
well as the fact that the pre-training and multi-task learning results are all
compatible, we can conclude that for stress detection, the nature of the
activation dimension of emotion is closer and more helpful than the valence
dimension. This potentially suggests that stress has a major effect on whether
people feel excited (activation), but a minor effect on their opinion toward
events and objects (valence).

We test our speed-based dynamic sampling algorithm using activation and
valence together as auxiliary tasks and it yields promising results with
history set to 5 and 10. It significantly outperforms both the uniform
sampling and our implementation of the loss-based strategy Liu et al. (2019)
(t-test, $p<0.05$), achieving state-of-the-art scores on MuSE stress detection
task with one single model and only two modalities.

Our model works the best with a history length between 5 and 10. If the
history is too short, the model takes longer to converge and has unstable
performance, while if the history is too long, it fails to capture the
dynamics.

####  5.3.3 Generalization

In real-world applications, stress detection data does not necessarily have
emotion labels. However, because of the intrinsic inter-dependence between
emotion and stress, any existing dataset with emotion labels can potentially
serve as an external auxiliary task for stress detection. However, this
requires our model and multi-task training algorithm to generalize beyond the
internal MuSE emotion tasks. We test our model on OMG-Emotions as an example
of external emotion datasets.

Table 9 shows results on MuSE stress detection using OMG-Emotion as an
auxiliary task. Comparing to Table 7, although the source and content of OMG-
Emotions are different from MuSE, multi-task learning still outperforms
single-task learning and pre-training (t-test, $p<0.05$). This reveals that
the connection between stress and emotion widely exists, and our multi-task
learning method works in general cases.

Additionally, Table 9 suggests that while using an external emotion dataset,
our speed-based sampling method still outperforms uniform sampling, as well as
our implementation of loss-based dynamic sampling Liu et al. (2019). This
supports the robustness and effectiveness of our speed-based strategy.

##  6 Conclusions

In this work, we uncovered the connections and differences between stress
detection and emotion recognition using textual and acoustic features, and
proposed to use emotion recognition as an auxiliary task for stress detection.
We proposed MUSER: a Transformer-based model structure, together with a novel
speed-based dynamic sampling strategy for multi-task learning. Experimental
results support the inter-dependence of stress and emotion
(activation/valence), and proves the effectiveness and robustness of our
methods. MUSER achieved state-of-the-art results on the MuSE stress detection
task both when internal (MuSE) and when external (OMG-Emotions) emotion data
and annotations were used.

Our code is publicly available at
https://lit.eecs.umich.edu/downloads.html#MUSER

## Acknowledgements

We would like to thank Cristian-Paul Bara and Mimansa Jaiswal for their
helpful discussion on the data processing and features of MuSE dataset. We
appreciate the insightful comments from all the reviewers and program
committee.

This material is based in part on work supported by the Toyota Research
Institute, the Precision Health initiative at the University of Michigan, the
National Science Foundation (grant #1815291), and the John Templeton
Foundation (grant #61156). Any opinions, findings, conclusions, or
recommendations in this material are those of the authors and do not
necessarily reflect the views of the Toyota Research Institute, Precision
Health initiative, the National Science Foundation, or the John Templeton
Foundation.

## References

  * Adolphs (2003) Ralph Adolphs. 2003.  Is the human amygdala specialized for processing social information?  _Annals of the New York Academy of Sciences_ , 985(1):326–340.
  * Aigrain et al. (2016) Jonathan Aigrain, Michel Spodenkiewicz, Séverine Dubuiss, Marcin Detyniecki, David Cohen, and Mohamed Chetouani. 2016.  Multimodal stress detection from multiple assessments.  _IEEE Transactions on Affective Computing_ , 9(4):491–506.
  * Alberdi et al. (2016) Ane Alberdi, Asier Aztiria, and Adrian Basarab. 2016.  Towards an automatic early stress recognition system for office environments based on multimodal measurements: A review.  _Journal of biomedical informatics_ , 59:49–75.
  * Bara et al. (2020) Cristian Paul Bara, Michalis Papakostas, and Rada Mihalcea. 2020.  A deep learning approach towards multimodal stress detection.  In _Proceedings of the AAAI-20 Workshop on Affective Content Analysis, New York, USA, AAAI_.
  * Barros et al. (2018) Pablo Barros, Nikhil Churamani, Egor Lakomkin, Henrique Siqueira, Alexander Sutherland, and Stefan Wermter. 2018.  The omg-emotion behavior dataset.  In _2018 International Joint Conference on Neural Networks (IJCNN)_ , pages 1–7. IEEE.
  * Botelho et al. (2019) M Catarina Botelho, Isabel Trancoso, Alberto Abad, and Teresa Paiva. 2019.  Speech as a biomarker for obstructive sleep apnea detection.  In _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 5851–5855. IEEE.
  * Bradley and Lang (1994) Margaret M Bradley and Peter J Lang. 1994.  Measuring emotion: the self-assessment manikin and the semantic differential.  _Journal of behavior therapy and experimental psychiatry_ , 25(1):49–59.
  * Busso et al. (2008) Carlos Busso, Murtaza Bulut, Chi-Chun Lee, Abe Kazemzadeh, Emily Mower, Samuel Kim, Jeannette N Chang, Sungbok Lee, and Shrikanth S Narayanan. 2008.  Iemocap: Interactive emotional dyadic motion capture database.  _Language resources and evaluation_ , 42(4):335.
  * Busso et al. (2016) Carlos Busso, Srinivas Parthasarathy, Alec Burmania, Mohammed AbdelWahab, Najmeh Sadoughi, and Emily Mower Provost. 2016.  Msp-improv: An acted corpus of dyadic interactions to study emotion perception.  _IEEE Transactions on Affective Computing_ , 8(1):67–80.
  * Chen et al. (2017) Lan-lan Chen, Yu Zhao, Peng-fei Ye, Jian Zhang, and Jun-zhong Zou. 2017.  Detecting driving stress in physiological signals based on multimodal feature analysis and kernel classifiers.  _Expert Systems with Applications_ , 85:279–291.
  * Chen et al. (2018a) Sheng-Yeh Chen, Chao-Chun Hsu, Chuan-Chun Kuo, Lun-Wei Ku, et al. 2018a.  Emotionlines: An emotion corpus of multi-party conversations.  _arXiv preprint arXiv:1802.08379_.
  * Chen et al. (2018b) Zhao Chen, Vijay Badrinarayanan, Chen-Yu Lee, and Andrew Rabinovich. 2018b.  Gradnorm: Gradient normalization for adaptive loss balancing in deep multitask networks.  In _International Conference on Machine Learning_ , pages 794–803. PMLR.
  * Cohen et al. (1991) Sheldon Cohen, David AJ Tyrrell, and Andrew P Smith. 1991.  Psychological stress and susceptibility to the common cold.  _New England journal of medicine_ , 325(9):606–612.
  * Devlin et al. (2019) Jacob Devlin, Ming-Wei Chang, Kenton Lee, and Kristina Toutanova. 2019.  Bert: Pre-training of deep bidirectional transformers for language understanding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 4171–4186.
  * Dobson and Smith (2000) Hilary Dobson and RF Smith. 2000.  What is stress, and how does it affect reproduction?  _Animal reproduction science_ , 60:743–752.
  * Eyben et al. (2015) Florian Eyben, Klaus R Scherer, Björn W Schuller, Johan Sundberg, Elisabeth André, Carlos Busso, Laurence Y Devillers, Julien Epps, Petri Laukka, Shrikanth S Narayanan, et al. 2015.  The geneva minimalistic acoustic parameter set (gemaps) for voice research and affective computing.  _IEEE transactions on affective computing_ , 7(2):190–202.
  * Eyben et al. (2010) Florian Eyben, Martin Wöllmer, and Björn Schuller. 2010.  Opensmile: the munich versatile and fast open-source audio feature extractor.  In _Proceedings of the 18th ACM international conference on Multimedia_ , pages 1459–1462.
  * Gao et al. (2014) Hua Gao, Anil Yüce, and Jean-Philippe Thiran. 2014.  Detecting emotional stress from facial expressions for driving safety.  In _2014 IEEE International Conference on Image Processing (ICIP)_ , pages 5961–5965. IEEE.
  * Gong et al. (2019) Ting Gong, Tyler Lee, Cory Stephenson, Venkata Renduchintala, Suchismita Padhy, Anthony Ndirango, Gokce Keskin, and Oguz H Elibol. 2019.  A comparison of loss weighting strategies for multi task learning in deep neural networks.  _IEEE Access_ , 7:141627–141632.
  * Gottumukkala et al. (2020) Ananth Gottumukkala, Dheeru Dua, Sameer Singh, and Matt Gardner. 2020.  Dynamic sampling strategies for multi-task reading comprehension.  In _Proceedings of the 58th Annual Meeting of the Association for Computational Linguistics_ , pages 920–924.
  * Guo et al. (2018a) Han Guo, Ramakanth Pasunuru, and Mohit Bansal. 2018a.  Dynamic multi-level multi-task learning for sentence simplification.  In _Proceedings of the 27th International Conference on Computational Linguistics_ , pages 462–476.
  * Guo et al. (2018b) Michelle Guo, Albert Haque, De-An Huang, Serena Yeung, and Li Fei-Fei. 2018b.  Dynamic task prioritization for multitask learning.  In _Proceedings of the European Conference on Computer Vision (ECCV)_ , pages 270–287.
  * Hayes et al. (2020) Tyler L Hayes, Kushal Kafle, Robik Shrestha, Manoj Acharya, and Christopher Kanan. 2020.  Remind your neural network to prevent catastrophic forgetting.  In _European Conference on Computer Vision_ , pages 466–483. Springer.
  * Horvath (1982) Frank Horvath. 1982.  Detecting deception: the promise and the reality of voice stress analysis.  _Journal of Forensic Science_ , 27(2):340–351.
  * Jaiswal et al. (2019) Mimansa Jaiswal, Zakaria Aldeneh, Cristian-Paul Bara, Yuanhang Luo, Mihai Burzo, Rada Mihalcea, and Emily Mower Provost. 2019.  Muse-ing on the impact of utterance ordering on crowdsourced emotion annotations.  In _ICASSP 2019-2019 IEEE International Conference on Acoustics, Speech and Signal Processing (ICASSP)_ , pages 7415–7419. IEEE.
  * Jaiswal et al. (2020) Mimansa Jaiswal, Cristian-Paul Bara, Yuanhang Luo, Mihai Burzo, Rada Mihalcea, and Emily Mower Provost. 2020.  Muse: a multimodal dataset of stressed emotion.  In _Proceedings of The 12th Language Resources and Evaluation Conference_ , pages 1499–1510.
  * Jaques et al. (2016) Natasha Jaques, Sara Taylor, Ehimwenma Nosakhare, Akane Sano, and Rosalind Picard. 2016.  Multi-task learning for predicting health, stress, and happiness.  In _NIPS Workshop on Machine Learning for Healthcare_.
  * Kendall et al. (2018) Alex Kendall, Yarin Gal, and Roberto Cipolla. 2018.  Multi-task learning using uncertainty to weigh losses for scene geometry and semantics.  In _Proceedings of the IEEE conference on computer vision and pattern recognition_ , pages 7482–7491.
  * Lample et al. (2017) Guillaume Lample, Alexis Conneau, Ludovic Denoyer, and Marc’Aurelio Ranzato. 2017\.  Unsupervised machine translation using monolingual corpora only.  _arXiv preprint arXiv:1711.00043_.
  * Lane et al. (2015) Nicholas D Lane, Petko Georgiev, and Lorena Qendro. 2015.  Deepear: robust smartphone audio sensing in unconstrained acoustic environments using deep learning.  In _Proceedings of the 2015 ACM International Joint Conference on Pervasive and Ubiquitous Computing_ , pages 283–294.
  * Lech and He (2014) Margaret Lech and Ling He. 2014.  Stress and emotion recognition using acoustic speech analysis.  In _Mental Health Informatics_ , pages 163–184. Springer.
  * Lin et al. (2014) Huijie Lin, Jia Jia, Quan Guo, Yuanyuan Xue, Qi Li, Jie Huang, Lianhong Cai, and Ling Feng. 2014.  User-level psychological stress detection from social media using deep neural network.  In _Proceedings of the 22nd ACM international conference on Multimedia_ , pages 507–516.
  * Liu et al. (2019) Shikun Liu, Edward Johns, and Andrew J Davison. 2019.  End-to-end multi-task learning with attention.  In _Proceedings of the IEEE Conference on Computer Vision and Pattern Recognition_ , pages 1871–1880.
  * Loshchilov and Hutter (2018) Ilya Loshchilov and Frank Hutter. 2018.  Decoupled weight decay regularization.  In _International Conference on Learning Representations_.
  * McKeown et al. (2011) Gary McKeown, Michel Valstar, Roddy Cowie, Maja Pantic, and Marc Schroder. 2011\.  The semaine database: Annotated multimodal records of emotionally colored conversations between a person and a limited agent.  _IEEE transactions on affective computing_ , 3(1):5–17.
  * Muthukumar and Nachiappan (2013) Kannan Muthukumar and Vasanthi Nachiappan. 2013.  Phosphatidylethanolamine from phosphatidylserine decarboxylase2 is essential for autophagy under cadmium stress in saccharomyces cerevisiae.  _Cell biochemistry and biophysics_ , 67(3):1353–1363.
  * Pasunuru and Bansal (2017) Ramakanth Pasunuru and Mohit Bansal. 2017.  Multi-task video captioning with video and entailment generation.  In _Proceedings of the 55th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 1273–1283.
  * Paulmann et al. (2016) Silke Paulmann, Desire Furnes, Anne Ming Bøkenes, and Philip J Cozzolino. 2016\.  How psychological stress affects emotional prosody.  _Plos one_ , 11(11):e0165022.
  * Pennebaker et al. (2001) James W Pennebaker, Martha E Francis, and Roger J Booth. 2001.  Linguistic inquiry and word count: Liwc 2001.  _Mahway: Lawrence Erlbaum Associates_ , 71(2001):2001.
  * Pruessner et al. (2008) Jens C Pruessner, Katarina Dedovic, Najmeh Khalili-Mahani, Veronika Engert, Marita Pruessner, Claudia Buss, Robert Renwick, Alain Dagher, Michael J Meaney, and Sonia Lupien. 2008.  Deactivation of the limbic system during acute psychosocial stress: evidence from positron emission tomography and functional magnetic resonance imaging studies.  _Biological psychiatry_ , 63(2):234–240.
  * Taylor and Stanton (2007) Shelley E Taylor and Annette L Stanton. 2007.  Coping resources, coping processes, and mental health.  _Annu. Rev. Clin. Psychol._ , 3:377–401.
  * Tull et al. (2007) Matthew T Tull, Heidi M Barrett, Elaine S McMillan, and Lizabeth Roemer. 2007.  A preliminary investigation of the relationship between emotion regulation difficulties and posttraumatic stress symptoms.  _Behavior Therapy_ , 38(3):303–313.
  * Vaswani et al. (2017) Ashish Vaswani, Noam Shazeer, Niki Parmar, Jakob Uszkoreit, Llion Jones, Aidan N Gomez, Łukasz Kaiser, and Illia Polosukhin. 2017.  Attention is all you need.  In _Advances in neural information processing systems_ , pages 5998–6008.
  * Wang and Saudino (2011) Manjie Wang and Kimberly J Saudino. 2011.  Emotion regulation and stress.  _Journal of Adult Development_ , 18(2):95–103.
  * Yao et al. (2020) Yiqun Yao, Verónica Pérez-Rosas, Mohamed Abouelenien, and Mihai Burzo. 2020\.  Morse: Multimodal sentiment analysis for real-life settings.  In _Proceedings of the 2020 International Conference on Multimodal Interaction_ , pages 387–396.
  * Yao et al. (2019) Yiqun Yao, Jiaming Xu, and Bo Xu. 2019.  The world in my mind: Visual dialog with adversarial multi-modal feature encoding.  In _Proceedings of the 2019 Conference of the North American Chapter of the Association for Computational Linguistics: Human Language Technologies, Volume 1 (Long and Short Papers)_ , pages 2588–2598.
  * Zadeh et al. (2016) Amir Zadeh, Rowan Zellers, Eli Pincus, and Louis-Philippe Morency. 2016.  Mosi: multimodal corpus of sentiment intensity and subjectivity analysis in online opinion videos.  _arXiv preprint arXiv:1606.06259_.
  * Zadeh et al. (2018) AmirAli Bagher Zadeh, Paul Pu Liang, Soujanya Poria, Erik Cambria, and Louis-Philippe Morency. 2018.  Multimodal language analysis in the wild: Cmu-mosei dataset and interpretable dynamic fusion graph.  In _Proceedings of the 56th Annual Meeting of the Association for Computational Linguistics (Volume 1: Long Papers)_ , pages 2236–2246.


∎

e1e-mail<EMAIL_ADDRESS>
11institutetext: Niels Bohr Institute, University of Copenhagen, 2200
Copenhagen, Denmark

# Generic multi-particle transverse momentum correlations as a new tool for
studying nuclear structure at the energy frontier

Emil Gorm Dahlbæk Nielsenaddr1  Frederik K. Rømeraddr1  Kristjan
Gulbrandsenaddr1  You Zhoue1,addr1

(Received: date / Accepted: date)

###### Abstract

The mean transverse momentum of produced particles, $[p_{\mathrm{T}}]$, and
its event-by-event fluctuations give direct access to the initial conditions
of ultra-relativistic heavy-ion collisions and help probe the colliding
nuclei’s structure. The $[p_{\mathrm{T}}]$ fluctuations can be studied via
multi-particle $p_{\mathrm{T}}$ correlations; so far, only the lowest four
orders have been studied. Higher-order fluctuations can provide stronger
constraints on the initial conditions and improved sensitivity to the detailed
nuclear structure; however, their direct implementation can be challenging and
is still lacking. In this paper, we apply a generic recursive algorithm for
the genuine multi-particle $p_{\mathrm{T}}$ correlations, which enables the
accurate study of higher-order $[p_{\mathrm{T}}]$ fluctuations without
computationally heavy processing for the first time. With this algorithm, we
will examine the power of multi-particle $p_{\mathrm{T}}$ correlations through
Monte Carlo model studies with different nuclear structures. The impact on the
nuclear structure studies, including the nuclear deformation and triaxial
structure, will be discussed. These results will demonstrate the usefulness of
multi-particle $p_{\mathrm{T}}$ correlations for studying nuclear structure in
high-energy nuclei collisions at RHIC and the LHC, which could serve as
complementary to existing low-energy nuclear structure studies.

###### Keywords:

High-energy nuclear collisions nuclear structure multi-particle correlations

††journal: Eur. Phys. J. A

##  1 Introduction

The primary goal of the heavy-ion collisions at ultra-relativistic energies,
such as at the Relativistic Heavy Ion Collider (RHIC) and the Large Hadron
Collider (LHC), is to create a quark-gluon plasma (QGP) and study its
properties in the laboratory Shuryak:1980tp ; Shuryak:1978ij . Experimental
measurements of collective flow ALICE:2011ab ; ALICE:2016ccg  and its
comparison to hydrodynamic model calculations Muller:2012zq ; Heinz:2013th ;
Song:2017wtw , in particular via global Bayesian analysis Bernhard:2019bmu ;
JETSCAPE:2020mzn , provide a powerful approach to extracting unique
information about the specific viscosity of the QGP and its time evolution
(temperature dependence). This represents the state-of-the-art understanding
of the QGP ALICE:2022wpn . However, the sizable uncertainty of the extracted
QGP transport coefficients originates from poor knowledge about the initial
conditions of heavy-ion collisions. A further reduction of this uncertainty
can be achieved via in-depth investigations of the initial conditions
Bernhard:2019bmu ; JETSCAPE:2020mzn . Among various collective flow
observables, the mean transverse momentum, $[p_{\mathrm{T}}]$, and its event-
by-event fluctuations trace their origins to the initial state of the heavy-
ion collision Broniowski:2009fm ; Bozek:2012fw  and are arguably a more direct
way to study the initial state fluctuations Giacalone:2020lbm . Several
initial state properties have been proposed as good estimators of event-by-
event mean transverse momentum fluctuations. These properties include the
initial size $R$ and the size fluctuations of the system that were introduced
in Broniowski:2009fm , the transverse area of the overlap region and the
initial eccentricity $\epsilon_{2}$, introduced as $A_{e}$ in Schenke:2020uqq
, as well as the initial energy per unit rapidity at the initial time
$\tau_{0}$, $E_{i}$ Gardim:2020sma . These estimators have been verified using
hydrodynamic calculations that show $E_{i}$ provides a strictly linear
correlation to the $[p_{\mathrm{T}}]$ Gardim:2019xjs . Essentially, the
$[p_{\mathrm{T}}]$ and its fluctuations reveal information about the
fluctuations of the energy deposited in the initial state on an event-by-event
basis.

Figure 1: Illustration of the surface of nuclei with different quadrupole
deformation. From left to right: spherical, prolate, triaxial, and oblate.

It has been proposed that an improved understanding of the initial conditions
of nuclei collisions will allow for a more precise description of the
geometrical shape of the overlapping colliding region Bally:2022vgo . In
particular, for central collisions where the two colliding nuclei fully
overlap, the spatial shape of the overlapping region could directly reflect
the nuclear structure that existed before the collision occurred (in the case
of body-body collisions). It has also been shown that the imprint of the
nuclear structure remains in the subsequent final state of the collective
expansion, i.e., $[p_{\mathrm{T}}]$ fluctuations Bhatta:2021qfk  and
anisotropic flow Lu:2023fqd , after the dynamic evolution of nucleus
collisions. The study of nuclear structure is an important research topic in
nuclear physics, as it reveals the fundamental properties and interactions of
nucleons and nuclei. These studies allow us to test and improve the
theoretical models and methods that describe the nuclear many-body problem,
which remains a critical challenge from the low-energy side Bally:2022vgo .
The imaging power of high-energy nucleus collisions provides a complementary
approach to the existing low-energy nuclear structure study Giacalone:2023hwk
.

In the initial conditions of heavy-ion collisions, the nucleon positions are
sampled from a modified Woods-Saxon distribution

$\displaystyle\rho^{\mathrm{WS}}(r)$
$\displaystyle=\frac{\rho_{0}}{1+\exp\left(\frac{r-R(\theta,\varphi)}{a_{0}}\right)}.$
(1)

The deformation of the nuclear surface is described by

$\displaystyle R(\theta,\varphi)$
$\displaystyle=R_{0}\left(1+\beta_{2}\left[\cos\gamma
Y_{20}(\theta,\varphi)+\sin\gamma Y_{22}(\theta,\varphi)\right]\right)$ (2)

where $R_{0}$ is the nuclear radius, $Y_{lm}$ are spherical harmonics,
$\beta_{2}$ quantifies the strength of the quadrupole deformation, and the
triaxial parameter $\gamma$ controls the relative order of the three axes in
the intrinsic frame of the nuclei. The impacts of $\beta_{2}$ and $\gamma$ on
the nuclear structure are illustrated in Fig. 1. Systematic studies of
anisotropic flow using multi-particle azimuthal angle correlations show that
several flow observables are sensitive to the nuclear quadrupole deformation
$\beta_{2}$ Lu:2023fqd . However, none of the flow observables (in terms of
azimuthal angle correlations) have the power to distinguish the triaxial
structures Lu:2023fqd , suggesting that different choices of observable may be
needed to accomplish this, such as $v_{n}^{2}$-$[p_{\mathrm{T}}]$ correlations
Giacalone:2019pca ; Bally:2021qys ; ALICE:2021gxt ; Bally:2023dxi , or multi-
particle $p_{\mathrm{T}}$ correlations Broniowski:2009fm ; Giacalone:2020lbm .

This paper will present a generic recursive algorithm for the genuine multi-
particle correlations (cumulants) of $p_{\mathrm{T}}$ in Section 2. The Monte
Carlo models, introduced in Section 3, will be used to examine the power of
multi-particle cumulants of $p_{\mathrm{T}}$, with the results presented in
Section 4. Finally, we will summarize in Section 5.

##  2 Algorithm for multi-particle $p_{\mathrm{T}}$ correlations

Multi-particle correlations have been a powerful tool in the study of heavy-
ion collisions Borghini:2000sa ; Ollitrault:2023wjk . One successful example
is using multi-particle azimuthal angle correlations to study the final state
anisotropic particle expansion, quantified by anisotropic flow $v_{n}$
Bilandzic:2013kga ; Moravcova:2020wnf . Systematic studies of anisotropic
flow, flow fluctuations, and flow correlations show that the selected flow
observables are sensitive to the nuclear quadrupole deformation $\beta_{2}$
Lu:2023fqd . Still, none of the flow observables has the power to distinguish
the triaxial structures Lu:2023fqd . For the study of multi-particle
$p_{\mathrm{T}}$ correlations, the lower orders have been available and used
in the investigations of the initial conditions in heavy-ion collisions and
also explored for their potential in the study of nuclear structure
Jia:2021qyu ; Bhatta:2021qfk . Meanwhile, it is challenging and
computationally demanding for the higher orders when one manually derives the
correlations. A generic recursive algorithm provides an ideal solution. It has
been successfully developed for multi-particle azimuthal correlations
Moravcova:2020wnf  and used in the above-mentioned anisotropic flow studies
Lu:2023fqd . Following a similar idea, we propose a generic recursive
algorithm for multi-particle $[p_{\mathrm{T}}]$ correlations in this paper.

In a single collision, many particles are produced. The averaged
$p_{\mathrm{T}}$ correlation in this event is defined as

$\displaystyle[p_{\mathrm{T}}^{(m)}]=\frac{\displaystyle\sum^{M}_{k_{1}\neq\ldots\neq
k_{m}}w_{k_{1}}\cdot\ldots\cdot w_{k_{m}}p_{\mathrm{T},k_{1}}\cdot\ldots\cdot
p_{\mathrm{T},k_{m}}}{\displaystyle\sum^{M}_{k_{1}\neq\ldots\neq
k_{m}}w_{k_{1}}\cdot\ldots\cdot w_{k_{m}}},$ (3)

where $M$ is the number of particles within a certain kinematic range in this
event and $w_{k_{i}}$ is a particle weight correcting for detector
inefficiencies. Applying the constraint that $k_{1}\neq\ldots\neq k_{m}$
removes any auto-correlations such that the averaged $p_{\mathrm{T}}$
correlation measures the dynamic correlations due to event-by-event
fluctuations and not the statistical correlations due to the finite number of
tracks. The notation $[p_{\mathrm{T}}]$ is used for the mean transverse
momentum calculated event-by-event within a certain kinematic range and is
distinct from the mean transverse momentum, $\langle p_{\mathrm{T}}\rangle$,
extrapolated to $p_{\mathrm{T}}=0$ from spectra analyses over a large event
ensemble ALICE:2010syw . Similar to the Generic Framework for azimuthal angle
correlations Bilandzic:2013kga ; Moravcova:2020wnf , the numerator and
denominator can be defined as

$\displaystyle N\langle m\rangle_{p_{\mathrm{T}}}$
$\displaystyle=\displaystyle\sum^{M}_{k_{1}\neq\ldots\neq
k_{m}}w_{k_{1}}\cdot\ldots\cdot w_{k_{m}}p_{\mathrm{T},k_{1}}\cdot\ldots\cdot
p_{\mathrm{T},k_{m}},$ (4) $\displaystyle D\langle m\rangle_{p_{\mathrm{T}}}$
$\displaystyle=\displaystyle\sum^{M}_{k_{1}\neq\ldots\neq
k_{m}}w_{k_{1}}\cdot\ldots\cdot w_{k_{m}},$ (5)

so that the averaged $m$-particle $p_{\mathrm{T}}$ correlation can be written
as

$\displaystyle[p_{\mathrm{T}}^{(m)}]=\frac{N\langle
m\rangle_{p_{\mathrm{T}}}}{D\langle m\rangle_{p_{\mathrm{T}}}}.$ (6)

The numerator and denominator in Eqs. (4) and (5), respectively, are
calculated by constructing sums of different powers of the transverse momentum
and particle weights

$\displaystyle P_{k}=\sum^{M}_{i=1}w_{i}^{k}p_{\mathrm{T},i}^{k}\quad,\quad
W_{k}=\sum^{M}_{i=1}w_{i}^{k}.$ (7)

Unlike the previously done analyses of azimuthal angle correlations
Bilandzic:2013kga ; Moravcova:2020wnf , the first-order $p_{\mathrm{T}}$
correlation does not vanish but is simply $[p_{\mathrm{T}}]$. The first few
transverse momentum correlations are given by

$\displaystyle[p_{\mathrm{T}}^{(1)}]$
$\displaystyle=\frac{P_{1}}{W_{1}}=\frac{\sum_{i=1}^{M}w_{i}p_{\mathrm{T},i}}{\sum_{i=1}^{M}w_{i}}=[p_{\mathrm{T}}],$
(8) $\displaystyle[p_{\mathrm{T}}^{(2)}]$
$\displaystyle=\frac{P_{1}^{2}-P_{2}}{W_{1}^{2}-W_{2}},$ (9)
$\displaystyle[p_{\mathrm{T}}^{(3)}]$
$\displaystyle=\frac{P_{1}^{3}-3P_{2}P_{1}+2P_{3}}{W_{1}^{3}-3W_{2}W_{1}+2W_{3}}$
(10) $\displaystyle[p_{\mathrm{T}}^{(4)}]$
$\displaystyle=\frac{P_{1}^{4}-6P_{2}P_{1}^{2}+3P_{2}^{2}+8P_{3}P_{1}-6P_{4}}{W_{1}^{4}-6W_{2}W_{1}^{2}+3W_{2}^{2}+8W_{3}W_{1}-6W_{4}}$
(11)

A general recursive formula can be constructed for Eqs. (4) and (5), allowing
one to calculate a $p_{\mathrm{T}}$ correlation of any order from powers of
$P_{k}$ and $W_{k}$

$\displaystyle N\langle
m\rangle_{p_{\mathrm{T}}}=\sum_{k=1}^{m}(-1)^{k-1}N\langle
m-k\rangle_{p_{\mathrm{T}}}\frac{(m-1)!}{(m-k)!}P_{k}$ (12) $\displaystyle
D\langle m\rangle_{p_{\mathrm{T}}}=\sum_{k=1}^{m}(-1)^{k-1}D\langle
m-k\rangle_{p_{\mathrm{T}}}\frac{(m-1)!}{(m-k)!}W_{k}$ (13)

where $N\langle 0\rangle_{p_{\mathrm{T}}}\equiv 1$ and $D\langle
0\rangle_{p_{\mathrm{T}}}\equiv 1$. Further written out applications of this
recursive formula up to eighth order can be found in A.

The event-averaged $p_{\mathrm{T}}$ correlations correspond to the raw sample
moments of the $p_{\mathrm{T}}$ distribution within some kinematic region

$\displaystyle\langle[p_{\mathrm{T}}^{(m)}]\rangle=\frac{\displaystyle\sum_{\mathrm{events}}W^{\prime}_{m}[p_{\mathrm{T}}^{(m)}]}{\displaystyle\sum_{\mathrm{events}}W^{\prime}_{m}},$
(14)

where $W^{\prime}_{m}$ is an event weight. The choice of this weight here is
$W^{\prime}_{m}=D\langle m\rangle_{p_{\mathrm{T}}}$, which is the number of
particle pairs similar to the weight used in anisotropic flow calculations to
minimize the effects of multiplicity fluctuations. The higher-order moments
contain contributions from the lower orders. The $p_{\mathrm{T}}$ cumulants,
$\kappa_{m}$, (calculated from the correlations,
$\langle[p_{\mathrm{T}}^{(m)}]\rangle$) reveal to which degree there exists a
genuine $p_{\mathrm{T}}$ correlation amongst many particles (m) which cannot
be factorized into correlations amongst fewer particles. A recursive formula
for the $m$th-order cumulant Smith:2012cum  is given by

$\displaystyle\kappa_{m}=\langle[p_{\mathrm{T}}^{(m)}]\rangle-\sum_{k=1}^{m-1}\begin{pmatrix}m-1\\\
k-1\end{pmatrix}\langle[p_{\mathrm{T}}^{(m-k)}]\rangle\kappa_{k}.$ (15)

One can expand this formula by grouping common terms to get a formula solely
dependent on powers of $\langle[p_{\mathrm{T}}^{(m)}]\rangle$. This has been
done in B up to eighth order. These genuine correlations will be shown to be
sensitive to the nuclear structure of the colliding nuclei.

##  3 Introduction to the theoretical models

The Heavy Ion Jet INteraction Generator, HIJING Wang:2000bf , is a Monte Carlo
event generator for parton and particle production in high-energy nuclear
collisions. Based on QCD-inspired models for multiple jet production, it
incorporates mechanisms such as multiple mini-jet production, soft excitation,
nuclear shadowing of parton distribution functions, and jet interactions in
dense hadronic matter. It has been reported that with the proper tuning, the
HIJING model can describe particle production in high-energy proton-proton,
proton-nucleus, and nucleus-nucleus collisions Wang:2000bf . Thus, the HIJING
model is expected to provide reliable baseline predictions for the mean
transverse momentum and its event-by-event fluctuations for heavy-ion
collisions at the LHC.

Besides the HIJING model, A Multi-Phase Transport model, AMPT, is another
popular Monte Carlo model that simulates nuclei collisions at relativistic
energies Lin:2004en . The AMPT model run with String Melting consists of four
main parts. First, the initial condition is based on the HIJING model
Wang:2000bf , which describes how the nucleon distributions are arranged and
moved inside the nuclei before they collide. Second, the Parton Cascade (ZPC)
Zhang:1997ej  simulates the partonic interactions inside the nuclei during the
collision. Third, the hadronization, which is done via quark-coalescence
Chen:2005mr , converts the partons into hadrons. Last, the hadronic
interactions are modeled using the ART model Li:1995pra . The AMPT model can
reproduce many experimental measurements from heavy-ion collisions, such as
the production of various particle types and their collective behavior (the
anisotropic flow). The AMPT model with default settings has difficulty
reproducing the centrality dependence of $[p_{\mathrm{T}}]$ measurements in
heavy-ion collisions at the LHC. Nevertheless, one can implement the nuclear
structures in the initial state of the AMPT model. This allows one to examine
if the imprint of the nuclear structure can be observed in the final state
observables despite undergoing both partonic and hadronic interactions.

Data set # | $\beta_{2}$ | $\gamma$ | $a_{0}$ | # of events
---|---|---|---|---
1 | 0 | 0 | 0.57 |  $\sim$1M
2 | 0.18 | 0 | 0.57 |  $\sim$1.5M
3 | 0.18 | 27∘ | 0.57 |  $\sim$2.5M
4 | 0.18 | 60∘ | 0.57 |  $\sim$2M
Table 1: Simulated Xe–Xe collisions data sets in AMPT with four different
configurations of the 129Xe nuclear structure.

This paper implements four different settings for the 129Xe structure. The
parameters relevant to the nuclear structure studies are listed in Table 1,
while other input parameters in the AMPT model are kept fixed between the data
sets. Like many previous AMPT model studies Lu:2023fqd , the Lund string
parameters are set to $a=0.3$ and $b=0.15$, the screening mass to
$\mu=2.2814$. The time-step is 0.2 fm, and the number of time steps for the
hadronic interactions within ART is set to 150, which gives a hadron cascade
time of 30 fm/$c$.

Final state |  Initial state |  Liquid-drop
---|---|---
cumulant |  cumulant |  model
$\kappa_{2}$ |  $\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle$ |  $\frac{1}{32\pi}\langle\beta_{2}^{2}\rangle$
$\kappa_{3}$ |  $\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{3}\right\rangle$ |  $\frac{\sqrt{5}}{896\pi^{3/2}}\langle\cos(3\gamma)\beta_{2}^{3}\rangle$
$\kappa_{4}$ |  $\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{4}\right\rangle-3\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle^{2}$ |  $-\frac{3}{14336\pi^{2}}(7\langle\beta_{2}^{2}\rangle-5\langle\beta_{2}^{4}\rangle)$
$\kappa_{5}$ |  $\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{5}\right\rangle-10\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{3}\right\rangle\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle$ |  $-\frac{5\sqrt{5}}{315392\pi^{5/2}}(11\langle\cos(3\gamma)\beta_{2}^{3}\rangle\langle\beta_{2}^{2}\rangle-5\langle\beta_{2}^{5}\rangle)$
$\kappa_{6}$ |  $\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{6}\right\rangle-15\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{4}\right\rangle\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle$ |  $\frac{5}{918412504\pi^{3}}(42042\langle\beta_{2}^{2}\rangle^{3}-5720\langle\cos(3\gamma)\beta_{2}^{3}\rangle^{2}$
$+30\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle^{3}-10\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{3}\right\rangle^{2}$ |  $-45045\langle\beta_{2}^{2}\rangle\langle\beta_{2}^{4}\rangle+8575\langle\beta_{2}^{6}\rangle+700\langle\cos(6\gamma)\beta_{2}^{6}\rangle)$
$\kappa_{7}$ |  $\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{7}\right\rangle-21\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{5}\right\rangle\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle$ |  $-\frac{15\sqrt{5}}{524812288}(2002\langle\beta_{2}^{2}\rangle^{2}\langle\cos(3\gamma)\beta_{2}^{3}\rangle$
$+210\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{3}\right\rangle\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle^{2}$ |  $+715\langle\cos(3\gamma)\beta_{2}^{3}\rangle\langle\beta_{2}^{4}\rangle$
$-35\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{3}\right\rangle\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{4}\right\rangle$ |  $+910\langle\cos(3\gamma)\beta_{2}^{5}\rangle\langle\beta_{2}^{2}\rangle-175\cos(3\gamma)\beta_{2}^{7}\rangle)$
$\kappa_{8}$ |  $\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{8}\right\rangle-28\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{6}\right\rangle\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle$ |  $\frac{5}{142748942336\pi^{4}}(2144142\langle\beta_{2}^{2}\rangle^{4}-3063060\langle\beta_{2}^{2}\rangle^{2}\langle\beta_{2}^{4}\rangle$
$+420\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{4}\right\rangle\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle^{2}$ |  $-340\langle\beta_{2}^{2}\rangle\biggl{(}2288\langle\cos(3\gamma)\beta_{2}^{3}\rangle^{2}-35\bigl{(}49\langle\beta_{2}^{6}\rangle$
$-35\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{4}\right\rangle^{2}-630\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle^{4}$ |  $+4\langle\cos(6\gamma)\beta_{2}^{6}\rangle\bigr{)}\biggr{)}+25\biggl{(}21879\langle\beta_{2}^{4}\rangle^{2}$
$+560\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{3}\right\rangle^{2}\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{2}\right\rangle$ |  $+14144\langle\cos(3\gamma)\beta_{2}^{3}\rangle\langle\cos(3\gamma)\beta_{2}^{5}\rangle$
$-56\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{5}\right\rangle\cdot\left\langle\left(\frac{\delta d_{\perp}}{d_{\perp}}\right)^{3}\right\rangle$ |  $-35\bigl{(}79\langle\beta_{2}^{8}\rangle+16\langle\cos(6\gamma)\beta_{2}^{8}\rangle\bigr{)}\biggr{)}$
Table 2: The cumulants of $d_{\perp}$ up to eighth order in a liquid-drop
model potential averaged over random orientations. The first three entries are
given in Jia:2021qyu .

Following the idea established in Jia:2021qyu , we consider the nuclei in a
liquid-drop model with a sharp surface potential and in head-on collisions
with a zero impact parameter. The mean transverse momentum is positively
correlated with the inverse transverse size, $d_{\perp}$, in the overlapping
region between the colliding nuclei

$\displaystyle d_{\perp}=\sqrt{N_{\mathrm{part}}/\langle
r_{\perp}^{2}\rangle},$ (16)

where $N_{\mathrm{part}}$ is the number of nucleons participating in the
collision and $r_{\perp}$ is the transverse radius. Within the liquid-drop
model, the dependence of the central moment of $d_{\perp}$ on the nuclear
structure parameters can be roughly approximated with the following relation
Jia:2021qyu

$\displaystyle\frac{\delta d_{\perp}}{d_{\perp}}$
$\displaystyle=\sqrt{\frac{5}{16\pi}}\beta_{2}\biggl{(}\cos(\gamma)D_{0,0}^{2}(\Omega)$
(17)
$\displaystyle+\frac{\sin(\gamma)}{\sqrt{2}}[D_{0,2}^{2}(\Omega)+D_{0,-2}^{2}(\Omega)]\biggr{)},$

where $D_{m,m^{\prime}}^{2}$ are the elements of the Wigner D-matrix and
$\Omega$ represents the Euler angles. The expressions for the higher-order
central moments are obtained by integrating over the Euler angles

$\displaystyle\biggl{\langle}\biggl{(}\frac{\delta
d_{\perp}}{d_{\perp}}\biggr{)}^{n}\biggr{\rangle}$
$\displaystyle=\beta_{2}^{n}\biggl{(}\frac{5}{16\pi}\biggr{)}^{n/2}\int\biggl{(}\cos(\gamma)D_{0,0}^{2}$
(18)
$\displaystyle+\frac{\sin(\gamma)}{\sqrt{2}}[D_{0,2}^{2}+D_{0,-2}^{2}]\biggr{)}^{n}\frac{d\Omega}{8\pi^{2}}$

from which the cumulants can be constructed, and the final cumulants of $n$-th
order are scaled by $1/2^{n-1}$ to account for the independent orientations of
the two nuclei. Table 2 shows the liquid-drop model estimates for the higher-
order cumulants.

##  4 Results

Figure 2: The mean transverse momentum, $\langle[p_{\mathrm{T}}]\rangle$, in
Xe–Xe collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.44 TeV as a function of
centrality for four configurations of the nuclear structure simulated with the
AMPT model.  Figure 3: The second-order mean transverse momentum cumulant,
$\kappa_{2}$, in Xe–Xe collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.44 TeV as a
function of centrality for four configurations of the nuclear structure
simulated with the AMPT model. Figure 4: The third-order mean transverse
momentum cumulant, $\kappa_{3}$, in Xe–Xe collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 5.44 TeV as a function of centrality for four
configurations of the nuclear structure simulated with the AMPT model. Figure
5: The fourth-order mean transverse momentum cumulant, $\kappa_{4}$, in Xe–Xe
collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.44 TeV as a function of centrality
for four configurations of the nuclear structure simulated with the AMPT
model.

Figure 2 shows the $\langle[p_{\mathrm{T}}]\rangle$ in Xe–Xe collisions at
$\sqrt{s_{\mathrm{NN}}}$ = 5.44 TeV as a function of centrality. The top panel
shows the $\langle[p_{\mathrm{T}}]\rangle$ for a spherical, prolate, triaxial,
and oblate configuration of 129Xe, and the bottom panel shows the ratio of
$\langle[p_{\mathrm{T}}]\rangle$ in collisions of deformed nuclei against the
$\langle[p_{\mathrm{T}}]\rangle$ calculated in collisions of spherical nuclei.
The different ratios are consistent with unity suggesting that to first order,
the transverse momentum cumulants are not a useful probe of the nuclear
structure as the $\langle[p_{\mathrm{T}}]\rangle$ does not differ between the
spherical ($\beta_{2}=0$), prolate ($\beta_{2}=0.18$), triaxial
($\beta_{2}=0.18$, $\gamma=27^{\circ}$), and oblate ($\beta_{2}=0.18$,
$\gamma=60^{\circ}$) nuclei. This is consistent with the expectation from the
liquid-drop model; the integral in (18) yields zero when evaluated for $n=1$.

The centrality dependence of $[p_{\mathrm{T}}]$ fluctuations, in the form of
the second-order cumulant or variance, $\kappa_{2}$, is shown in Fig. 3. The
$\kappa_{2}$ result increases by about 10-15% in the central collisions for
$\beta_{2}\neq 0$ compared to the $\kappa_{2}$ for spherical nuclei. However,
the results show no distinction between $\kappa_{2}$ calculations from
prolate, triaxial, and oblate nuclei. This agrees with the expectation from
the liquid-drop model, which predicts a positive contribution (an increase) in
the second-order cumulant that is proportional to
$\langle\beta_{2}^{2}\rangle$ and no sensitivity to the $\gamma$ parameter,
shown in Table 2. Therefore, we recommend that $\kappa_{2}$ should be
considered as an additional probe to the anisotropic flow measurements
ALICE:2018lao  for constraining the parameter of 129Xe. In particular, it
should be included in the Bayesian fits on the experimental measurements to
extract reliable deformation parameters from the high-energy heavy-ion
collisions.

The centrality dependence of the third-order cumulant, $\kappa_{3}$, is shown
in Figure 4, where three distinct trends in central collisions, depending on
the triaxiality parameter $\gamma$, are observed. Compared to the calculations
from spherical nuclei, the magnitude of $\kappa_{3}$ is around 30% larger for
oblate nuclei ($\gamma=60^{\circ}$) and around 30% smaller for prolate nuclei
($\gamma=0^{\circ}$) in central collisions albeit with large statistical
uncertainties. The triaxial 129Xe ($\gamma=27^{\circ}$) leads to a value of
$\kappa_{3}$ consistent with the spherical case ($\beta_{2}=0$) in the 10%
most central collisions. The emergence of the sensitivity to the triaxiality
in $\kappa_{3}$ is unsurprising as triaxiality is a three-point structure, so
a three-particle correlation is needed to measure it. This also explains why
$\kappa_{2}$ can distinguish the overall quadrupole deformation
($\beta_{2}=0.18$) from the spherical case but offers no additional
information: the two-particle $\kappa_{2}$ can only probe two axes of the
nuclei at a time. In addition, the liquid-drop model predicts a positive
contribution from the $\gamma$ parameter in the prolate case
($\cos(3\gamma)=1$), zero contribution from fully triaxial
($\cos(3\gamma)=0$), and negative contribution from the oblate case
($\cos(3\gamma)=-1$). The AMPT calculations follow the same trend. However, as
the $\kappa_{3}$ values are negative, the ordering of the points with
different $\gamma$ parameters is inverted in the ratio plot, shown in the
bottom panel of Fig. 4. The above results further underline the importance of
$\kappa_{3}$ in exploring triaxial nuclear structure. It suggests that not
only the correlations between $p_{\mathrm{T}}$ and $v_{2}^{2}$, which have
been widely used Bozek:2016yoj ; Giacalone:2020dln ; ALICE:2021gxt  but also
the third order $[p_{\mathrm{T}}]$ fluctuations in central collisions should
be considered in the comparisons between model calculations and experimental
measurements from the LHC experiments.

Figure 5 shows the centrality dependence of the fourth-order cumulant,
$\kappa_{4}$, and the comparison between the different nuclear structures.
Similar to $\kappa_{3}$, $\kappa_{4}$ results can distinguish between
different triaxial shapes in central collisions despite sizable uncertainties.
The collisions of oblate nuclei yield the largest $\kappa_{4}$ in central
collisions. The $\kappa_{4}$ for the triaxial case is smaller than the oblate
case but larger than for spherical nuclei. The prolate nuclei yield the
smallest $\kappa_{4}$ in central collisions, almost a 100% difference with
respect to the spherical case. This does not match the prediction from the
simple liquid-drop model, where the fourth-order cumulant is expected to be
proportional to some combination of the moments of $\beta_{2}$ only. However,
a dependence of the initial state cumulant $(\delta
d_{\perp}/d_{\perp})^{4}-3(\delta d_{\perp}/d_{\perp})^{2}$ on the $\gamma$
parameter is also seen in MC-Glauber simulations of various deformed collision
systems Jia:2021qyu . The ordering of the calculations in the bottom plot of
Fig. 5 is similar to that in Fig. 4 (bottom), but in this case, the ordering
is not inverted, suggesting that $\gamma$ either has an opposite contribution
to $\kappa_{4}$ compared to $\kappa_{3}$ or that the $\gamma$ affects the
magnitude of the cumulants. That $\kappa_{4}$ is sensitive to the $\gamma$
parameter is surprising based on the leading-order approximation of the
liquid-drop model, and it would be interesting to see whether $\kappa_{4}$
exhibits similar sensitivity within hydrodynamic models.

Figure 6: The cumulants of the mean transverse momentum, $\kappa_{n}$ ($n\leq
8$), in Xe–Xe collisions at $\sqrt{s_{\mathrm{NN}}}$ = 5.44 TeV as a function
of centrality simulated with the HIJING model.

A baseline prediction for the various orders of $\kappa_{n}$ $(n\leq 8)$ is
given with the HIJING model in Fig. 6. Within the HIJING model, the
fluctuations of the transverse momentum follow that of a superposition of
independent sources, which can be both short- and long-range in nature
Bhatta:2021qfk ; Wang:1991hta . The first three orders of cumulants,
$\kappa_{1}$-$\kappa_{3}$, are predicted to be strictly positive across the
centrality range. $\kappa_{4}$ is positive in central collisions, changes
signs in semi-central collisions and then becomes positive again above 30%
centrality. The higher-order cumulants ($n>4$) fluctuate around zero in
central and semi-central collisions, except $\kappa_{6}$, which is strictly
negative. The non-zero higher-order cumulants suggest that the fluctuations of
the mean transverse momentum are not trivial even in heavy-ion collisions
viewed as a superposition of pp collisions.

##  5 Summary

Multi-particle correlations and cumulants of transverse momentum constrain the
size and size fluctuations in the initial conditions of heavy-ion collisions,
making them an ideal probe of the nuclear structure of colliding nuclei at
ultra-relativistic colliders. This paper proposes a new generic algorithm of
multi-particle $p_{\mathrm{T}}$ correlations that enables an efficient and
precise study of $p_{\mathrm{T}}$ fluctuations up to arbitrary orders. Using
the AMPT transport model with various settings of nuclear structure, including
spherical, prolate, triaxial, and oblate shapes of 129Xe, we investigate the
impact of nuclear structure on the first four orders of $[p_{\mathrm{T}}]$
fluctuations. We find that at first order, the $p_{\mathrm{T}}$ cumulants are
insensitive to the nuclear structure and that at least second order is
required to probe the deformation of the colliding nuclei. The variance,
$\kappa_{2}$, shows a $\beta_{2}$-dependent increase in central collisions but
cannot distinguish between the prolate, triaxial, or oblate cases consistent
with the expectation from the simple liquid-drop model. The third- and fourth-
order cumulants, $\kappa_{3}$ and $\kappa_{4}$, show a $\gamma$-dependent
splitting of the calculated values in the most central collisions with
opposite contribution. The $\gamma$-dependence of $\kappa_{4}$ is not expected
from the liquid-drop model prediction, suggesting that a more complex model of
the interplay between nuclear structure and size fluctuations is necessary to
understand this effect.

These studies show the impact of nuclear structure on the $[p_{\mathrm{T}}]$
fluctuations and how they can be utilized to probe the shape of the colliding
nuclei at the energy frontier (at the level of TeV). The $p_{\mathrm{T}}$
cumulants are, through their relation to the size fluctuations, coupled to the
moments of the deformation parameters. Further study of these relations and
measurements of the cumulants at experiments such as those at the LHC can
introduce a new avenue for studying nuclear structure in high-energy heavy-ion
collisions, which is complementary to existing low-energy measurements and can
serve as experimental validation of low-energy theoretical predictions.
Furthermore, including the $p_{\mathrm{T}}$ cumulants as experimental inputs
in Bayesian analyses can provide independent constraints on the deformation
parameters.

The study of nuclear structure with heavy-ion collisions is an evolving field,
and the presented method and observables can provide a useful approach in
potential future runs of deformed nuclei species at the LHC to constrain the
shape of nuclei across the nuclide chart.

###### Acknowledgements.

We thank Chunjian Zhang for his help on the AMPT code and Anna Ingmer Boye for
her initial contributions to the work. The authors are funded by the European
Union (ERC, InitialConditions), VILLUM FONDEN (grant number 00025462), and
Danmarks Frie Forskningsfond (Independent Research Fund Denmark).

Data availability: The manuscript has associated data in a data repository
[Authors’ comment: Data are available upon request.]

## References

  * (1) E. V. Shuryak, “Quantum Chromodynamics and the Theory of Superdense Matter”, Phys. Rept. 61 (1980) 71–158.
  * (2) E. V. Shuryak, “Quark-Gluon Plasma and Hadronic Production of Leptons, Photons and Psions”, Phys. Lett. 78B (1978) 150. [Yad. Fiz.28,796(1978)].
  * (3) ALICE Collaboration, K. Aamodt et al., “Higher harmonic anisotropic flow measurements of charged particles in Pb–Pb collisions at $\sqrt{s_{\mathrm{NN}}}$ = 2.76 TeV”, Phys. Rev. Lett. 107 (2011) 032301, arXiv:1105.3865 [nucl-ex].
  * (4) ALICE Collaboration, J. Adam et al., “Anisotropic flow of charged particles in Pb-Pb collisions at $\sqrt{s_{\rm NN}}=5.02$ TeV”, Phys. Rev. Lett. 116 no. 13, (2016) 132302, arXiv:1602.01119 [nucl-ex].
  * (5) B. Muller, J. Schukraft, and B. Wyslouch, “First Results from Pb+Pb collisions at the LHC”, Ann. Rev. Nucl. Part. Sci. 62 (2012) 361–386, arXiv:1202.3233 [hep-ex].
  * (6) U. Heinz and R. Snellings, “Collective flow and viscosity in relativistic heavy-ion collisions”, Ann. Rev. Nucl. Part. Sci. 63 (2013) 123–151, arXiv:1301.2826 [nucl-th].
  * (7) H. Song, Y. Zhou, and K. Gajdosova, “Collective flow and hydrodynamics in large and small systems at the LHC”, Nucl. Sci. Tech. 28 no. 7, (2017) 99, arXiv:1703.00670 [nucl-th].
  * (8) J. E. Bernhard, J. S. Moreland, and S. A. Bass, “Bayesian estimation of the specific shear and bulk viscosity of quark–gluon plasma”, Nature Phys. 15 no. 11, (2019) 1113–1117.
  * (9) JETSCAPE Collaboration, D. Everett et al., “Multisystem Bayesian constraints on the transport coefficients of QCD matter”, Phys. Rev. C 103 no. 5, (2021) 054904, arXiv:2011.01430 [hep-ph].
  * (10) ALICE Collaboration, “The ALICE experiment – A journey through QCD”, arXiv:2211.04384 [nucl-ex].
  * (11) W. Broniowski, M. Chojnacki, and L. Obara, “Size fluctuations of the initial source and the event-by-event transverse momentum fluctuations in relativistic heavy-ion collisions”, Phys. Rev. C 80 (2009) 051902, arXiv:0907.3216 [nucl-th].
  * (12) P. Bozek and W. Broniowski, “Transverse-momentum fluctuations in relativistic heavy-ion collisions from event-by-event viscous hydrodynamics”, Phys. Rev. C 85 (2012) 044910, arXiv:1203.1810 [nucl-th].
  * (13) G. Giacalone, F. G. Gardim, J. Noronha-Hostler, and J.-Y. Ollitrault, “Skewness of mean transverse momentum fluctuations in heavy-ion collisions”, Phys. Rev. C 103 no. 2, (2021) 024910, arXiv:2004.09799 [nucl-th].
  * (14) B. Schenke, C. Shen, and D. Teaney, “Transverse momentum fluctuations and their correlation with elliptic flow in nuclear collision”, Phys. Rev. C 102 no. 3, (2020) 034905, arXiv:2004.00690 [nucl-th].
  * (15) F. G. Gardim, G. Giacalone, M. Luzum, and J.-Y. Ollitrault, “Effects of initial state fluctuations on the mean transverse momentum”, Nucl. Phys. A 1005 (2021) 121999, arXiv:2002.07008 [nucl-th].
  * (16) F. G. Gardim, G. Giacalone, M. Luzum, and J.-Y. Ollitrault, “Thermodynamics of hot strong-interaction matter from ultrarelativistic nuclear collisions”, Nature Phys. 16 no. 6, (2020) 615–619, arXiv:1908.09728 [nucl-th].
  * (17) B. Bally et al., “Imaging the initial condition of heavy-ion collisions and nuclear structure across the nuclide chart”, arXiv:2209.11042 [nucl-ex].
  * (18) S. Bhatta, C. Zhang, and J. Jia, “Higher-order transverse momentum fluctuations in heavy-ion collisions”, Phys. Rev. C 105 no. 2, (2022) 024904, arXiv:2112.03397 [nucl-th].
  * (19) Z. Lu, M. Zhao, J. Jia, and Y. Zhou, “Probe nuclear structure using the anisotropic flow at the Large Hadron Collider”, arXiv:2309.09663 [nucl-th].
  * (20) G. Giacalone, “Many-body correlations for nuclear physics across scales: from nuclei to quark-gluon plasmas to hadron distributions”, arXiv:2305.19843 [nucl-th].
  * (21) G. Giacalone, “Observing the deformation of nuclei with relativistic nuclear collisions”, Phys. Rev. Lett. 124 no. 20, (2020) 202301, arXiv:1910.04673 [nucl-th].
  * (22) B. Bally, M. Bender, G. Giacalone, and V. Somà, “Evidence of the triaxial structure of 129Xe at the Large Hadron Collider”, arXiv:2108.09578 [nucl-th].
  * (23) ALICE Collaboration, S. Acharya et al., “Characterizing the initial conditions of heavy-ion collisions at the LHC with mean transverse momentum and anisotropic flow correlations”, Phys. Lett. B 834 (2022) 137393, arXiv:2111.06106 [nucl-ex].
  * (24) B. Bally, G. Giacalone, and M. Bender, “The shape of gold”, Eur. Phys. J. A 59 no. 3, (2023) 58, arXiv:2301.02420 [nucl-th].
  * (25) N. Borghini, P. M. Dinh, and J.-Y. Ollitrault, “A New method for measuring azimuthal distributions in nucleus-nucleus collisions”, Phys. Rev. C 63 (2001) 054906, arXiv:nucl-th/0007063.
  * (26) J.-Y. Ollitrault, “Measures of azimuthal anisotropy in high-energy collisions”, Eur. Phys. J. A 59 no. 10, (2023) 236, arXiv:2308.11674 [nucl-ex].
  * (27) A. Bilandzic, C. H. Christensen, K. Gulbrandsen, A. Hansen, and Y. Zhou, “Generic framework for anisotropic flow analyses with multiparticle azimuthal correlations”, Phys. Rev. C89 no. 6, (2014) 064904, arXiv:1312.3572 [nucl-ex].
  * (28) Z. Moravcova, K. Gulbrandsen, and Y. Zhou, “Generic algorithm for multiparticle cumulants of azimuthal correlations in high energy nucleus collisions”, Phys. Rev. C 103 no. 2, (2021) 024913, arXiv:2005.07974 [nucl-th].
  * (29) J. Jia, “Probing triaxial deformation of atomic nuclei in high-energy heavy ion collisions”, Phys. Rev. C 105 no. 4, (2022) 044905, arXiv:2109.00604 [nucl-th].
  * (30) ALICE Collaboration, K. Aamodt et al., “Transverse momentum spectra of charged particles in proton-proton collisions at $\sqrt{s}=900$~GeV with ALICE at the LHC”, Phys. Lett. B 693 (2010) 53–68, arXiv:1007.0719 [hep-ex].
  * (31) P. J. Smith, “A recursive formulation of the old problem of obtaining moments from cumulants and vice versa”, The American Statistician 49 no. 2, (1995) 217–218. http://www.jstor.org/stable/2684642.
  * (32) X.-N. Wang and M. Gyulassy, “Energy and centrality dependence of rapidity densities at RHIC”, Phys. Rev. Lett. 86 (2001) 3496–3499, arXiv:nucl-th/0008014.
  * (33) Z.-W. Lin, C. M. Ko, B.-A. Li, B. Zhang, and S. Pal, “A Multi-phase transport model for relativistic heavy ion collisions”, Phys. Rev. C72 (2005) 064901, arXiv:nucl-th/0411110 [nucl-th].
  * (34) B. Zhang, “ZPC 1.0.1: A Parton cascade for ultrarelativistic heavy ion collisions”, Comput. Phys. Commun. 109 (1998) 193–206, arXiv:nucl-th/9709009.
  * (35) L.-W. Chen and C. M. Ko, “System size dependence of elliptic flows in relativistic heavy-ion collisions”, Phys. Lett. B634 (2006) 205–209, arXiv:nucl-th/0505044 [nucl-th].
  * (36) B.-A. Li and C. M. Ko, “Formation of superdense hadronic matter in high-energy heavy ion collisions”, Phys. Rev. C52 (1995) 2037–2063, arXiv:nucl-th/9505016 [nucl-th].
  * (37) ALICE Collaboration, S. Acharya et al., “Anisotropic flow in Xe-Xe collisions at $\mathbf{\sqrt{s_{\mathrm{NN}}}=5.44}$ TeV”, Phys. Lett. B 784 (2018) 82–95, arXiv:1805.01832 [nucl-ex].
  * (38) P. Bozek, “Transverse-momentum–flow correlations in relativistic heavy-ion collisions”, Phys. Rev. C 93 no. 4, (2016) 044908, arXiv:1601.04513 [nucl-th].
  * (39) G. Giacalone, F. G. Gardim, J. Noronha-Hostler, and J.-Y. Ollitrault, “Correlation between mean transverse momentum and anisotropic flow in heavy-ion collisions”, Phys. Rev. C 103 no. 2, (2021) 024909, arXiv:2004.01765 [nucl-th].
  * (40) X.-N. Wang and M. Gyulassy, “HIJING: A Monte Carlo model for multiple jet production in p p, p A and A A collisions”, Phys. Rev. D 44 (1991) 3501–3516.

##  Appendix A Multi-particle $p_{\mathrm{T}}$ correlations

$\displaystyle[p_{\mathrm{T}}^{(1)}]$ $\displaystyle=\frac{P_{1}}{W_{1}}$ (19)
$\displaystyle[p_{\mathrm{T}}^{(2)}]$
$\displaystyle=\frac{P_{1}^{2}-P_{2}}{W_{1}^{2}-W_{2}}$ (20)
$\displaystyle[p_{\mathrm{T}}^{(3)}]$
$\displaystyle=\frac{P_{1}^{3}-3P_{2}P_{1}+2P_{3}}{W_{1}^{3}-3W_{2}W_{1}+2W_{3}}$
(21) $\displaystyle[p_{\mathrm{T}}^{(4)}]$
$\displaystyle=\frac{P_{1}^{4}-6P_{2}P_{1}^{2}+3P_{2}^{2}+8P_{3}P_{1}-6P_{4}}{W_{1}^{4}-6W_{2}W_{1}^{2}+3W_{2}^{2}+8W_{3}W_{1}-6W_{4}}$
(22) $\displaystyle[p_{\mathrm{T}}^{(5)}]$
$\displaystyle=\left[P_{1}^{5}-10P_{2}P_{1}^{3}+15P_{2}^{2}P_{1}\right.$ (23)
$\displaystyle\left.+20P_{3}P_{1}^{2}-20P_{3}P_{2}-30P_{4}P_{1}+24P_{5}\right]$
$\displaystyle/\left[W_{1}^{5}-10W_{2}W_{1}^{3}+15W_{2}^{2}W_{1}+20W_{3}W_{1}^{2}\right.$
$\displaystyle\left.-20W_{3}W_{2}-30W_{4}W_{1}+24W_{5}\right]$
$\displaystyle[p_{\mathrm{T}}^{(6)}]$
$\displaystyle=\left[P_{1}^{6}-15P_{2}P_{1}^{4}+45P_{1}^{2}P_{2}^{2}\right.$
(24) $\displaystyle-15P_{2}^{3}+40P_{3}P_{1}^{3}-120P_{3}P_{2}P_{1}$
$\displaystyle+40P_{3}^{2}-90P_{4}P_{1}^{2}+90P_{4}P_{2}$
$\displaystyle\left.+144P_{5}P_{1}-120P_{6}\right]/\left[W_{1}^{6}-\right.$
$\displaystyle 15W_{2}W_{1}^{4}+45W_{1}^{2}W_{2}^{2}-15W_{2}^{3}$
$\displaystyle+40W_{3}W_{1}^{3}-120W_{3}W_{2}W_{1}$
$\displaystyle+40W_{3}^{2}-90W_{4}W_{1}^{2}+90W_{4}W_{2}$
$\displaystyle\left.+144W_{5}W_{1}-120W_{6}\right]$
$\displaystyle[p_{\mathrm{T}}^{(7)}]$
$\displaystyle=\left[P_{1}^{7}-21P_{2}P_{1}^{5}+105P_{1}^{3}P_{2}^{2}\right.$
(25) $\displaystyle-105P_{2}^{3}P_{1}+70P_{3}P_{1}^{4}-420P_{3}P_{2}P_{1}^{2}$
$\displaystyle+210P_{3}P_{2}^{2}+280P_{3}^{2}P_{1}-210P_{4}P_{1}^{3}$
$\displaystyle+630P_{4}P_{2}P_{1}-420P_{4}P_{3}+504P_{5}P_{1}^{2}$
$\displaystyle\left.-504P_{5}P_{2}-840P_{6}P_{1}+720P_{7}\right]/\left[W_{1}^{7}\right.$
$\displaystyle-21W_{2}W_{1}^{5}+105W_{1}^{3}W_{2}^{2}$
$\displaystyle-105W_{2}^{3}W_{1}+70W_{3}W_{1}^{4}-420W_{3}W_{2}W_{1}^{2}$
$\displaystyle+210W_{3}W_{2}^{2}+280W_{3}^{2}W_{1}-210W_{4}W_{1}^{3}$
$\displaystyle+630W_{4}W_{2}W_{1}-420W_{4}W_{3}+504W_{5}W_{1}^{2}$
$\displaystyle\left.-504W_{5}W_{2}-840W_{6}W_{1}+720W_{7}\right]$
$\displaystyle[p_{\mathrm{T}}^{(8)}]$
$\displaystyle=\left[P_{1}^{8}-28P_{2}P_{1}^{6}+210P_{2}^{2}P_{1}^{4}\right.$
(26) $\displaystyle-420P_{2}^{3}P_{1}^{2}+105P_{2}^{4}+112P_{3}P_{1}^{5}$
$\displaystyle-1120P_{3}P_{2}P_{1}^{3}+1680P_{3}P_{2}^{2}P_{1}+1120P_{3}^{2}P_{1}^{2}$
$\displaystyle-1120P_{3}^{2}P_{2}-420P_{4}P_{1}^{4}+2520P_{4}P_{2}P_{1}^{2}$
$\displaystyle-1260P_{4}P_{2}^{2}-3360P_{4}P_{3}P_{1}+1260P_{4}^{2}$
$\displaystyle+1344P_{5}P_{1}^{3}-4032P_{5}P_{2}P_{1}+2688P_{5}P_{3}$
$\displaystyle-3360P_{6}P_{1}^{2}+3360P_{6}P_{2}+5760P_{7}P_{1}$
$\displaystyle\left.-5040P_{8}\right]/\left[W_{1}^{8}-28W_{2}W_{1}^{6}+210W_{2}^{2}W_{1}^{4}\right.$
$\displaystyle-420W_{2}^{3}W_{1}^{2}+105W_{2}^{4}+112W_{3}W_{1}^{5}$
$\displaystyle-1120W_{3}W_{2}W_{1}^{3}+1680W_{3}W_{2}^{2}W_{1}+1120W_{3}^{2}W_{1}^{2}$
$\displaystyle-1120W_{3}^{2}W_{2}-420W_{4}W_{1}^{4}+2520W_{4}W_{2}W_{1}^{2}$
$\displaystyle-1260W_{4}W_{2}^{2}-3360W_{4}W_{3}W_{1}+1260W_{4}^{2}$
$\displaystyle+1344W_{5}W_{1}^{3}-4032W_{5}W_{2}W_{1}+2688W_{5}W_{3}$
$\displaystyle-3360W_{6}W_{1}^{2}+3360W_{6}W_{2}+5760W_{7}W_{1}$
$\displaystyle\left.-5040W_{8}\right]$

##  Appendix B Multi-particle $p_{\mathrm{T}}$ cumulants

$\displaystyle\kappa_{1}$ $\displaystyle=\langle[p_{\mathrm{T}}]\rangle$ (27)
$\displaystyle\kappa_{2}$
$\displaystyle=\langle[p_{\mathrm{T}}^{(2)}]\rangle-\langle[p_{\mathrm{T}}]\rangle^{2}$
(28) $\displaystyle\kappa_{3}$
$\displaystyle=\langle[p_{\mathrm{T}}^{(3)}]\rangle-3\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle+2\langle[p_{\mathrm{T}}]\rangle^{3}$
(29) $\displaystyle\kappa_{4}$
$\displaystyle=\langle[p_{\mathrm{T}}^{(4)}]\rangle-4\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}]\rangle-3\langle[p_{\mathrm{T}}^{(2)}]\rangle^{2}$
(30)
$\displaystyle+12\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{2}-6\langle[p_{\mathrm{T}}]\rangle^{4}$
$\displaystyle\kappa_{5}$
$\displaystyle=\langle[p_{\mathrm{T}}^{(5)}]\rangle-5\langle[p_{\mathrm{T}}^{(4)}]\rangle\langle[p_{\mathrm{T}}]\rangle-10\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle$
(31)
$\displaystyle+30\langle[p_{\mathrm{T}}^{(2)}]\rangle^{2}\langle[p_{\mathrm{T}}]\rangle+20\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{2}$
$\displaystyle-60\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{3}+24\langle[p_{\mathrm{T}}]\rangle^{5}$
$\displaystyle\kappa_{6}$
$\displaystyle=\langle[p_{\mathrm{T}}^{(6)}]\rangle-6\langle[p_{\mathrm{T}}^{(5)}]\rangle\langle[p_{\mathrm{T}}]\rangle-15\langle[p_{\mathrm{T}}^{(4)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle$
(32)
$\displaystyle-10\langle[p_{\mathrm{T}}^{(3)}]\rangle^{2}+30\langle[p_{\mathrm{T}}^{(2)}]\rangle^{3}+30\langle[p_{\mathrm{T}}^{(4)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{2}$
$\displaystyle+120\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle-270\langle[p_{\mathrm{T}}^{(2)}]\rangle^{2}\langle[p_{\mathrm{T}}]\rangle^{2}$
$\displaystyle-120\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{3}+360\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{4}$
$\displaystyle-120\langle[p_{\mathrm{T}}]\rangle^{6}$
$\displaystyle\kappa_{7}$
$\displaystyle=\langle[p_{\mathrm{T}}^{(7)}]\rangle-7\langle[p_{\mathrm{T}}^{(6)}]\rangle\langle[p_{\mathrm{T}}]\rangle-21\langle[p_{\mathrm{T}}^{(5)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle$
(33)
$\displaystyle+42\langle[p_{\mathrm{T}}^{(5)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{2}-35\langle[p_{\mathrm{T}}^{(4)}]\rangle\langle[p_{\mathrm{T}}^{(3)}]\rangle$
$\displaystyle+210\langle[p_{\mathrm{T}}^{(4)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle-210\langle[p_{\mathrm{T}}^{(4)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{3}$
$\displaystyle+140\langle[p_{\mathrm{T}}^{(3)}]\rangle^{2}\langle[p_{\mathrm{T}}]\rangle+210\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle^{2}$
$\displaystyle-1260\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{2}+840\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{4}$
$\displaystyle-630\langle[p_{\mathrm{T}}^{(2)}]\rangle^{3}\langle[p_{\mathrm{T}}]\rangle+2520\langle[p_{\mathrm{T}}^{(2)}]\rangle^{2}\langle[p_{\mathrm{T}}]\rangle^{3}$
$\displaystyle-2520\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{5}+720\langle[p_{\mathrm{T}}]\rangle^{7}$
$\displaystyle\kappa_{8}$
$\displaystyle=\langle[p_{\mathrm{T}}^{(8)}]\rangle-8\langle[p_{\mathrm{T}}^{(7)}]\rangle\langle[p_{\mathrm{T}}]\rangle$
(34)
$\displaystyle-28\langle[p_{\mathrm{T}}^{(6)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle+56\langle[p_{\mathrm{T}}^{(6)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{2}$
$\displaystyle-56\langle[p_{\mathrm{T}}^{(5)}]\rangle\langle[p_{\mathrm{T}}^{(3)}]\rangle+336\langle[p_{\mathrm{T}}^{(5)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle$
$\displaystyle-336\langle[p_{\mathrm{T}}^{(5)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{3}-35\langle[p_{\mathrm{T}}^{(4)}]\rangle^{2}$
$\displaystyle+560\langle[p_{\mathrm{T}}^{(4)}]\rangle\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}]\rangle+420\langle[p_{\mathrm{T}}^{(2)}]\rangle^{2}\langle[p_{\mathrm{T}}^{(4)}]\rangle$
$\displaystyle-2520\langle[p_{\mathrm{T}}^{(4)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{2}+1680\langle[p_{\mathrm{T}}^{(4)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{4}$
$\displaystyle+560\langle[p_{\mathrm{T}}^{(3)}]\rangle^{2}\langle[p_{\mathrm{T}}^{(2)}]\rangle-1680\langle[p_{\mathrm{T}}^{(3)}]\rangle^{2}\langle[p_{\mathrm{T}}]\rangle^{2}$
$\displaystyle-5040\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle^{2}\langle[p_{\mathrm{T}}]\rangle+13440\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{3}$
$\displaystyle-6720\langle[p_{\mathrm{T}}^{(3)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{5}-630\langle[p_{\mathrm{T}}^{(2)}]\rangle^{4}+10080\langle[p_{\mathrm{T}}^{(2)}]\rangle^{3}\langle[p_{\mathrm{T}}]\rangle^{2}$
$\displaystyle-25200\langle[p_{\mathrm{T}}^{(2)}]\rangle^{2}\langle[p_{\mathrm{T}}]\rangle^{4}+20160\langle[p_{\mathrm{T}}^{(2)}]\rangle\langle[p_{\mathrm{T}}]\rangle^{6}-5040\langle[p_{\mathrm{T}}]\rangle^{8}$