# Brain Signals to Rescue Aphasia, Apraxia and Dysarthria Speech Recognition

Gautam Krishna1, Mason Carnahan1, Shilpa Shamapant2, Yashitha Surendranath1,∗,
Saumya Jain1,∗, Arundhati Ghosh1,∗,
Co Tran1, Jose del R Millan3 and Ahmed H Tewfik1 1Author is with Brain Machine
Interface Lab, The University of Texas at Austin, Austin, Texas, USA 2Author
is with Austin Speech Labs, Austin, Texas, USA 3Author is with Clinical
Neuroprosthetics and Brain Interaction Lab, The University of Texas at Austin,
Austin, Texas, USA ∗Equal author contribution

###### Abstract

In this paper, we propose a deep learning-based algorithm to improve the
performance of automatic speech recognition (ASR) systems for aphasia,
apraxia, and dysarthria speech by utilizing electroencephalography (EEG)
features recorded synchronously with aphasia, apraxia, and dysarthria speech.
We demonstrate a significant decoding performance improvement by more than 50%
during test time for isolated speech recognition task and we also provide
preliminary results indicating performance improvement for the more
challenging continuous speech recognition task by utilizing EEG features. The
results presented in this paper show the first step towards demonstrating the
possibility of utilizing non-invasive neural signals to design a real-time
robust speech prosthetic for stroke survivors recovering from aphasia,
apraxia, and dysarthria. Our aphasia, apraxia, and dysarthria speech-EEG data
set will be released to the public to help further advance this interesting
and crucial research.

##  I INTRODUCTION

Automatic speech recognition (ASR) system converts speech to text and it forms
the back-end in many state-of-the-art virtual voice assistants like Apple’s
Siri, Amazon’s Alexa, Samsung’s Bixby, etc. These voice assistants are trained
to recognize the uniform speech of users with no speech disorders. The
performance of ASR systems degrades in presence of incomplete, distorted, or
broken speech. This limits technology accessibility to users with speech
disorders. The three most common speech, language disorders are aphasia,
apraxia, and dysarthria. Aphasia is a disturbance of the comprehension and
formulation of language caused by dysfunction in specific brain regions. The
major causes are a stroke or head trauma[1, 2]. Apraxia is a speech disorder
caused due to the impairment of motor planning of speech [3]. Dysarthria is
also a speech disorder caused due to neurological damage to the motor
component of the motor–speech system and it is closely related to Apraxia [4].
People recovering from these speech disorders produce distorted and incomplete
speech. The work described by authors in [5, 6] demonstrate that
electrophysiological monitoring of neural signals like electroencephalography
(EEG) and electrocorticography (ECoG) carry important information about speech
articulation and speech perception. They demonstrated the results using neural
signals recorded from subjects with no speech disorders. In [7] authors
demonstrated that EEG features can be used to enhance the performance of
isolated speech recognition systems trained to decode speech of users with no
speech disorders. In their work, they demonstrated results on an English
vocabulary consisting of four words and five vowels. EEG is a non-invasive way
of measuring the electrical activity of the human brain. The EEG sensors are
placed on the scalp of the subject to obtain EEG recordings. The EEG signals
offer a very high temporal resolution. The non-invasive nature of EEG signals
makes it safe and easy to deploy even-though EEG signals offer poor spatial
resolution and signal-to-noise ratio compared to invasive ECoG neural activity
recording techniques. The high temporal resolution property of EEG signals
also allows capturing the human speech-related neural activities as normal
human speech occurs at a high rate of 150 words per minute. In [8] authors
explored speech recognition using aphasia speech and reported a very high word
error rate (WER) during test time. For a reduced vocabulary, they reported a
WER as high as 97.5 %. In [9] authors demonstrated aphasia speech recognition
by training their acoustic models on a large scale aphasia speech data-set
named AphasiaBank but they reported a high phoneme error rate (PER) in the
range of 75% to 89% for severe cases of aphasia. A high PER indicates an even
higher WER. In a very recent work described in [10] authors explored the
possibility of using ASR systems as a feedback tool while providing speech
therapy to aphasia patients. Their results demonstrated an increase in the
effectiveness of the speech therapy when coupled with ASR technology.
References [11, 12, 13] investigated speech recognition for apraxia and
dysarthria speech and reported low accuracy on a word-level vocabulary. In
[14] authors carried out an EEG study to analyze the EEG delta wavebands to
understand the brain damage on patients recovering from aphasia. In related
studies described in references [15, 16] authors investigated EEG activity in
the left-hemisphere of the brain of subjects recovering from aphasia and an
EEG sleep study to understand the brain activity of the aphasia patients.
These studies demonstrated that EEG signals carried useful information about
brain function recovery in aphasia patients. In this paper, we propose an
algorithm to train a deep learning-based speech recognizer using acoustic
features along with acoustic representations derived from EEG features to
significantly improve the test time decoding performance of aphasia + apraxia
+ dysarthria isolated speech recognizer. We were able to achieve a performance
improvement of more than 50% during test time for the task of isolated speech
recognition and a slight improvement in performance for the more challenging
task of continuous speech recognition using our proposed algorithm. The
results presented in this paper demonstrate how non-invasive neural signals
can be utilized to improve the performance of speech recognizers used to
decode aphasia, apraxia, and dysarthria speech. Designing a speech recognizer
that can decode aphasia, apraxia, and dysarthria speech with high accuracy has
the potentiality to lead to a design of a speech prosthetic and a better
speech therapy tool for stroke survivors.

Our main contributions and major highlights of our proposed algorithm are
listed below:

  * •

We developed a deep learning-based algorithm to improve the performance of
speech recognition for aphasia, apraxia, and dysarthria speech by utilizing
EEG features.

  * •

We collected large-scale aphasia, apraxia and dysarthria Speech-EEG data set
that will be released to the public to help further advance this research.

  * •

Our proposed algorithm can be used with any type of speech recognition model.

  * •

Our proposed algorithm can be easily extended to other tasks like speaker
identification, speech detection, etc for aphasia, apraxia, and dysarthria
speech by replacing the ASR model in our algorithm with the new objective
model.

##  II Proposed Deep learning algorithm to improve Speech Recognition

Figure 1 describes the architecture of our proposed deep learning training
strategy to improve the ASR performance of aphasia, apraxia, and dysarthria
speech by utilizing EEG features. As seen from the figure, we make use of an
EEG to acoustic feature mapping, regression model to generate additional
features that are provided to the ASR model to improve its training. We first
train the regression model described on the right-hand side of the figure to
predict acoustic features or Mel frequency cepstral coefficients (MFCC) [17]
of dimension 13 from EEG features. The regression model consists of a single
layer of gated recurrent unit (GRU) [18] with 128 hidden units connected to a
time distributed dense layer consisting of 13 hidden units with a linear
activation function. The regression model was trained for 70 epochs with mean
square error (MSE) as the loss function and with adam [19] as the optimizer.
The batch size was set to 100. The GRU layer in the regression model learns
the acoustic representation present in the input EEG features. We then
concatenate these acoustic representations or outputs of the GRU layer of the
regression model with the input acoustic or MFCC features of dimension 13
which are then used to train the ASR model to produce the text output during
training time. The ASR model is trained after completing the training of the
regression model. During test time, the EEG features from the test set are
provided as input to the trained regression model, and the output of the GRU
layer of the regression model is concatenated with the simultaneously recorded
MFCC features from the test set to produce text output from the trained ASR
model. The output of the GRU layer of the regression model or the acoustic
representations present in EEG features is of dimension 128. The choice of the
ASR model architecture depends on the task. We investigated both the tasks of
isolated and continuous speech recognition in this paper. Isolated speech
recognition refers to a sentence or sequence classification task, where the
model decodes closed vocabulary and directly learns the input feature to
sentence mapping. Here the model predicts the complete sentence or label token
as output per decoding stage. On the other hand, continuous speech recognition
refers to the task where the model is predicting the text by predicting the
character or word or phoneme at every time-step and these models are capable
of performing open vocabulary decoding. Continuous speech recognition is a
more challenging task than isolated speech recognition due to the increase in
the number of model parameters and learning alignments.

Figure 1: Proposed Training Algorithm Figure 2: Isolated Speech Recognition
Model

Next, we briefly describe the architecture of the isolated and continuous
speech recognition models used in this work. Our isolated speech recognition
model consists of a single layer of GRU with 512 hidden units connected to a
dropout regularization [20] with a drop-out rate of 0.2. The drop-out
regularization is followed by a dense layer consisting of 57 hidden units and
a linear activation function. The dense layer contained 57 hidden units since
our vocabulary contained 57 unique English sentences. The last time-step
output of the GRU layer is passed to dropout regularization and dense layer.
Finally, the dense layer output or logits are passed through a softmax
activation function to obtain the label prediction probabilities. Each label
token corresponds to a complete English sentence text. The labels were one-hot
vector encoded and the model was trained for 10 epochs with batch size set to
50. We used early stopping during training to prevent over-fitting. We used
categorical cross-entropy as the loss function and adam was used as the
optimizer. The model architecture is described in Figure 2. Our continuous
speech recognition model consists of a GRU layer with 512 hidden units acting
as an encoder and the decoder consists of a combination of a dense layer with
linear activation function and softmax activation function. The output of the
encoder is passed to the decoder at every time-step. The model was trained for
100 epochs with batch size set to 50 to optimize connectionist temporal
classification (CTC) loss function [21, 22]. We used adam as the optimizer.
For this work, a character-based CTC model was used. The model was predicting
a character at every time-step. We used an external 4-gram language model
along with a CTC beam search decoder during inference time [23].

Figure 3: Continuous Speech Recognition Model

##  III Design of experiments for building the data base

Nine subjects with a known diagnosis of aphasia or apraxia or dysarthria or a
combination of any of these disorders volunteered to take part in our data
collection experiments. All experimental procedures were approved by the
Institutional Review Board at the University of Texas at Austin. The
demographic information of the subjects is shown below in Table<EMAIL_ADDRESS>Each
subject was asked to perform two different tasks while they were receiving
speech therapy at Austin Speech Labs. The first task involved subjects reading
out loud English sentences shown to them on a computer screen and their
simultaneous EEG, electromyography (EMG), and speech signals were recorded.
The second task involved subjects listening to the recorded audio of English
sentences and they were then asked to speak out loud what they listened to and
their simultaneous EEG, EMG, and speech signals were recorded. We collected a
total of 8854 data samples from the 9 subjects for both the tasks combined.
The vocabulary consisted of 57 unique daily used common English sentences. We
used brain products wet EEG sensors for this data collection. We used 29 EEG
sensors in total. The placement of 29 sensors was based on the standard 10-20
EEG sensor placement guidelines. Figure 4 shows a subject wearing our EEG cap
during the experiment. We used the brain product’s Actchamp amplifier as the
EEG amplifier. We further used two EMG sensors to keep track of EMG artifacts
generated during articulation. The EMG sensor placement location is shown in
Figure 5. The speech signals were recorded using a mono-channel microphone. We
used 70% of the data as the training set, 10% as the validation set, and the
remaining 20% as the test set. The data set split was done randomly using the
scikit-learn train-test split python function. There was no overlap between
training, validation, and test set data points.

ID |  | Male/
---
Female
| Aphasia
---
Quotient
| Aphasia
---
Type
Severity | Speech Disorders
1 | M | 48.4 | global | severe | aphasia, apraxia
2 | F | 53.2 | global | severe | aphasia
3 | M | 74.2 | Broca’s | moderate | aphasia
4 | M | 87.6 | Anomia | mild | aphasia
5 | M | 86.2 | Broca’s | moderate | aphasia, apraxia
6 | M | 94.8 | Anomic | mild | aphasia, dysarthria
7 | M | 48 | global | severe | aphasia
8 | M | 87.4 | Broca’s | moderate | aphasia, apraxia
9 | M | 72.4 | mixed | mild | aphasia
TABLE I: Data-set Demographics Figure 4: EEG sensor placement Figure 5: EMG
sensor placement

##  IV EEG and Speech feature extraction details

The recorded EEG signals were sampled at a sampling frequency of 1000Hz and a
fourth-order IIR bandpass filter with cut-off frequencies 0.1Hz and 70Hz was
applied. A notch filter with a cut off frequency of 60 Hz was used to remove
the power line noise. We then used the linear regression technique to remove
EMG artifacts from EEG signals.
${Corrected}_{EEG}$ = ${Recorded}_{EEG}-\alpha\ast{Recorded}_{EMG}$, where
$\alpha$ is the regression coefficient computed by Ordinary Least Squares
method. We then extracted five features per EEG channel. The five features
extracted were root mean square, zero-crossing rate, moving window average,
kurtosis, and power spectral entropy [7, 24]. This EEG feature set was first
introduced by authors in [7] where they demonstrated that these features carry
neural information about speech perception and production. The EEG features
were extracted at a sampling frequency of 100 Hz per channel. The speech
signal was recorded at a sampling frequency of 16KHz. We extracted Mel
frequency cepstral coefficients (MFCC) [17] of dimension 13 as features for
speech signal. The MFCC features were also extracted at the same sampling
frequency 100Hz as that of EEG feature extraction.

Figure 6: Explained Variance vs No of Components

##  V EEG Feature Dimension Reduction Algorithm Details

Since the local structure of our EEG feature space was not linear, we used
non-linear dimension reduction technique to perform dimension reduction on EEG
features. We plotted cumulative explained variance vs the number of components
as shown in Figure 6 to identify the optimal EEG feature space dimension. We
used kernel principal component analysis (KPCA) [25] with a polynomial kernel
of degree 3 to reduce our EEG feature space of dimension 145 (five features
per each of the 29 channels) to a final dimension of 10. Before applying KPCA,
the EEG features were normalized by removing the mean and scaling to unit
variance.

##  VI Results and Discussion

We used percentage accuracy, F1-score, precision, and recall [26] as
performance metrics to evaluate the performance of the isolated speech
recognition model. The higher the accuracy, F1-score, precision, and recall
values the better the performance of the model. For computing F1-score,
precision and recall we added a small value e-07 called epsilon to the
denominator of F1-score, precision and recall expressions to prevent a divide
by zero error. We used word error rate (WER) as the performance metric to
evaluate the performance of the continuous speech recognition model. The lower
the WER value, the better the speech recognition system performance. For
obtaining baseline results, the speech recognition models were trained and
tested using only acoustic or MFCC features. Table ii@ shows the test time
results obtained for isolated speech recognition task for various EEG
frequency bands. We compared results obtained using low-frequency EEG signals
( 0.1 Hz to 15 Hz), high-frequency EEG signals ( 15 Hz to 70 Hz), and all
frequency EEG signals ( 0.1 Hz to 70 Hz). The results shown in Table ii@
demonstrate that choice of EEG frequency range had less effect on decoding
performance for the isolated speech recognition task. The work carried out by
authors in [6] demonstrated that both high and low-frequency neural signals
carry important information about speech production. Table iii@ shows test
times results for isolated speech recognition task with and without EMG
artifact removal and obtained results demonstrate that even though removing
EMG artifacts improved the test-time performance of the model, the improvement
was not that significant. Table iv@ shows the test time results for isolated
speech recognition task with and without EEG dimension reduction. The results
demonstrate that EEG dimension reduction using KPCA had resulted in
significant performance improvement of the model during test time. Table v@
shows the test time results for isolated speech recognition task when we used
only temporal lobe EEG sensor features, frontal lobe EEG sensor features, and
concatenation of temporal and frontal lobe EEG sensor features. The temporal
and frontal lobe contains brain regions responsible for speech perception and
production [27, 28]. EEG features from frontal lobe sensors Fp1, Fz, F3, F7,
FT9, FC5, FT10 , FC6 , FC2 , F4 , F8 and Fp2 were extracted and then reduced
to a dimension of 10 using KPCA. Similarly, EEG features were extracted from
temporal lobe sensors T7, TP9, TP10, and T8 and then reduced to a dimension of
10 using KPCA. The results shown in Table v@ demonstrate that it is possible
to achieve comparable decoding performance for isolated speech recognition
task using EEG sensors from just frontal and temporal lobe regions instead of
using all the EEG sensors. Table vi@ shows the test time results for isolated
speech recognition task when we used only the first half-length of the input
EEG and MFCC features instead of the complete length of EEG or MFCC features
for decoding text. The motivation here was to see if the model can decode text
if we provide incomplete input as most of the aphasia or apraxia speech
involves a lot of pauses in between. As seen from the Table vi@ results we
observed that when half the length of the input signal is used, the baseline
results improved significantly but adding acoustic representation in EEG to
MFCC features still outperformed the baseline for all the test-time
performance metrics. We believe the baseline results improved when shorter
sequences were used as input signal due to the fact that GRU can process
shorter sequences more efficiently than longer input sequences [18, 29]. The
overall results from Tables ii@,iii@,v@ and vi@ show that adding acoustic
representation in EEG features with MFCC features significantly outperform the
baseline for all the test-time performance metrics for the task of isolated
speech recognition using aphasia, apraxia, and dysarthria speech. Figure 7
shows the training and validation loss convergence for the regression model
and Figure 8 shows the training and validation accuracy of the isolated speech
recognition model. The training, validation loss values were comparable as
well as the training and validation accuracy values, indicating the models
didn’t overfit. Figure 9 shows the confusion matrix obtained during test time
for the isolated speech recognition task when tested using MFCC+ High-
frequency EEG representation. Each token in the confusion matrix represents a
complete English sentence from the test set. Table vii@ shows test time
results for isolated speech recognition task when acoustic features were
concatenated with acoustic representation in EMG features of dimension 10
compared to acoustic representations from EEG features of dimension 10. We
extracted the same set of 5 features that we extracted for EEG for each EMG
channel. The results show that the acoustic representations present in EMG is
not significant compared to acoustic representation features present in EEG
signals for boosting the performance of the speech recognizer. Table viii@
shows the test time average WER obtained for the continuous speech recognition
task. The obtained results demonstrate that adding acoustic representation in
High-frequency EEG features to MFCC outperformed the baseline for a test set
vocabulary consisting of 1771 English sentences. We obtained a p value [30] of
0.0000213, demonstrating high statistical significance for our result. We
further computed the test time WER’s with 95 % confidence level value and
observed that for the baseline, the WER range was between 48.58% and 51.1%
where as using our proposed method, the WER range for the same confidence
level value was between 44.25% and 47.13%. Therefore a thorough statistical
analysis of our test time continuous speech recognition results demonstrate
that our proposed method outperformed the baseline result.

| Performance
---
Metric(%)
MFCC |  | MFCC
---
+
Acoustic Representation
in Low Freq EEG dim 10
| MFCC
---
+
Acoustic Representation
in High Freq EEG dim 10
| MFCC
---
+
Acoustic Representation
in All Freq EEG dim 10
Accuracy | 28.40 | 78.93 | 81.02 | 80.74
F1 - score | 34.93 | 81.65 | 82.86 | 83.25
precision | 75.88 | 86.54 | 87.23 | 88.27
recall | 23.06 | 77.40 | 79.02 | 78.91
TABLE II: Test time results for isolated speech recognition for various EEG frequency bands | Performance
---
Metric(%)
MFCC |  | MFCC
---
+
Acoustic Representation
in All Freq EEG dim 10
EMG artifacts
Removed
| MFCC
---
+
Acoustic Representation
in All Freq EEG dim 10
No EMG artifacts
Removed
Accuracy | 28.40 | 80.74 | 80.46
F1 - score | 34.93 | 83.25 | 82.24
precision | 75.88 | 88.27 | 86.00
recall | 23.06 | 78.91 | 78.91
TABLE III: Isolated Speech recognition Test Time results - effect of EMG artifact removal Figure 7: Training and validation loss convergence of the regression model Figure 8: Training and validation accuracy of the isolated speech recognition model when trained using MFCC+ acoustic representation in High frequency EEG Figure 9: Confusion matrix of the isolated speech recognition model during test time when tested using MFCC+ acoustic representation in High frequency EEG | Performance
---
Metric(%)
MFCC |  | MFCC
---
+
Acoustic Representation
in High Freq EEG dim 10
| MFCC
---
+
Acoustic Representation
in High Freq EEG dim 145
No KPCA
Accuracy | 28.40 | 81.02 | 28.79
F1 - score | 34.93 | 82.86 | 34.56
precision | 75.88 | 87.23 | 73.66
recall | 23.06 | 79.02 | 22.98
TABLE IV: Isolated Speech recognition test time results-effect of KPCA dimension reduction | Performance
---
Metric(%)
MFCC |  | MFCC
---
+
Acoustic Representation
in High Freq EEG dim 10
Temporal Lobe
| MFCC
---
+
Acoustic Representation
in High Freq EEG dim 10
Frontal Lobe
| MFCC
---
+
Acoustic Representation
in High Freq EEG dim 20
Temporal Lobe and Frontal
Lobe
Accuracy | 28.40 | 81.25 | 81.14 | 80.68
F1 - score | 34.93 | 82.97 | 83.78 | 83.77
precision | 75.88 | 86.91 | 89.20 | 90.68
recall | 23.06 | 79.53 | 79.13 | 78.08
TABLE V: Isolated Speech recognition test time results-effect of EEG sensor reduction | Performance
---
Metric(%)
MFCC |  | MFCC
---
+
Acoustic Representation
in High Freq EEG dim 10
Accuracy | 78.09 | 79.84
F1 - score | 80.69 | 82.03
precision | 86.02 | 87.29
recall | 76.12 | 77.46
TABLE VI: Isolated Speech recognition test time results when first half length of the complete speech and EEG signals are used as input | Performance
---
Metric(%)
MFCC |  | MFCC
---
+
Acoustic Representation
in High Freq EMG dim 10
| MFCC
---
+
Acoustic Representation
in High Freq EEG dim 10
Accuracy | 28.40 | 30.15 | 81.02
F1-score | 34.93 | 37.02 | 82.86
precision | 75.88 | 80.71 | 87.23
recall | 23.06 | 24.42 | 79.02
TABLE VII: Isolated Speech recognition test time results-effect of EMG vs EEG sensors Number of sentences |  | MFCC
---
(WER %)
| MFCC
---
+
Acoustic Representation
in High Freq EEG dim 10
(WER %)
1771 | 49.84 | 45.69
TABLE VIII: Continuous Speech Recognition Test time results

##  VII Conclusion, Limitation and Future work

In this paper, we proposed a deep learning based algorithm to improve the
performance of isolated and continuous speech recognition systems for aphasia,
apraxia, and dysarthria speech by utilizing non-invasive neural EEG signals
recorded synchronously with the speech. Our proposed algorithm outperformed
the baseline results for the task of isolated speech recognition during test
time by more than 50% and at the same time outperforming the baseline results
for the more challenging task of continuous speech recognition by a small
margin. To the best of our knowledge, this is the first work that demonstrates
how to utilize non-invasive neural signals to improve the decoding performance
of speech recognition systems for aphasia, apraxia, and dysarthria speech. One
major limitation of the proposed algorithm is the latency that might be
observed when this system is deployed in real-time as the first step is to
obtain the acoustic representations in EEG using the trained regression model
before it is concatenated with the acoustic features to decode text. All the
results presented in this paper are based on the offline analysis. The latency
will be a function of the input sequence length, model size, and computational
resources (GPU memory and RAM). Our future work will focus on validating these
results on larger data set as we make progress in our data collection efforts.
Future work will also focus on performing more experiments for the task of
continuous speech recognition and developing tools to improve the performance
of our proposed algorithm. Our aphasia, apraxia, and dysarthria speech-EEG
data set will be released to the public to help further advance this
interesting and crucial research.

## References

  * [1] A. R. Damasio, “Aphasia,” _New England Journal of Medicine_ , vol. 326, no. 8, pp. 531–539, 1992.
  * [2] D. F. Benson and A. Ardila, _Aphasia: A clinical perspective_. Oxford University Press on Demand, 1996.
  * [3] R. D. Kent and J. C. Rosenbek, “Acoustic patterns of apraxia of speech,” _Journal of Speech, Language, and Hearing Research_ , vol. 26, no. 2, pp. 231–249, 1983.
  * [4] F. L. Darley, A. E. Aronson, and J. R. Brown, “Differential diagnostic patterns of dysarthria,” _Journal of speech and hearing research_ , vol. 12, no. 2, pp. 246–269, 1969.
  * [5] G. Krishna, C. Tran, Y. Han, M. Carnahan, and A. Tewfik, “Speech synthesis using eeg,” in _Acoustics, Speech and Signal Processing (ICASSP), 2020 IEEE International Conference on_. IEEE, 2020.
  * [6] G. K. Anumanchipalli, J. Chartier, and E. F. Chang, “Speech synthesis from neural decoding of spoken sentences,” _Nature_ , vol. 568, no. 7753, p. 493, 2019.
  * [7] G. Krishna, C. Tran, J. Yu, and A. Tewfik, “Speech recognition with no speech or with noisy speech,” in _Acoustics, Speech and Signal Processing (ICASSP), 2019 IEEE International Conference on_. IEEE, 2019.
  * [8] K. C. Fraser, F. Rudzicz, N. Graham, and E. Rochon, “Automatic speech recognition in the diagnosis of primary progressive aphasia,” in _Proceedings of the fourth workshop on speech and language processing for assistive technologies_ , 2013, pp. 47–54.
  * [9] D. Le and E. M. Provost, “Improving automatic recognition of aphasic speech with aphasiabank.” in _Interspeech_ , 2016, pp. 2681–2685.
  * [10] K. J. Ballard, N. M. Etter, S. Shen, P. Monroe, and C. Tien Tan, “Feasibility of automatic speech recognition for providing feedback during tablet-based treatment for apraxia of speech plus aphasia,” _American journal of speech-language pathology_ , vol. 28, no. 2S, pp. 818–834, 2019.
  * [11] A. Jacks, K. L. Haley, G. Bishop, and T. G. Harmon, “Automated speech recognition in adult stroke survivors: Comparing human and computer transcriptions,” _Folia Phoniatrica et Logopaedica_ , vol. 71, no. 5-6, pp. 286–296, 2019.
  * [12] P. Green, J. Carmichael, A. Hatzis, P. Enderby, M. Hawley, and M. Parker, “Automatic speech recognition with sparse training data for dysarthric speakers,” in _Eighth European conference on speech communication and technology_ , 2003.
  * [13] L. Ferrier, H. Shane, H. Ballard, T. Carpenter, and A. Benoit, “Dysarthric speakers’ intelligibility and speech characteristics in relation to computer speech recognition,” _Augmentative and Alternative Communication_ , vol. 11, no. 3, pp. 165–175, 1995.
  * [14] C. Spironelli and A. Angrilli, “Eeg delta band as a marker of brain damage in aphasic patients after recovery of language,” _Neuropsychologia_ , vol. 47, no. 4, pp. 988–994, 2009.
  * [15] S. Hensel, B. Rockstroh, P. Berg, T. Elbert, and P. W. Schönle, “Left-hemispheric abnormal eeg activity in relation to impairment and recovery in aphasic patients,” _Psychophysiology_ , vol. 41, no. 3, pp. 394–400, 2004.
  * [16] S. Sarasso, S. Määttä, F. Ferrarelli, R. Poryazova, G. Tononi, and S. L. Small, “Plastic changes following imitation-based speech and language therapy for aphasia: a high-density sleep eeg study,” _Neurorehabilitation and neural repair_ , vol. 28, no. 2, pp. 129–138, 2014\.
  * [17] R. Vergin, D. O’Shaughnessy, and A. Farhat, “Generalized mel frequency cepstral coefficients for large-vocabulary speaker-independent continuous-speech recognition,” _IEEE Transactions on speech and audio processing_ , vol. 7, no. 5, pp. 525–532, 1999.
  * [18] J. Chung, C. Gulcehre, K. Cho, and Y. Bengio, “Empirical evaluation of gated recurrent neural networks on sequence modeling,” _arXiv preprint arXiv:1412.3555_ , 2014.
  * [19] D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,” _arXiv preprint arXiv:1412.6980_ , 2014.
  * [20] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: a simple way to prevent neural networks from overfitting,” _The journal of machine learning research_ , vol. 15, no. 1, pp. 1929–1958, 2014.
  * [21] A. Graves, S. Fernández, F. Gomez, and J. Schmidhuber, “Connectionist temporal classification: labelling unsegmented sequence data with recurrent neural networks,” in _Proceedings of the 23rd international conference on Machine learning_. ACM, 2006, pp. 369–376.
  * [22] A. Graves and N. Jaitly, “Towards end-to-end speech recognition with recurrent neural networks,” in _International Conference on Machine Learning_ , 2014, pp. 1764–1772.
  * [23] S. Toshniwal, A. Kannan, C.-C. Chiu, Y. Wu, T. N. Sainath, and K. Livescu, “A comparison of techniques for language model integration in encoder-decoder speech recognition,” in _2018 IEEE Spoken Language Technology Workshop (SLT)_. IEEE, 2018, pp. 369–375.
  * [24] G. Krishna, C. Tran, M. Carnahan, and A. Tewfik, “Advancing speech recognition with no speech or with noisy speech,” in _2019 27th European Signal Processing Conference (EUSIPCO)_. IEEE, 2019.
  * [25] S. Mika, B. Schölkopf, A. J. Smola, K.-R. Müller, M. Scholz, and G. Rätsch, “Kernel pca and de-noising in feature spaces,” in _Advances in neural information processing systems_ , 1999, pp. 536–542.
  * [26] C. Goutte and E. Gaussier, “A probabilistic interpretation of precision, recall and f-score, with implication for evaluation,” in _European conference on information retrieval_. Springer, 2005, pp. 345–359.
  * [27] A. Flinker, A. Korzeniewska, A. Y. Shestyuk, P. J. Franaszczuk, N. F. Dronkers, R. T. Knight, and N. E. Crone, “Redefining the role of broca’s area in speech,” _Proceedings of the National Academy of Sciences_ , vol. 112, no. 9, pp. 2871–2875, 2015.
  * [28] J. Chartier, G. K. Anumanchipalli, K. Johnson, and E. F. Chang, “Encoding of articulatory kinematic trajectories in human speech sensorimotor cortex,” _Neuron_ , vol. 98, no. 5, pp. 1042–1054, 2018.
  * [29] S. Bai, J. Z. Kolter, and V. Koltun, “An empirical evaluation of generic convolutional and recurrent networks for sequence modeling,” _arXiv preprint arXiv:1803.01271_ , 2018.
  * [30] P. Griffiths and J. Needleman, “Statistical significance testing and p-values: Defending the indefensible? a discussion paper and position statement,” _International journal of nursing studies_ , vol. 99, p. 103384, 2019.


<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
11institutetext: Max Planck Institute for Physics, Föhringer Ring 6, 80805
München, Germany 22institutetext: Institute for Theoretical Physics, ETH,
CH-8093 Zürich, Switzerland 33institutetext: Department of Physics, University
of Zürich, CH-8057 Zürich, Switzerland 44institutetext: Institute for Particle
Physics Phenomenology, Durham University, Durham DH1 3LE, UK 55institutetext:
Department of Physics, Durham University, Durham DH1 3LE, UK 66institutetext:
Theoretical Physics Department, CERN, CH-1211 Geneva 23, Switzerland

# NNLO QCD predictions for Z-boson production in association with a charm jet
within the LHCb fiducial region

R. GauldMPI,rg  A. Gehrmann–De RidderETH,UZH,ag  E. W. N. GloverDUR1,DUR2,ng
A. HussCERN,ah  A. Rodriguez GarciaETH,ar  G. StagnittoUZH,gs

(MPP-2023-37, ZU-TH 12/23, IPPP/23/09, CERN-TH-2023-034)

###### Abstract

We compute next-to-next-to-leading order (NNLO) QCD corrections to neutral
vector boson production in association with a charm jet at the LHC. This
process is studied in the forward kinematics at $\sqrt{s}=13$ TeV, which may
provide valuable constraints on the intrinsic charm component of the proton. A
comparison is performed between fixed order and NLO predictions matched to a
parton shower showing mutual compatibility within the respective
uncertainties. NNLO corrections typically lead to a reduction of theoretical
uncertainties by a factor of two and the perturbative convergence is further
improved through the introduction of a theory-inspired constraint on the
transverse momentum of the vector boson plus jet system. A comparison between
these predictions with data will require an alignment of a flavour-tagging
procedure in theory and experiment that is infrared and collinear safe.

††journal: Eur. Phys. J. C

##  1 Introduction

The study of scattering processes that involve the direct production of
(heavy)-flavoured jets, i.e. those consistent with originating from charm
(${\mathrm{c}}$) or bottom (${\mathrm{b}}$) quarks, in association with a
leptonically decaying vector boson is essential for collider physics
phenomenology. They form a major background for several Standard Model (SM)
physics processes, including the production of a Higgs boson in association
with a gauge boson where the Higgs boson decays into heavy-flavoured jets
ATLAS:2018kot ; CMS:2018nsn ; ATLAS:2019yhn ; ATLAS:2020fcp ; ATLAS:2020jwz ,
as well as signals expected in models of physics beyond the SM (BSM)
CMS:2019zmd ; ATLAS:2021yij . Furthermore, they can provide unique information
on the distribution of flavoured partons inside the proton ATLAS:2014jkm ;
LHCb:2021stx ; CMS:2022bjk . Focussing on the process of ${\mathrm{Z}}$ plus
flavoured jet at the Large Hadron Collider (LHC), several measurements have
been performed by the ATLAS, CMS and the LHCb collaborations at 7 and 8 TeV
proton-proton collision energies ATLAS:2014rjv ; CMS:2014jqj ; LHCb:2014ydc ;
CMS:2017snu ; CMS:2016gmz ; ATLAS:2020juj . Recent studies at 13 TeV by the
CMS collaboration CMS:2020hmf  presented measurements of observables related
to the production of ${\mathrm{c}}$ and/or ${\mathrm{b}}$-quark jets in a
sample containing a ${\mathrm{Z}}$-boson produced in association with at least
one jet.

The production of a leptonically decaying ${\mathrm{Z}}$-boson in association
with a charm jet, and particularly at forward kinematics Boettcher:2015sqn
which is the focus of this work, could yield a unique probe of the charm
content of the proton Boettcher:2015sqn ; Bailas:2015jlc ; Lipatov:2016feu ,
provided that precise predictions and measurements of flavour-sensitive
${\mathrm{Z}}+c\text{-jet}$ observables are available and can be compared at a
similar level of accuracy. The LHCb collaboration has recently analysed events
containing a ${\mathrm{Z}}$-boson and a charm jet in the forward region of
phase space in proton-proton collisions LHCb:2021stx . These measurements
simultaneously provide direct access to the small- and large-$x$ regions of
the ${\mathrm{c}}$-quark parton distribution function (PDF) that is not well
explored by other experiments. Specifically, LHCb has presented results
LHCb:2021stx  for the ratio of production cross sections
$R^{c}_{j}=\sigma({\mathrm{Z}}+c\text{-jet})/\sigma({\mathrm{Z}}+\text{jet})$.
This ratio is measured differentially as a function of the rapidity of the
${\mathrm{Z}}$-boson $y^{{\mathrm{Z}}}$ in the range
$2.0<y^{{\mathrm{Z}}}<4.5$. The experimental result for the ratio $R^{c}_{j}$
has been compared with several SM predictions obtained at NLO QCD accuracy
interfaced with a parton shower (NLO+PS), each using different input PDF sets.
It is demonstrated that the most forward $y^{{\mathrm{Z}}}$ region is
particularly sensitive to the theoretical modelling of the charm quark PDF in
these sets, with the best agreement between theory and data obtained by
choosing a PDF set with a valence-like intrinsic (non-perturbative) charm
quark component. The presence or absence of this component is a long standing
theoretical issue Brodsky:1980pb ; Brodsky:2015fna . Recently, the NNPDF
collaboration has claimed evidence for an intrinsic charm quark component in
the proton Ball:2022qks , with a local significance at the 2.5$\sigma$ level
for momentum fractions in the region $0.3\lesssim x\lesssim 0.6$. By including
the LHCb data for $R^{c}_{j}$ in a reweighting of their fit, and adopting a
theory prediction based on NLO+PS, the local significance further increases to
about 3.0$\sigma$. Other PDF fitting collaborations have independently
investigated the possible presence of intrinsic charm in the proton
Hou:2017khm : for instance, a recent analysis by the CTEQ-TEA collaboration
Guzzi:2022rca  concludes that finding evidence for nonperturbative charm
continues to be elusive, by highlighting challenging aspects that must be
confronted in extracting nonperturbative charm in PDF fits.

State-of-the-art predictions for this kind of processes featuring the
associated production of a vector boson with one or more flavoured jets has
reached next-to-next-to-leading order (NNLO) accuracy in QCD calculations with
massless quarks Gauld:2020deh ; Czakon:2020coa ; Hartanto:2022qhh . Given the
importance of the LHCb data for PDF extractions, it is highly desirable to
have fixed-order predictions for the ${\mathrm{Z}}+c\text{-jet}$ process in
the forward region, in order to incorporate data for $R^{c}_{j}$ or other
flavour-sensitive observables in global PDF analyses based on collinear
factorisation. In this paper we focus on a description of
${\mathrm{Z}}+c\text{-jet}$ production within the LHCb fiducial region
following the approach of Gauld:2022lem  to define flavoured jets, and provide
a detailed comparison of (new) fixed-order predictions up to NNLO QCD accuracy
as well as those based on NLO+PS for a variety of differential distributions.

Despite our focus on the forward kinematics relevant to the LHCb detector, in
this work we refrain from performing a comparison to the available LHCb data
LHCb:2021stx . This is due to a significant contamination of the observable
$R^{c}_{j}$ measured by LHCb from Multiple Particle Interactions (MPI), which
should be removed/subtracted before considering this data in a (single parton
scattering) collinear PDF fit. Moreover, the experimental definition of jet
flavour in LHCb:2021stx  is not infrared and collinear (IRC) safe, rendering a
massless fixed-order calculation ill defined for the experimental set-up. With
the goal in mind to provide constraining information on the potential presence
of an intrinsic charm quark component within the proton, it is critical that
the definition of the presented data is IRC safe such that a massless
calculation (where collinear factorisation for the charm-quark PDF has been
performed) can be applied. We elaborate more on these issues in Appendix A and
B.

The structure of the paper is as follows. In Section 2, we present the main
ingredients entering the computation of observables associated with
${\mathrm{Z}}+c\text{-jet}$ production both at pure fixed-order perturbative
QCD as well as using the NLO+PS framework. We further comment on the proposal
of Gauld:2022lem , the flavour dressing algorithm, which allows flavour to be
assigned to anti-$k_{\mathrm{T}}$ jets in an IRC safe way. In Section 3, after
having described the numerical set-up and defined the scale variation
prescriptions, we present for the first time fixed-order predictions up to
next-to-next-to-leading order (NNLO) in QCD for several observables related to
the ${\mathrm{Z}}+c\text{-jet}$ process computed for the LHCb experimental
fiducial region at 13 TeV. We compare these fixed-order predictions with NLO
predictions matched to a parton shower at the parton level using the flavour
dressing procedure to define flavoured jets. We further investigate the impact
of an additional constraint on the transverse momentum of the vector boson
plus jet system in the computations. We shall find that the inclusion of this
theoretically motivated cut, brings NNLO and NLO+PS predictions closer and
improves the perturbative convergence of the fixed-order results for a large
fraction of the flavour-sensitive observables and for most of the kinematical
range studied. We also present predictions for the ratio $R^{c}_{j}$. In
Section 4, we summarise our findings and discuss the prospects of a direct
comparison between theory and data in the future. The Appendices contain a
discussion of the IRC safety of the current experimental definition (Appendix
A) and the role of MPI in the current experimental set-up (Appendix B).

##  2 Details of the calculation

In this paper, one of the main goals is to present fixed-order predictions
including QCD corrections up to ${\cal O}(\alpha_{s}^{3})$ for observables
related to the process
${\mathrm{p}}{\mathrm{\bar{p}}}\to{\ell}{\bar{\ell}}+c\text{-jet}+X$, yielding
a final state that contains a pair of charged leptons and at least one
identified charm jet.

The computation of higher order corrections to observables with identified
flavour poses several challenges as compared to flavour-blind cross sections:
it requires a complete flavour tracking of the particles in all subprocesses
(which are inputs to the flavour-dependent jet reconstruction/tagging), and
additionally to those appearing in all subtraction terms (for a calculation
based on subtraction). A flavour tracking procedure at parton level has been
pioneered for the computation of $b$-jet flavoured observables in
Gauld:2019yng  and Gauld:2020deh  and implemented within the NNLO  JET parton-
level generator that is used here. Within this framework, which employs the
antenna subtraction method to capture the infrared behaviour of matrix-
elements yielding fully differential cross section predictions at NNLO level,
the flavour tracking procedure has the crucial property that it can be applied
to any flavour-blind computation already present in the NNLO  JET code. For
example, the computation of ${\mathrm{Z}}+b\text{-jet}$ observables in
Gauld:2020deh  relied on the use of the existing flavour blind
${\mathrm{Z}}+\text{jet}$ computation presented in Gehrmann-DeRidder:2015wbt
and used the flavour-$k_{\mathrm{T}}$ algorithm to select $b$-jets. To compute
observables related to the ${\mathrm{Z}}+c\text{-jet}$ production process in
this work, we adopt a similar strategy by using the ${\mathrm{Z}}+\text{jet}$
computation including up to ${\cal O}(\alpha_{\mathrm{s}}^{3})$ corrections,
and then apply the flavour tracking procedure as was done for the
${\mathrm{Z}}+b\text{-jet}$ process.

As was the case in Gauld:2020deh , the prediction of scattering processes with
heavy-flavour jets can be further improved by exactly including the
contribution from a massive heavy-flavour quark at fixed order—i.e. the
resultant prediction is made in a general mass variable flavour number scheme.
We follow a similar procedure here and include mass corrections up to ${\cal
O}(\alpha_{s}^{2}m_{{\mathrm{c}}})$ in the fixed-order distributions which are
labelled as NLO and NNLO.

A further complication which is encountered in a calculation of the type
presented here—the calculation of QCD corrections to a scattering process with
flavoured jets based on massless quarks—is that an IRC safe definition of jet
flavour must be used. A first solution to this problem was introduced in
Banfi:2006hf , with the formulation of the IRC safe flavour-$k_{\mathrm{T}}$
algorithm. This algorithm features a $k_{\mathrm{T}}$-like clustering
sequence, and introduces a specific flavour-dependent clustering sequence to
achieve IRC flavour safety. However, since the algorithm requires the
knowledge of the flavour of all the particles in the event at each step of the
clustering, it is challenging to realise experimentally, and so far has not
been implemented in experimental analyses. Therefore, the kinematics and
flavour of the jets obtained with the flavour-$k_{\mathrm{T}}$ algorithm are
not compatible with those obtained in experiment. Various alternatives to the
use of flavour $k_{T}$ in theoretical predictions have recently been proposed
Caletti:2022glq ; Caletti:2022hnc ; Czakon:2022wam .

In the present analysis, we will adopt the flavour dressing algorithm of
Gauld:2022lem . This approach is particularly suitable as it enables to assign
flavour quantum numbers to a set of flavour blind jets, obtained with any jet
clustering algorithm. This allows us to apply it to jets reconstructed with
the anti-$k_{T}$ algorithm Cacciari:2008gp , i.e. the same algorithm which is
used in experiment to define the kinematics of the jets (although the flavour
assignment procedure is different, as detailed in Appendix A). Here, we stress
the fact that in the flavour dressing algorithm the flavour assignment is
entirely factorised from the initial jet reconstruction, hence the kinematics
of the jets is not affected. Such a key property is relevant in the present
context, since it ensures that for a ratio observable such as
$\sigma({\mathrm{Z}}+c\text{-jet})/\sigma({\mathrm{Z}}+\text{jet})$ both the
numerator and the denominator feature the same sample of anti-$k_{\mathrm{T}}$
jets.

A comparison of the fixed-order predictions as described above will also be
made to several NLO predictions matched with a parton shower. Those NLO+PS
predictions are obtained either with the MadGraph5_aMC@NLO (v. 2.7.3)
Alwall:2014hca  framework interfaced to Pythia8 (v. 8.243) Sjostrand:2014zea
(default $p_{\mathrm{T}}$-ordered parton shower) or Herwig7 (v. 7.2.2)
Bahr:2008pv ; Bellm:2015jjp ; Bellm:2019zci  (default angular-ordered parton
shower), and using the same flavour dressing procedure to define flavoured
jets as for the fixed-order predictions. To allow for a more direct
comparison, those predictions are obtained at the parton level where neither
the impact of MPI or hadronisation effects are included. A discussion on the
important role of MPI for the ${\mathrm{Z}}+c\text{-jet}$ process in the
forward region is provided in Appendix B. Hadronisation effects, on the other
hand, were found not to impact the considered observables in any significant
manner.

##  3 Numerical results

###  3.1 Numerical set-up and scale variation prescription

In this section, we review the calculational set-up as well as the kinematical
constraints imposed to obtain the fiducial cross sections for Z+$c\text{-jet}$
production. To select our final-state events, we focus on the forward region
with fiducial cuts mirroring those of the LHCb measurement LHCb:2021stx  at
$\sqrt{s}=13$ TeV.

In particular, the following fiducial cuts for jets and charged leptons are
applied: $20\leavevmode\nobreak\
\mathrm{GeV}<p_{\mathrm{T},j}<100\leavevmode\nobreak\ \mathrm{GeV}$,
$2.2<\eta_{j}<4.2$, $p_{\mathrm{T},\ell}>20\leavevmode\nobreak\ \mathrm{GeV}$,
$2.0<y_{\ell}<4.5$, $M_{\ell\bar{\ell}}\in[60,120]\leavevmode\nobreak\
\mathrm{GeV}$ and $\Delta R(j,\ell)>0.5$. The jets are reconstructed with the
anti-$k_{T}$ algorithm Cacciari:2008gp  with $R=0.5$. As discussed in Section
2, the selection of $c\text{-jets}$ is performed using the flavour dressing
procedure described in Gauld:2022lem . The algorithm proceeds in two stages
with internal parameters that control the overall flavour-tagging procedure: a
flavour clustering stage that employs a Soft-Drop-inspired criterion
Larkoski:2014wba  with parameters $z_{\rm cut}$, $R_{\rm cut}$, and $\beta$,
followed by a flavour dressing stage based on the flavour-$k_{T}$ distance
measure Banfi:2006hf  with a parameter $\alpha$. In the present calculation,
we set the parameters to their default values Gauld:2022lem : $z_{\rm
cut}=0.1$, $R_{\rm cut}=0.1$, $\beta=2$, $\alpha=2$. In addition, events are
only retained if the flavour-tagged $c\text{-jet}$ is the jet carrying the
largest transverse momentum of reconstructed jets passing the selection cuts.

We provide predictions for proton-proton collisions at $\sqrt{s}=13$ TeV and
use the PDF4LHC21 Monte Carlo PDF set PDF4LHCWorkingGroup:2022cjn , with
$\alpha_{s}(M_{\mathrm{Z}})=0.118$ and $n_{f}^{\rm max}=5$, where both the PDF
and $\alpha_{s}$ values are accessed via LHAPDF Buckley:2014ana . For the
electroweak input parameters, the results are obtained in the
$G_{\mu}$-scheme, using a complex mass scheme for the unstable internal
particles, and we adopt the following values for the input parameters:
$M_{{\mathrm{Z}}}^{\mathrm{os}}=91.1876\leavevmode\nobreak\ \mathrm{GeV}$,
$\Gamma_{{\mathrm{Z}}}^{\mathrm{os}}=2.4952\leavevmode\nobreak\ \mathrm{GeV}$,
$M_{{\mathrm{W}}}^{\mathrm{os}}=80.379\leavevmode\nobreak\ \mathrm{GeV}$,
$\Gamma_{{\mathrm{W}}}^{\mathrm{os}}=2.085\leavevmode\nobreak\ \mathrm{GeV},$
and $G_{\mu}=1.1663787\times 10^{-5}\leavevmode\nobreak\ \mathrm{GeV}^{-2}$.

For differential distributions, the impact of missing higher-order corrections
is assessed using the conventional 7-point scale variation prescription: the
values of factorisation ($\mu_{F}$) and renormalisation ($\mu_{R}$) scales are
varied independently by a factor of two around the central scale
$\mu_{0}\equiv E_{\mathrm{T},{\mathrm{Z}}}$, with the additional constraint
that $\frac{1}{2}\leq\mu_{F}/\mu_{R}\leq 2$.

When considering theoretical predictions for the ratio of distributions, we
estimate the uncertainties in an uncorrelated way between the numerator and
denominator i.e. by considering

$R^{c}_{j}(\mu_{R},\mu_{F};\mu^{\prime}_{R},\mu^{\prime}_{F})=\frac{\sigma^{{\mathrm{Z}}+c\text{-jet}}(\mu_{R},\mu_{F})}{\sigma^{{\mathrm{Z}}+\text{jet}}(\mu^{\prime}_{R},\mu^{\prime}_{F})}\,,$
(1)

providing a total of 31-points when dropping the extreme variations in any
pair of scales.

###  3.2 $Z+c$-jet distributions

We here present results for the ${\mathrm{Z}}+c\text{-jet}$ process at
$\sqrt{s}=13$ TeV and choose to focus on the following observables: the
leading flavoured jet transverse momentum $p_{\mathrm{T}}^{c\text{-jet}}$
(Figure 1), the leading flavoured jet pseudorapidity $\eta^{c\text{-jet}}$
(Figure 2), and the rapidity of the ${\mathrm{Z}}$-boson $y^{{\mathrm{Z}}}$,
reconstructed from the two final-state opposite-charge leptons (Figure 3).
Besides presenting results for the LHCb kinematical set-up as indicated in
Section 3.1, we also explore the impact of the introduction of a cut on the
transverse momentum of the ${\mathrm{Z}}+\text{jet}$ system:

$p_{\mathrm{T}}({\mathrm{Z}}+\text{jet})<p_{\mathrm{T}}^{\text{jet}}\,,$ (2)

with the leading jet in the acceptance region. The theoretical motivation
behind this cut is to discard those contributions where the flavoured jet is
not the jet with the largest transverse momentum in the event, i.e. cases
where the hardest jet was disregarded because it fell outside of the LHCb
acceptance. At Born level, the $p_{\mathrm{T}}$ of the
${\mathrm{Z}}+\text{jet}$ system vanishes, hence the cut in (2) limits the
hard QCD radiation outside the LHCb acceptance in a dynamical way.

For each of the figures presented in the remainder of this Section, the left
sides present results without the additional cut of (2), whereas the right
sides present results where the additional cut has been imposed, such that the
impact of the cut can be seen comparing the left and right hand sides of the
figures. To highlight the size and shape of the fixed-order results at each
perturbative order, and to best compare fixed-order with and without matching
to PS results, all figures illustrating the results in this Section are
composed of three panels: the top-panel shows the absolute predictions at
fixed-order (LO, NLO, NNLO) and for NLO+PS where PS is modelled by Pythia8;
the middle panel shows the ratio of NNLO to the fixed-order NLO result while
the lower panel shows the ratio to NNLO of the NLO+PS results where PS is
modeled by either Pythia8 or Herwig7. As noted in Section 2, fixed-order
predictions labelled as NLO and NNLO in all the figures include QCD
corrections up to ${\cal O}(\alpha_{s}^{2})$ (NLO) and ${\cal
O}(\alpha_{s}^{3})$ (NNLO) obtained with massless c-quarks in the
computations, and both additionally include the exact charm-quark mass
corrections up to ${\cal O}(\alpha_{s}^{2})$.

(a)

(b)

Figure 1:  Comparison of parton-level predictions for the leading flavoured
jet transverse momentum $p_{\mathrm{T}}^{c\text{-jet}}$ in the
${\mathrm{Z}}+c\text{-jet}$ process: fixed-order predictions at LO (green),
NLO (blue) and NNLO (red); NLO+PS predictions with Pythia8 (orange) or Herwig7
(purple) as parton showers. A dynamical cut on the transverse momentum of the
${\mathrm{Z}}+\text{jet}$ system is further applied in (1(b)).

(a)

(b)

Figure 2: Comparison of parton-level predictions for the leading flavoured jet
pseudo-rapidity $\eta^{c\text{-jet}}$ in the ${\mathrm{Z}}+c\text{-jet}$
process: fixed-order predictions at LO (green), NLO (blue) and NNLO (red);
NLO+PS predictions with Pythia8 (orange) or Herwig7 (purple) as parton
showers. A dynamical cut on the transverse momentum of the
${\mathrm{Z}}+\text{jet}$ system is further applied in (2(b)).

(a)

(b)

Figure 3: Comparison of parton-level predictions for the rapidity distribution
of the lepton pair $y^{Z}$ in the ${\mathrm{Z}}+c\text{-jet}$ process: fixed-
order predictions at LO (green), NLO (blue) and NNLO (red); NLO+PS predictions
with Pythia8 (orange) or Herwig7 (purple) as parton showers. A dynamical cut
on the transverse momentum of the ${\mathrm{Z}}+\text{jet}$ system is further
applied in (3(b)).

The theory predictions which _do not_ include the additional kinematic
constraint of Eq. (2) are shown in Figs. 1(a), 2(a), and 3(a). The NNLO
corrections are observed to be of the order (10–20)% as compared to the NLO
prediction, and typically outside of the scale variation band of the NLO
result. In particular, the shape of the distribution is modified in the
large-$\eta^{c\text{-jet}}$ and large-$y^{{\mathrm{Z}}}$ bins in Fig. 2(a) in
Fig. 3(a) respectively, and in the small-$p_{\mathrm{T}}^{c\text{-jet}}$
region in Fig. 1(a), featuring a larger cross section at NNLO. Those features
are made apparent in the second panel of those figures. In the third panel,
the NNLO and both NLO+PS predictions are shown normalised to the central NNLO
result. It is found that the NNLO result always lies between the two different
NLO+PS results. We observe that the NNLO result is more consistent with the
NLO+Herwig7 prediction for lower values of $p_{\mathrm{T}}^{c\text{-jet}}$,
$\eta^{c\text{-jet}}$, and $y^{{\mathrm{Z}}}$, and instead agrees better with
NLO+Pythia8 predictions at larger values. Overall, the NLO+Herwig7 prediction
seems to be rather similar to the NLO prediction, and the angular-ordered
parton shower does not appear to impact flavour-sensitive observables in a
significant way. Instead, the NLO+Pythia8 prediction does seem to capture some
of the NNLO higher order corrections: at
large-$p_{\mathrm{T}}^{c\text{-jet}}$,$\eta^{c\text{-jet}}$ and
$y^{{\mathrm{Z}}}$ values it tends to reproduce the shape of the fixed-order
NNLO corrections.

The theory predictions that _do_ include the additional kinematic constraint
of Eq. (2) are shown in Figs. 1(b), 2(b), and 3(b). We observe that the
constraint leads to a slight reduction of the fiducial cross section and
produces no significant change in the shape of the distributions with the
exception of the low-$p_{\mathrm{T}}^{c\text{-jet}}$ region. The LO, NLO, and
NNLO results display an improved mutual compatibility across all considered
observables, indicating a better perturbative convergence in the presence of
this kinematic constraint. Qualitatively, the comparison between the NNLO and
NLO+PS results in the third panel of these figures is similar to that of the
case without the kinematic constraint which was already discussed.

Overall, for all considered set-ups, we find that the NNLO corrections bring a
new level of precision for the considered ${\mathrm{Z}}+c\text{-jet}$
observables. As compared to the corresponding NLO(+PS) results, the
uncertainties of the NNLO predictions are reduced by a factor of two or more.
It is also reassuring that the NNLO corrections lead to predictions that tend
to lie between the two different NLO+PS predictions—the latter differ in the
treatment of ${\cal O}(\alpha_{s}^{3})$ terms, which begin at the NNLO level.
As the perturbative convergence of the predictions appears to be improved by
applying the kinematic constraint of Eq. (2), with only a small reduction in
the cross section, it seems to be well motivated when directly considering
${\mathrm{Z}}+c\text{-jet}$ observables.

###  3.3 The $\sigma({\mathrm{Z}}+\text{jet})$ cross-section and the ratio
$R^{c}_{j}$

In Section 3.2 we have presented IRC safe predictions for the rates of
${\mathrm{Z}}+c\text{-jet}$ production within the LHCb fiducial region.
Instead, in this subsection we will consider the ${\mathrm{Z}}+\text{jet}$
process (i.e. the flavour inclusive one) and then subsequently the ratio
observable
$R^{c}_{j}=\sigma({\mathrm{Z}}+c\text{-jet})/\sigma({\mathrm{Z}}+\text{jet})$.
As the experimental measurement of $R^{c}_{j}$ is performed differentially in
the rapidity of the ${\mathrm{Z}}$-boson $y^{{\mathrm{Z}}}$ LHCb:2021stx , our
theory predictions will also focus on this same quantity. We again note that
the we do not perform any comparison to the available data for reasons of
consistency, as detailed in Appendix A and B.

The theory predictions for ${\mathrm{Z}}+\text{jet}$ production are presented
in Fig. 4, with the same structure for the plots as shown in Section 3.2, i.e.
with the “no-cut” and “with-cut” cases shown on the left and right parts of
the figure respectively, and with three panels for each sub-figure. As
compared to the ${\mathrm{Z}}+c\text{-jet}$ predictions, heavy-flavour mass
corrections have not been included for these predictions. While this could be
achieved following the procedure outlined in Gauld:2021zmq , the numerical
impact of such corrections for the flavour inclusive process is negligible
(sub-percent).

The theory predictions without the extra kinematic cut of Eq. (2) are shown in
Fig. 4(a). The NNLO result is observed to be contained within the scale
variation band of the NLO one, except in the large-$y^{{\mathrm{Z}}}$ region,
where it also features a different behaviour compared to the NLO+PS results.
We also find very good agreement between the two NLO+PS results with Pythia8
and Herwig7.

The impact of applying the additional kinematic cut is shown in Fig. 4(b). The
cut leads to negative NNLO corrections in the (relatively) low
$y^{{\mathrm{Z}}}$ region, resulting in a NNLO prediction lying outside the
scale variation band of the NLO result, except in the large-$y^{{\mathrm{Z}}}$
region. Examination of the upper panel of Fig 4(b) shows that the cut has the
effect of moving all the curves closer to the LO result, with positive NLO
corrections and negative NNLO corrections. The inclusion of this cut thus
appears to degrade the perturbative stability for the flavour-inclusive set-
up.

(a)

(b)

Figure 4: Comparison of parton-level predictions for the rapidity distribution
of the lepton pair $y^{Z}$ for the unflavoured process
${\mathrm{Z}}+\text{jet}$: fixed-order predictions at LO (green), NLO (blue)
and NNLO (red); NLO+PS predictions with Pythia8 (orange) or Herwig7 (purple)
as parton showers. A dynamical cut on the transverse momentum of the
${\mathrm{Z}}+\text{jet}$ system is further applied in (4(b)).

We now consider the ratio observable $R^{c}_{j}$, differential in
$y^{{\mathrm{Z}}}$. The result is constructed using the same inputs which lead
to the distributions for ${\mathrm{Z}}+c\text{-jet}$ results in Fig. 3 and
${\mathrm{Z}}+\text{jet}$ in Fig. 4, but including the uncorrelated
uncertainty prescription defined in Eq. (1). The predictions for $R^{c}_{j}$
are displayed in Fig. 5, again the “no-cut” and “with-cut” cases are shown on
the left and right sides of the figure respectively. Focussing first on the
fixed-order results, the NNLO corrections are observed to be positive and of
the order (10–20)% with the largest value observed at large-$y^{{\mathrm{Z}}}$
values. That behaviour is observed for both the “no-cut” (left) and “with-cut”
(right) cases. Overall, the inclusion of the cut on the
${\mathrm{Z}}+\text{jet}$ system does not significantly impact the
perturbative behaviour of the fixed-order prediction for the ratio. This is a
consequence of the fact that the reduction of uncertainties in the numerator
for ${\mathrm{Z}}+c\text{-jet}$ is then compensated by increased uncertainties
in the denominator, when the cut is applied. However, by inspecting the lowest
panels of Fig. 5, we observe slightly better agreement between NNLO and the
two NLO+PS predictions when the cut has been applied.

Finally, as a result of the uncorrelated prescription for uncertainties, we
note that the relative theory uncertainties for the NNLO predictions of the
ratio are increased as compared to the individual predictions for
${\mathrm{Z}}+c\text{-jet}$ and ${\mathrm{Z}}+\text{jet}$. The sensitivity to
the input PDFs is also typically reduced for such an observable due to
correlations between PDF-dependence of the numerator and denominator. With the
aim of reducing the theory uncertainties due to missing higher corrections,
one could therefore consider to include absolute ${\mathrm{Z}}+c\text{-jet}$
cross-section data rather than that for the ratio $R^{c}_{j}$ in a collinear
PDF fit. Given however that several experimental uncertainties are correlated
between numerator and denominator (and therefore cancel in the ratio), and
that a treatment of the MPI contribution to the observable should also be
considered, overall it is not clear which observable is the most sensitive in
constraining PDFs.

(a)

(b)

Figure 5:  Comparison of parton-level predictions for the ratio
$R^{c}_{j}=\sigma^{{\mathrm{Z}}+c\text{-jet}}/\sigma^{{\mathrm{Z}}+\text{jet}}$
differential in the rapidity $y^{{\mathrm{Z}}}$ of the system of the final-
state leptons: fixed-order predictions at LO (green), NLO (blue) and NNLO
(red); NLO+PS predictions with Pythia8 (orange) or Herwig7 (purple) as parton
showers. A dynamical cut on the transverse momentum of the
${\mathrm{Z}}+\text{jet}$ system is further applied in (5(b)).

##  4 Conclusions

In this paper, we have studied the associated production of a
${\mathrm{Z}}$-boson with a charm-jet at the LHC at 13 TeV in the forward
region. We computed NNLO predictions at $\mathcal{O}(\alpha_{s}^{3})$ for a
set of differential observables related to the ${\mathrm{Z}}+c\text{-jet}$
process using the flavour dressing procedure to define charm-tagged jets. NNLO
corrections are found to be at the level of (10–20)%; they can impact the
shapes of distributions with the high-$y^{{\mathrm{Z}}}$ and
low-$p_{\mathrm{T}}^{c\text{-jet}}$ regions receiving enhanced corrections.
The residual uncertainties as estimated through scale variations are typically
$\pm 5$% or smaller, a factor two reduction compared to the respective NLO
uncertainty estimate. Additionally, comparisons to two different NLO+PS
predictions based on the Herwig7 (angular ordering) and Pythia8
($p_{\mathrm{T}}$ ordering) showers have been performed. The two predictions
can differ by up to 10% but remain mutually compatible within their respective
uncertainties that are NLO-like and thus at the level of $\pm 10$%. The NNLO
distributions are found to lie in between the two NLO+PS predictions, hinting
to an insensitivity to missing higher-order effects as modelled by the
showers. Moreover, we have found that a theory-inspired constraint on the
transverse momentum of the ${\mathrm{Z}}+\text{jet}$ system improves the
perturbative convergence of all considered distributions.

We further considered the ratio
$R^{c}_{j}=\sigma^{{\mathrm{Z}}+c\text{-jet}}/\sigma^{{\mathrm{Z}}+\text{jet}}$,
which has been measured by the LHCb experiment. In this context, the usage of
the flavour dressing procedure ensures the same kinematic reconstruction of
jets entering the numerator and denominator, thus allowing for a faithful
theoretical definition. The pattern of higher-order corrections mimics those
of the ${\mathrm{Z}}+c\text{-jet}$ process with enhanced NNLO corrections of
up to 20% in the high-$y^{{\mathrm{Z}}}$ region, albeit with larger
uncertainties due to the de-correlated scale prescription for the ratio.

A direct comparison to the available LHCb data was not performed due to IRC
unsafety issues and an unexpectedly large contamination from MPI; both of
which are discussed in the Appendices. In the future, a fair comparison
between experimental measurements and theory predictions will require a
detailed study of the experimental feasibility of the flavour dressing
algorithm (or other IRC safe variants). In the case a direct application of an
IRC-safe flavour definition is prohibitively challenging, it would be highly
desirable for experimental measurements to carry out an unfolding to a IRC-
safe definition of jet flavour. Only a joint effort of both communities,
theory and experimental, will enable to exploit in the best way the huge
amount of data that LHC will provide us in the next decades, better enabling
the use flavour signatures as a powerful window into short-distance
interactions from $\mathrm{GeV}$ to $\mathrm{TeV}$ energy scales.

###### Acknowledgements.

We thank Marco Zaro for his help with the usage of MadGraph5 aMC@NLO and
Christian Preuss for discussions about the treatment of heavy quarks in
Pythia8. RG and AH thank Philip Ilten for several discussions related to the
LHCb, and more generally on the experimental definition of flavoured jets. GS
thanks the Institute for Theoretical Physics at ETH for hospitality during the
course of this work. This research was supported in part by the UK Science and
Technology Facilities Council under contract ST/T001011/1 and by the Swiss
National Science Foundation (SNF) under contracts 200021-197130 and
200020-204200. Numerical simulations were facilitated by the High Performance
Computing group at ETH Zürich and the Swiss National Supercomputing Centre
(CSCS) under project ID ETH5f.

##  Appendix A IRC safety

According to Ref. LHCb:2021stx , as part of the LHCb measurement, charm jets
are defined in the following way. First, jets are reconstructed using the
anti-$k_{T}$ algorithm Cacciari:2008gp  with $R=0.5$, and the leading jet
passing the fiducial selection is considered (see also Section 3.1 for the
definition of the fiducial selection). From this, the leading jet is
considered to be a charm jet (at truth/unfolded level) if it additionally
satisfies the criterion: $p_{\mathrm{T},c{\leavevmode\nobreak\ \rm
hadron}}>5\leavevmode\nobreak\ \mathrm{GeV}$, and $\Delta
R(j,c{\leavevmode\nobreak\ \rm hadron})<0.5$. That is to say, the jet is
considered to be charm tagged if there is at least one $c{\leavevmode\nobreak\
\rm hadron}$ satisfying these selections.

Such an experimental definition of jet flavour is collinear unsafe. This is
not an LHCb specific issue (or even particular to this measurement), but is
also relevant in one way or another for the definitions of jet flavour
commonly adopted by LHC collaborations LHCb:2015tna ; ATLAS:2015thz ;
CMS:2017wtu ; LHCb:2021dlw . We briefly explain why the observable definition
taken in Ref. LHCb:2021stx  is IRC unsafe, then discuss the implications of
this for data interpretation.

(a)

(b)

(c)

Figure 6: Examples of configurations which lead to an IRC sensitivity: (6(a))
a double (even) tagged jet; (6(b)) a tag with a $p_{\mathrm{T},q}^{\rm min}$
requirement; (6(c)) a jet tag with a soft sensitivity (in the absence of a
$p_{\mathrm{T},q}^{\rm min}$ requirement).

In Fig. 6 we depict three specific kinematic configurations which lead to
different sources of IRC sensitivity (i.e. those which are IRC unsafe). In
each case, the configurations correspond to those encountered in a parton-
level fixed-order prediction for the ${\mathrm{Z}}+c\text{-jet}$ process. The
flavoured quarks are depicted as red lines, i.e. representing charm quarks for
the process under consideration here.

  1. (a)

The first configuration depicts the production of the lepton pair, recoiling
against a hard $q\bar{q}$ pair produced in a collinear configuration (i.e. at
least one, or both, of the quarks are hard but $p_{q}\cdot p_{\bar{q}}\to 0$).
When the anti-$k_{\mathrm{T}}$ algorithm is applied, the $q$ and $\bar{q}$ are
reconstructed inside the same jet. According to the prescription to assign jet
flavour, as this jet contains at least one quark it will be assigned a quark
flavour tag (e.g. charm for $q=c$). This introduces a collinear sensitivity as
it is distinguished from the case where the hard gluon does not split into a
collinear $q\bar{q}$ pair. In that case, the jet (composed of a single hard
gluon) clearly carries zero quark flavour. This could be overcome by
accounting for the total quark flavour in the jet, such as assigning quantum
numbers $q(\bar{q})=+(-)1$ and summing them to obtain the net flavour
(alternatively one can consider flavoured jets as those with an overall odd
number of $q$ and $\bar{q}$).

  2. (b)

The second configuration depicts the production of the lepton pair, recoiling
against a hard $qg$ pair in a collinear configuration (i.e. at least one, or
both, of the quark and gluon is hard but $p_{q}\cdot p_{g}\to 0$). Again, when
the anti-$k_{\mathrm{T}}$ algorithm is applied, both $q$ and $g$ are
reconstructed inside the same jet. As the tagging prescription requires the
presence of a $c$-hadron with $p_{\mathrm{T},c}>5\leavevmode\nobreak\
\mathrm{GeV}$, it is possible that the outgoing quark does not satisfy this
criterion (depending on the momentum sharing with the gluon). This introduces
a collinear sensitivity as the $p_{\mathrm{T},c}$ requirement may distinguish
from the case where the hard initial quark does not split to the collinear
$qg$ pair.

  3. (c)

The collinear sensitivity discussed above can be overcome by removing the
$p_{\mathrm{T},c}^{\rm min}$ requirement. However, this would introduce a new
problem which is depicted in the third configuration. In this case, the
culprit is a soft gluon which subsequently splits to a $q\bar{q}$ pair at wide
angles. It is possible that one of the quarks ($p_{q}$ in the figure) is
produced close in $\Delta R$ to a hard parton ($p_{j}$), i.e. $\Delta
R(j,q)<0.5$. This would introduces a soft sensitivity as the flavour of the
jet would be altered by the presence of the soft quark.

From a purely experimental point of view, it is clear that the current
definition of tagging heavy-flavour jets is a sensible one. Identifying those
jets with multiple tags (as oppose to at least one) requires to more carefully
account for the experimental (in)efficiency and mistag rates. It is also
extremely to difficult distinguish between the signature of one or two
collinear heavy-flavour objects (e.g. a bunch of displaced tracks appearing to
originate from a single displaced vertex). Furthermore, removing the
$p_{\mathrm{T},c}^{\rm min}$ would mean accounting for a region where it is
experimentally challenging to identify displaced vertices. However, this
choice has serious ramifications for the theory predictions, and importantly
for the interpretation of the data.

Theoretical predictions of charm jet observables which are not IRC safe are
logarithmically sensitive to the mass of the charm quark $m_{c}$. The
corresponding fixed-order predictions for such observables therefore include
corrections which depend logarithmically on the charm quark mass. If the
observable/process under consideration involves energy scales which are large
compared to the quark mass (e.g. the transverse energy/momentum of a jet or a
boson), the logarithmic corrections become large (due to the separation of
scales) and thus limit the theory precision/perturbative stability. The
$m_{c}\to 0$ limit of such predictions is not well defined (it is divergent),
meaning that a calculation based on massless quarks of such observables does
not exist. The implications of this are that fixed-order predictions must be
performed in a scheme where the charm quark is massive, i.e. in a fixed-
flavour-number scheme with $n_{f}^{\rm max}=3$, where mass factorisation is
not performed for the charm quark, and it is decoupled from the running of
$\alpha_{\mathrm{s}}$. Note that the perturbative charm-quark PDF does still
exist in the massive scheme (where a logarithmic sensitivity to the charm-
quark mass exists, see for example Laenen:1992zk ; Riemersma:1994hv ;
Harris:1995tu ; Buza:1996wv ; Laenen:1998kp ; Kawamura:2012cr ). Practically,
it is generated numerically after integration over the phase-space of the
massive quark during the calculation.

The requirement that a fixed-order prediction must be carried out with a
massive calculation is problematic for observables which are designed to be
sensitive to the nature of the charm quark PDF. Such observables contain a
logarithmic sensitivity on the charm-quark mass as a result of the IRC unsafe
configurations which were highlighted above (which limit the theory
precision/predictability), at the same time the observables are (by design)
directly sensitive to the large logarithmic corrections associated to the
perturbative charm quark PDF which have not been resummed, and the non-
perturbative component of the charm quark PDF (which is aiming to be probed)
does not exist.

Instead, the definition of jet flavour used in this work, given its IRC
safety, removes all logarithmic sensitivity on the charm quark mass which
results from the jet flavour assignment procedure. This also allows for the
application of a massless calculation, where mass factorisation for the charm
quark is performed and where remaining logarithmic sensitivity to the quark
mass is resummed.

##  Appendix B Multiple Particle Interactions

During the high-energy scattering of two protons, there is a probability for
multiple hard-interactions to occur (i.e. more than one).

(a)

(b)

(c)

Figure 7: Effect of MPI contributions on the ${\mathrm{Z}}$ rapidity
distribution $y^{\mathrm{Z}}$ in the ${\mathrm{Z}}+c\text{-jet}$ process
(7(a)), in the ${\mathrm{Z}}+\text{jet}$ process (7(b)) and in the ratio of
the two (7(c)). NLO+PS predictions are obtained with Pythia8 (orange) or
Herwig7 (purple) as parton showers. In the upper panels predictions including
(excluding) MPI contributions are depicted in darker (lighter) colours. The
lower panels show the ratios of curves with and without MPI effects.

For the LHCb kinematics defined at the beginning of Section 3.1, and also
applying the (IRC unsafe) definition of jet flavour as in LHCb:2021stx , we
observe a large contribution to the production of a ${\mathrm{Z}}$ boson and a
$c\text{-jet}$ from MPI. In Fig. 7(a) we show the cross-section for
${\mathrm{Z}}+c\text{-jet}$ production after fiducial cuts, which is plotted
differentially with respect to the ${\mathrm{Z}}$-rapidity $y^{{\mathrm{Z}}}$.
The predictions are provided at NLO+PS accuracy for ${\mathrm{Z}}+1j$ events
generated with MadGraph5_aMC@NLO interfaced with Pythia8 and Herwig7, where
the role of MPI is subsequently modelled by the two different Monte Carlo
generators. We show the predictions obtained when including/excluding the MPI
contributions, which lead to a large (and rapidity dependent) effect on the
resultant distribution. In the central rapidity region, the effect is of order
10%, increasing up to 20% at larger rapidities, and the effect is very similar
within Pythia8 and Herwig7. In Fig. 7(b) we further show the same plot for the
unflavoured process ${\mathrm{Z}}+\text{jet}$. We note a constant shift of the
$y^{{\mathrm{Z}}}$ distribution of order 10% due to MPI effects. The 10%–20%
effects observed in Fig. 7(a) for ${\mathrm{Z}}+c\text{-jet}$ could be
explained with an interplay of a constant 10% effect (impacting both
${\mathrm{Z}}+c\text{-jet}$ and ${\mathrm{Z}}+\text{jet}$ processes) and an
additional 10% flavour-dependent effect. Finally, we consider the ratio
$R^{c}_{j}$, even if its behaviour w.r.t. MPI effects can be straightforwardly
inferred from Fig. 7(a) and Fig. 7(b). In particular, we observe a partial
cancellation of MPI effects, leading to a difference between predictions which
is negligible at $y^{\mathrm{Z}}\sim 2.0$ and that increases monotonically up
to 10% at larger rapidities. Hence, given such a partial cancellation in the
ratio, the inclusion of MPI effects seems to be particularly relevant for
$3.5\lesssim y^{\mathrm{Z}}\lesssim 4.5$. Such a region is hugely important
for the determination of the intrinsic charm in the proton, given that the
LHCb measurement in the bin $y^{\mathrm{Z}}\in[3.5,4.5]$ is highly correlated
to the charm PDF in the region of the valence peak $x\simeq 0.45$ Ball:2022qks
.

As these MPI contributions arise from the overlap of multiple hard
interactions, it is not taken into account by a theoretical description based
solely on single parton scattering contributions (i.e. those which form the
basis of a collinear PDF fit). If the data is to be considered for such a fit,
it would be necessary to first remove/subtract the MPI component. Due to the
complicated (overlapping) nature of multiple hard interactions, the removal of
this component may rely on theoretical modelling. As the theoretical
description of inclusive charm-quark production with the LHCb acceptance
suffers from uncertainties far in excess of $50\%$ (see for example, Fig. 6 of
Ref. Gauld:2015yia  for $p_{\mathrm{T}}^{D}\sim 5\leavevmode\nobreak\
\mathrm{GeV}$), the modelling of this contribution (and the associated
uncertainty) will require quite some care. We note that such an uncertainty
(optimistically, 50% of the MPI contribution as shown in Fig. 7) would already
be in excess of the systematic uncertainty quoted for the ratio observable in
LHCb:2021stx . It may be possible to overcome this issue if a data-driven
based approach can be achieved.

## References

  * (1) ATLAS Collaboration, M. Aaboud et al., Observation of $H\rightarrow b\bar{b}$ decays and $VH$ production with the ATLAS detector, Phys. Lett. B 786 (2018) 59–86, [arXiv:1808.08238].
  * (2) CMS Collaboration, A. M. Sirunyan et al., Observation of Higgs boson decay to bottom quarks, Phys. Rev. Lett. 121 (2018), no. 12 121801, [arXiv:1808.08242].
  * (3) ATLAS Collaboration, M. Aaboud et al., Measurement of VH, $\mathrm{H}\to\mathrm{b}\overline{\mathrm{b}}$ production as a function of the vector-boson transverse momentum in 13 TeV pp collisions with the ATLAS detector, JHEP 05 (2019) 141, [arXiv:1903.04618].
  * (4) ATLAS Collaboration, G. Aad et al., Measurements of $WH$ and $ZH$ production in the $H\rightarrow b\bar{b}$ decay channel in $pp$ collisions at 13 TeV with the ATLAS detector, Eur. Phys. J. C 81 (2021), no. 2 178, [arXiv:2007.02873].
  * (5) ATLAS Collaboration, G. Aad et al., Measurement of the associated production of a Higgs boson decaying into $b$-quarks with a vector boson at high transverse momentum in $pp$ collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, Phys. Lett. B 816 (2021) 136204, [arXiv:2008.02508].
  * (6) CMS Collaboration, A. M. Sirunyan et al., Search for supersymmetry in proton-proton collisions at 13 TeV in final states with jets and missing transverse momentum, JHEP 10 (2019) 244, [arXiv:1908.04722].
  * (7) ATLAS Collaboration, G. Aad et al., Search for new phenomena in final states with $b$-jets and missing transverse momentum in $\sqrt{s}=13$ TeV $pp$ collisions with the ATLAS detector, JHEP 05 (2021) 093, [arXiv:2101.12527].
  * (8) ATLAS Collaboration, G. Aad et al., Measurement of the production of a $W$ boson in association with a charm quark in $pp$ collisions at $\sqrt{s}=$ 7 TeV with the ATLAS detector, JHEP 05 (2014) 068, [arXiv:1402.6263].
  * (9) LHCb Collaboration, R. Aaij et al., Study of Z Bosons Produced in Association with Charm in the Forward Region, Phys. Rev. Lett. 128 (2022), no. 8 082001, [arXiv:2109.08084].
  * (10) CMS Collaboration, Measurement of the production cross section of a W boson in association with a charm quark in proton-proton collisions at $\sqrt{s}=13\leavevmode\nobreak\ \mathrm{TeV}$, Tech. Rep. CMS-PAS-SMP-21-005, CERN, Geneva, 2022\.
  * (11) ATLAS Collaboration, G. Aad et al., Measurement of differential production cross-sections for a $Z$ boson in association with $b$-jets in 7 TeV proton-proton collisions with the ATLAS detector, JHEP 10 (2014) 141, [arXiv:1407.3643].
  * (12) CMS Collaboration, S. Chatrchyan et al., Measurement of the production cross sections for a Z boson and one or more b jets in pp collisions at sqrt(s) = 7 TeV, JHEP 06 (2014) 120, [arXiv:1402.1521].
  * (13) LHCb Collaboration, R. Aaij et al., Measurement of the Z+b-jet cross-section in pp collisions at $\sqrt{s}$ = 7 TeV in the forward region, JHEP 01 (2015) 064, [arXiv:1411.1264].
  * (14) CMS Collaboration, A. M. Sirunyan et al., Measurement of associated Z + charm production in proton-proton collisions at $\sqrt{s}=$ 8 TeV, Eur. Phys. J. C 78 (2018), no. 4 287, [arXiv:1711.02143].
  * (15) CMS Collaboration, V. Khachatryan et al., Measurements of the associated production of a Z boson and b jets in pp collisions at ${\sqrt{s}}=8\,\text{TeV}$, Eur. Phys. J. C 77 (2017), no. 11 751, [arXiv:1611.06507].
  * (16) ATLAS Collaboration, G. Aad et al., Measurements of the production cross-section for a $Z$ boson in association with $b$-jets in proton-proton collisions at $\sqrt{s}=13$ TeV with the ATLAS detector, JHEP 07 (2020) 044, [arXiv:2003.11960].
  * (17) CMS Collaboration, A. M. Sirunyan et al., Measurement of the associated production of a $Z$ boson with charm or bottom quark jets in proton-proton collisions at $\sqrt{s}$=13 TeV, Phys. Rev. D 102 (2020), no. 3 032007, [arXiv:2001.06899].
  * (18) T. Boettcher, P. Ilten, and M. Williams, Direct probe of the intrinsic charm content of the proton, Phys. Rev. D 93 (2016), no. 7 074008, [arXiv:1512.06666].
  * (19) G. Bailas and V. P. Goncalves, Phenomenological implications of the intrinsic charm in the $Z$ boson production at the LHC, Eur. Phys. J. C 76 (2016), no. 3 105, [arXiv:1512.06007].
  * (20) A. V. Lipatov, G. I. Lykasov, Y. Y. Stepanenko, and V. A. Bednyakov, Probing proton intrinsic charm in photon or Z boson production accompanied by heavy jets at the LHC, Phys. Rev. D 94 (2016), no. 5 053011, [arXiv:1606.04882].
  * (21) S. J. Brodsky, P. Hoyer, C. Peterson, and N. Sakai, The Intrinsic Charm of the Proton, Phys. Lett. B 93 (1980) 451–455.
  * (22) S. J. Brodsky, A. Kusina, F. Lyonnet, I. Schienbein, H. Spiesberger, and R. Vogt, A review of the intrinsic heavy quark content of the nucleon, Adv. High Energy Phys. 2015 (2015) 231547, [arXiv:1504.06287].
  * (23) NNPDF Collaboration, R. D. Ball, A. Candido, J. Cruz-Martinez, S. Forte, T. Giani, F. Hekhorn, K. Kudashkin, G. Magni, and J. Rojo, Evidence for intrinsic charm quarks in the proton, Nature 608 (2022), no. 7923 483–487, [arXiv:2208.08372].
  * (24) T.-J. Hou, S. Dulat, J. Gao, M. Guzzi, J. Huston, P. Nadolsky, C. Schmidt, J. Winter, K. Xie, and C. P. Yuan, CT14 Intrinsic Charm Parton Distribution Functions from CTEQ-TEA Global Analysis, JHEP 02 (2018) 059, [arXiv:1707.00657].
  * (25) M. Guzzi, T. J. Hobbs, K. Xie, J. Huston, P. Nadolsky, and C. P. Yuan, The persistent nonperturbative charm enigma, arXiv:2211.01387.
  * (26) R. Gauld, A. Gehrmann-De Ridder, E. W. N. Glover, A. Huss, and I. Majer, Predictions for $Z$ -Boson Production in Association with a $b$-Jet at $\mathcal{O}(\alpha_{s}^{3})$, Phys. Rev. Lett. 125 (2020), no. 22 222002, [arXiv:2005.03016].
  * (27) M. Czakon, A. Mitov, M. Pellen, and R. Poncelet, NNLO QCD predictions for W+c-jet production at the LHC, JHEP 06 (2021) 100, [arXiv:2011.01011].
  * (28) H. B. Hartanto, R. Poncelet, A. Popescu, and S. Zoia, Next-to-next-to-leading order QCD corrections to Wbb¯ production at the LHC, Phys. Rev. D 106 (2022), no. 7 074016, [arXiv:2205.01687].
  * (29) R. Gauld, A. Huss, and G. Stagnitto, A dress of flavour to suit any jet, arXiv:2208.11138.
  * (30) R. Gauld, A. Gehrmann-De Ridder, E. W. N. Glover, A. Huss, and I. Majer, Associated production of a Higgs boson decaying into bottom quarks and a weak vector boson decaying leptonically at NNLO in QCD, JHEP 10 (2019) 002, [arXiv:1907.05836].
  * (31) A. Gehrmann-De Ridder, T. Gehrmann, E. W. N. Glover, A. Huss, and T. A. Morgan, Precise QCD predictions for the production of a Z boson in association with a hadronic jet, Phys. Rev. Lett. 117 (2016), no. 2 022001, [arXiv:1507.02850].
  * (32) A. Banfi, G. P. Salam, and G. Zanderighi, Infrared safe definition of jet flavor, Eur. Phys. J. C 47 (2006) 113–124, [hep-ph/0601139].
  * (33) S. Caletti, A. J. Larkoski, S. Marzani, and D. Reichelt, A fragmentation approach to jet flavor, JHEP 10 (2022) 158, [arXiv:2205.01117].
  * (34) S. Caletti, A. J. Larkoski, S. Marzani, and D. Reichelt, Practical jet flavour through NNLO, Eur. Phys. J. C 82 (2022), no. 7 632, [arXiv:2205.01109].
  * (35) M. Czakon, A. Mitov, and R. Poncelet, Infrared-safe flavoured anti-$k_{T}$ jets, arXiv:2205.11879.
  * (36) M. Cacciari, G. P. Salam, and G. Soyez, The anti-$k_{t}$ jet clustering algorithm, JHEP 04 (2008) 063, [arXiv:0802.1189].
  * (37) J. Alwall, R. Frederix, S. Frixione, V. Hirschi, F. Maltoni, O. Mattelaer, H. S. Shao, T. Stelzer, P. Torrielli, and M. Zaro, The automated computation of tree-level and next-to-leading order differential cross sections, and their matching to parton shower simulations, JHEP 07 (2014) 079, [arXiv:1405.0301].
  * (38) T. Sjöstrand, S. Ask, J. R. Christiansen, R. Corke, N. Desai, P. Ilten, S. Mrenna, S. Prestel, C. O. Rasmussen, and P. Z. Skands, An introduction to PYTHIA 8.2, Comput. Phys. Commun. 191 (2015) 159–177, [arXiv:1410.3012].
  * (39) M. Bahr et al., Herwig++ Physics and Manual, Eur. Phys. J. C 58 (2008) 639–707, [arXiv:0803.0883].
  * (40) J. Bellm et al., Herwig 7.0/Herwig++ 3.0 release note, Eur. Phys. J. C 76 (2016), no. 4 196, [arXiv:1512.01178].
  * (41) J. Bellm et al., Herwig 7.2 release note, Eur. Phys. J. C 80 (2020), no. 5 452, [arXiv:1912.06509].
  * (42) A. J. Larkoski, S. Marzani, G. Soyez, and J. Thaler, Soft Drop, JHEP 05 (2014) 146, [arXiv:1402.2657].
  * (43) PDF4LHC Working Group Collaboration, R. D. Ball et al., The PDF4LHC21 combination of global PDF fits for the LHC Run III, J. Phys. G 49 (2022), no. 8 080501, [arXiv:2203.05506].
  * (44) A. Buckley, J. Ferrando, S. Lloyd, K. Nordström, B. Page, M. Rüfenacht, M. Schönherr, and G. Watt, LHAPDF6: parton density access in the LHC precision era, Eur. Phys. J. C75 (2015) 132, [arXiv:1412.7420].
  * (45) R. Gauld, A massive variable flavour number scheme for the Drell-Yan process, SciPost Phys. 12 (2022), no. 1 024, [arXiv:2107.01226].
  * (46) LHCb Collaboration, R. Aaij et al., Identification of beauty and charm quark jets at LHCb, JINST 10 (2015), no. 06 P06013, [arXiv:1504.07670].
  * (47) ATLAS Collaboration, G. Aad et al., Performance of $b$-Jet Identification in the ATLAS Experiment, JINST 11 (2016), no. 04 P04008, [arXiv:1512.01094].
  * (48) CMS Collaboration, A. M. Sirunyan et al., Identification of heavy-flavour jets with the CMS detector in pp collisions at 13 TeV, JINST 13 (2018), no. 05 P05011, [arXiv:1712.07158].
  * (49) LHCb Collaboration, R. Aaij et al., Identification of charm jets at LHCb, JINST 17 (2022), no. 02 P02028, [arXiv:2112.08435].
  * (50) E. Laenen, S. Riemersma, J. Smith, and W. L. van Neerven, Complete O (alpha-s) corrections to heavy flavor structure functions in electroproduction, Nucl. Phys. B 392 (1993) 162–228.
  * (51) S. Riemersma, J. Smith, and W. L. van Neerven, Rates for inclusive deep inelastic electroproduction of charm quarks at HERA, Phys. Lett. B 347 (1995) 143–151, [hep-ph/9411431].
  * (52) B. W. Harris and J. Smith, Heavy quark correlations in deep inelastic electroproduction, Nucl. Phys. B 452 (1995) 109–160, [hep-ph/9503484].
  * (53) M. Buza, Y. Matiounine, J. Smith, and W. L. van Neerven, Charm electroproduction viewed in the variable flavor number scheme versus fixed order perturbation theory, Eur. Phys. J. C 1 (1998) 301–320, [hep-ph/9612398].
  * (54) E. Laenen and S.-O. Moch, Soft gluon resummation for heavy quark electroproduction, Phys. Rev. D 59 (1999) 034027, [hep-ph/9809550].
  * (55) H. Kawamura, N. A. Lo Presti, S. Moch, and A. Vogt, On the next-to-next-to-leading order QCD corrections to heavy-quark production in deep-inelastic scattering, Nucl. Phys. B 864 (2012) 399–468, [arXiv:1205.5727].
  * (56) R. Gauld, J. Rojo, L. Rottoli, and J. Talbert, Charm production in the forward region: constraints on the small-x gluon and backgrounds for neutrino astronomy, JHEP 11 (2015) 009, [arXiv:1506.08025].