capbtabboxtable[][]

# Understanding Cognitive Fatigue from fMRI Scans with Self-supervised
Learning

Ashish Jaiswal
The University of Texas at Arlington
Arlington, TX 76013
<EMAIL_ADDRESS>
Ashwin Ramesh Babu
The University of Texas at Arlington
Arlington, TX 76013
<EMAIL_ADDRESS>
Mohammad Zaki Zadeh
The University of Texas at Arlington
Arlington, TX 76013
<EMAIL_ADDRESS>
Fillia Makedon
The University of Texas at Arlington
Arlington, TX 76013
<EMAIL_ADDRESS>
Glenn Wylie
Kessler Foundation
East Hanover, New Jersey, 07936
<EMAIL_ADDRESS>

###### Abstract

Functional magnetic resonance imaging (fMRI) is a neuroimaging technique that
records neural activations in the brain by capturing the blood oxygen level in
different regions based on the task performed by a subject. Given fMRI data,
the problem of predicting the state of cognitive fatigue in a person has not
been investigated to its full extent. This paper proposes tackling this issue
as a multi-class classification problem by dividing the state of cognitive
fatigue into six different levels, ranging from no-fatigue to extreme fatigue
conditions. We built a spatio-temporal model that uses convolutional neural
networks (CNN) for spatial feature extraction and a long short-term memory
(LSTM) network for temporal modeling of 4D fMRI scans. We also applied a self-
supervised method called MoCo (Momentum Contrast) to pre-train our model on a
public dataset BOLD5000 and fine-tuned it on our labeled dataset to predict
cognitive fatigue. Our novel dataset contains fMRI scans from Traumatic Brain
Injury (TBI) patients and healthy controls (HCs) while performing a series of
N-back cognitive tasks. This method establishes a state-of-the-art technique
to analyze cognitive fatigue from fMRI data and beats previous approaches to
solve this problem.

_Keywords_ fMRI $\cdot$ brain imaging $\cdot$ deep learning $\cdot$ self-
supervised learning $\cdot$ contrastive learning $\cdot$ cognitive fatigue

##  1 Introduction

Functional magnetic resonance imaging (fMRI) measures slight changes in blood
flow that occur with activity in different brain regions. This imaging
technique is completely safe and non-intrusive to the human brain. It is used
to identify parts of the brain that handle critical functions and evaluate the
effects of conditions such as stroke and other diseases. Some abnormalities
can only be found with fMRI scans as it provides detailed access to activity
patterns in a human brain. Traumatic Brain Injury (TBI) is one of the most
prevalent causes of neurological disorders in the US Faul et al. (2010). It is
a condition that has been shown to affect working memory Christodoulou et al.
(2001), and cognitive fatigue Kohl et al. (2009). In this work, we focus on
understanding cognitive fatigue that results from performing standardized
cognitive tasks as it is one of the primary indicators of moderate-to-severe
TBI.

Cognitive Fatigue (CF) is a subjective lack of mental energy that is perceived
by an individual which interferes with everyday activities DeLuca et al.
(2008). It is a common condition among people suffering from moderate to
severe brain injury. There are multiple approaches to assess CF through
various cognitive tasks and assessment tests by using objective measures such
as response time (RT), error rate (ER) DeLuca et al. (2008). However, these
measures have certain limitations and do not correlate well with the self-
reported scores during the tasks Wylie et al. (2017). The inability to relate
objective measures to self-reported cognitive fatigue led us to study the
blood-oxygen-level-dependent (BOLD) signal associated with neural activation
changes. The increased BOLD activation represents an increase in the cognitive
work that individuals with TBI must expend to perform a certain task Wylie et
al. (2017).

With the advancements in deep learning techniques that can efficiently extract
meaningful information from images/videos, we attempted to build a model that
predicts self-reported cognitive fatigue scores based on neural activations
captured through the fMRI scans. The main contribution of this work can be
listed as follows:

  * •

To the best of our knowledge, there has not been enough work done to analyse
cognitive fatigue from fMRI data.

  * •

We worked with a novel dataset that consists of fMRI scans from both TBI
patients and healthy subjects.

  * •

A self-supervised model finetuned on labeled data to understand and predict
the correlation between neural activation and self-reported fatigue scores.

The rest of the paper is structured as follows: we start by discussing some of
the previous works that examine cognitive fatigue from fMRI data and results
that use deep learning for brain imaging. Then we explain the data collection,
the preprocessing stage, and the system architecture. Finally, we show our
experiments and results, followed by our conclusion and future directions at
the end.

##  2 Related Work

Previous researches have demonstrated that in people with Traumatic Brain
Injury (TBI), the caudate nucleus of the basal ganglia shows a distinct
pattern of activation over time than in healthy controls Kohl et al. (2009).
This finding was consistent with Chauduri and Behan’s fatigue model Chaudhuri
and Behan (2000), in which the basal ganglia played a crucial part in fatigue
experience. On the other hand, Kohl et al. Kohl et al. (2009) inferred the
presence of fatigue based on the pattern of brain activations across time. The
following study by Wylie et al. Wylie et al. (2017) was the first to look into
state fatigue in people who have had a moderate-to-severe TBI. Furthermore,
they investigated the involvement of the Caudate nucleus in fatigue by seeing
if it changes its activity in direct proportion to the patients’ instantaneous
(state) fatigue experience Genova et al. (2013).

While the Caudate nucleus and the Striatum as a whole were previously thought
to be solely responsible for motor behavior control Middleton and Strick
(2000), recent evidence from animal and human research shows that this region
is involved in a wide range of cognitive behaviors, including learning Tricomi
and Fiez (2008); Dobryakova and Tricomi (2013); Schultz et al. (1993), outcome
processing Delgado et al. (2001, 2003), and working memory Lewis et al.
(2004); Akhlaghpour et al. (2016). Recent data indicate that fatigue caused by
such cognitive tasks may manifest itself in Caudate nucleus activation
Dobryakova et al. (2015). In children who have suffered a TBI, cognitive
fatigue has been linked to a network of areas in the Striatum and PFC,
including the vmPFC, nucleus accumbens, and Anterior Cingulate Cortex (ACC)
Ryan et al. (2016).

With the rapid increase in the availability of medical imaging datasets, deep
learning has been adopted efficiently to process the data for diagnosing
various diseases Ramesh Babu et al. (2020) and rehabilitation purposes
Farahanipad et al. (2020). Specifically, works have been done in predicting
diseases and subject traits using fMRI data with machine learning techniques
Pereira et al. (2009); Khosla et al. (2019). Further, convolutional neural
networks (CNNs) LeCun et al. (1995), an approach that has been known to be
very successful in solving computer vision tasks have been widely used to
analyze the spatial features in fMRI images as well. A 4-layer convolutional
neural network was proposed in Wang et al. (2018) for classification from raw
fMRI voxel values.

In Shen et al. (2019), the authors used deep convolutional networks (DCNs) to
encode fMRI images into low-dimensional feature space and decode them back for
image reconstruction. Similarly, the authors in Mozafari et al. (2020)
proposed a large-scale bi-directional generative adversarial network called
BigBiGAN to decode and reconstruct natural scenes from fMRI patterns.
Furthermore, an architecture based on sparse convolutional autoencoder was
used in Huang et al. (2017) to learn high-level features from handcrafted time
series derived from raw fMRI data.

There has also been a recent surge in the use of sequence models to process
temporal fMRI data. Mao et al. Mao et al. (2019) applied a specific type of
RNN known as Long Short-Term Memory (LSTM) to process spatial features
extracted from a CNN network. Another similar work in Thomas et al. (2018)
used bi-directional LSTM along with a CNN.

##  3 Methodology

To analyze cognitive fatigue from fMRI brain scans, data was recorded from
both TBI patients and healthy controls to study the differences in the extent
of cognitive fatigue between the two groups. The recorded data was stored in
NIfTI format (Neuroimaging Informatics Technology Initiative – a file format
for storing neuroimaging data). As raw data in NIfTI format contains noise, a
standard preprocessing pipeline was implemented to normalize and smoothen the
data as shown in figure 4. Both raw and pre-processed data were used to train
different deep learning models, and their results were compared. The following
sections explain data collection, pre-processing, and system architecture to
process the collected fMRI data.

Figure 1: A flow diagram of a series of N-back tasks performed during data
collection (SR Score: Self-Reported Fatigue Score)

###  3.1 Data Collection and Pre-processing

During data collection, fMRI scans of the brain were recorded over a period
where each subject was asked to perform a series of standardized cognitive
N-back tasks as shown in figure 1. Thirty participants with moderate-severe
TBI and 24 healthy controls (HCs) were involved in data collection. The
average age of the subjects was 41 years (SD=12.7). Each participant performed
four rounds of both 0-back and 2-back tasks. A resting-state fatigue score was
reported at the beginning, followed by scores being reported after each round.
Functional images were collected in 32 contiguous slices during eight blocks
(four at each of two difficulty levels), resulting in 140 acquisitions per
block (echo time = 30 ms; repetition time = 2000 ms; field of view = 22 cm;
flip angle = 80°; slice thickness = 4 mm, matrix = 64 × 64, in-plane
resolution = 3.438 $mm^{2}$). Using the Visual Analog Scale of Fatigue
(VAS-F), the subjects were asked to rate the amount of fatigue they
experienced (in the range 0-100) after each round of the N-back task. The
self-reported scores were mapped to six different classes to make it a multi-
class classification problem as represented in table 1.

Figure 2: Distribution of self-reported cognitive fatigue scores after every
N-back session from healthy subjects. The score is a difference between the
reported session score and the resting-state fatigue score recorded at the
beginning of the first session. Figure 3: Distribution of self-reported
cognitive fatigue scores after every N-back session from TBI subjects. The
score is a difference between the reported session score and the resting-state
fatigue score recorded at the beginning of the first session.

The visual inspection of the distribution (in fig 2 & 3) showed that six
categories struck a good balance/compromise between adequately describing the
distribution of VAS-F scores in a limited set of categories while also
maintaining sufficient complexity in the VAS-F data to allow for accurate
modeling. Reducing the number of categories to five or increasing it to seven
did not materially change the performance of the model. Additionally, the
cognitive fatigue levels shown in table 1 are for reference only so as to
quantify different levels of fatigue corresponding to the class label. The
size of the final 4D tensor acquired in NIfTI format was 140 x 32 x 64 x 64.
The raw fMRI images were preprocessed using Analysis of Functional NeuroImages
(AFNI) Cox (1996) and other standard techniques as discussed in previous works
Wylie et al. (2017), shown in figure 4.

###  3.2 System Architecture

Encoders are vital in any self-supervised pipeline as they map high-
dimensional input to a low-dimensional latent space ensuring that significant
features of the input are preserved. The choice of encoder architecture
determines how well a model learns the feature representations. In our
encoder, the output from a specific layer is pooled to get a single-
dimensional feature vector for every input sample.

Table 1: Mapping self-reported (SR) Cognitive Fatigue scores to respective class labels. The fatigue levels are for references only and are not of any clinical significance. Fatigue Score (SR) | Fatigue Level (Reference) | Class
---|---|---
0-10 | No Fatigue | 0
10-20 | Very Low Fatigue | 1
20-40 | Mild Fatigue | 2
40-60 | Fatigue | 3
60-80 | High Fatigue | 4
80-100 | Extreme Fatigue | 5

fMRI scans are 4D in shape and are represented as (t, x, y, z), where ’t’
represents the time-steps of individual 3D brain scans, and the other three
dimensions represent the intensity of voxels (1 x 1 x 1 mm) in the brain. The
temporal relation between the scans recorded at different time steps is
captured using a Recurrent Neural Network (RNN) based architecture. We
combined a CNN architecture with an LSTM network for the encoder as shown in
figure 5. We used three layers of 2D convolution and batch normalization to
learn the spatial features of the images.

The encoder was pre-trained on a public dataset called BOLD5000 Chang et al.
(2019) in a self-supervised approach and was fine-tuned on our labeled dataset
by adding a linear classifier layer at the end. BOLD5000 is a large-scale,
slow event-related human fMRI study incorporating 5,000 real-word images as
stimuli. It also accounts for image diversity, overlapping with standard
computer vision datasets, which makes it ideal for transfer learning tasks.
The image representations learned by the encoder by first pre-training on the
BOLD5000 dataset turned out to be more effective than being directly trained
on the supervised dataset.

Figure 4: Pre-processing pipeline for fMRI scans to convert from raw NIfTI
format to normalized and smoothed version Figure 5: Spatio-temporal Encoder
Architecture: CNN layers extract spatial features while LSTM layers model the
temporal aspect of the fMRI images followed by attention-based averaging over
time

####  3.2.1 Self-supervised Pre-training

The performance of any supervised learning method highly depends on the choice
and quality of data features and annotations available for the data.
Unfortunately, traditional supervised methods are hitting a bottleneck in
terms of performance because they rely on manually expensive annotated labels.
But, with the recent advancements in self-supervised approaches as discussed
in Jaiswal et al. (2021), unlabeled samples can be efficiently utilized to
learn data representations. They utilize meta-data generated from the dataset
that acts as pseudo-labels during the training process. One such approach is
contrastive learning which has proven to show great results in recent times to
learn image and video representations effectively. In this case, 4D fMRI data
can be treated as a series of videos such that self-supervised methods Zadeh
et al. (2020a) can be applied to it.

Figure 6: Self-supervised Pre-training Framework: MoCo algorithm for pre-
training on BOLD5000 dataset. The green arrows represent positive pairs and
red arrows represent negative pairs.

For pre-training the encoder, we used a contrastive-based approach where for
every sample in a batch containing $N$ samples, two augmented versions are
generated, resulting in 2N number of samples in the batch. For every sample,
its augmented version is considered as the positive candidate whose similarity
is encouraged to be maximum. In contrast, the similarity between the positive
and negative pair is encouraged to be as low as possible. This condition is
represented in figure 6 with green and red double-headed arrows. We use cosine
similarity to measure the closeness between two samples in the batch. Cosine
similarity between two vectors is the cosine of the angle between them and is
defined as,

$cos\\_sim(A,B)=\frac{A.B}{\|A\|\|B\|}$ (1)

During training, we applied extensive spatial and temporal augmentation for
the image data. As part of spatial transformation, methods such as random
affine, z-normalization, re-scale intensity, and one of random blur, random
gamma, random motion, and random noise were used. Similarly, a random starting
time t is selected for temporal augmentation, and ’n’ consecutive scans are
extracted. Finally, the loss is calculated using a variant of the Noise
Contrastive Estimation function (NCE) called InfoNCE, which is used when there
is more than one negative sample present during the learning process and is
defined by equation 2.

$L_{infoNCE}=-log\frac{exp({sim(q,{k_{+}})}/\tau)}{exp({sim(q,{k_{+}})}/\tau)+\sum_{i=0}^{K}exp({sim(q,{k_{i}})}/\tau)}$

(2)

In the equation 2, $q$ represents the current sample, $k_{+}$ represents the
positive sample (augmented version of $q$), and $k_{i}$ represents the
negative samples (other samples in the batch). $\tau$ represents the
temperature coefficient.

In this experiment, we used a variant of the contrastive approach called MoCo
He et al. (2020) that includes negative samples as a dictionary queue and has
been proven to be effective compared to other methods. Two encoders with the
same architectural configuration are used; the main encoder Q (Query Encoder)
is trained end-to-end on the sample pairs, while the second encoder (Momentum
Encoder) shares the same parameters as Q. The momentum encoder generates a
dictionary as a queue of encoded keys with the current mini-batch enqueued and
the oldest mini-batch dequeued. It gets updated based on the parameters of the
query encoder using an update parameter called momentum coefficient as
represented by equation 3. In equation 3, $m\in[0,1)$ is the momentum
coefficient. Only the parameters $\theta_{q}$ are updated by back-propagation.

$\theta_{k}\leftarrow m\theta_{k}+(1-m)\theta_{q}$ (3)

##  4 Experiments and Discussion

For fMRI images, most of the publicly available datasets contain data in NIfTI
format. We used two different data formats to train the models: one using the
raw NIfTI version and the other using pre-processed normalized version as
obtained from the preprocessing pipeline in fig 4. For self-supervised pre-
training, we only used NIfTI data from all four subjects in the BOLD5000
dataset. As shown in figure 6, the encoder was trained using MoCo He et al.
(2020) algorithm and Adam optimizer. The pre-training was carried out for a
total of 200 epochs. The starting learning rate was set to $0.03$ with a
weight decay factor of $10^{-4}$ and momentum parameter as 0.9. The learning
rate was decayed by a factor of 10 at 120 and 160 epochs, respectively.

Figure 7: Confusion matrix for self-supervised (fine-tuned) method for classes
0-5 (left-to-right) evaluated on the test set. The diagonal represents the
number of correctly classified instances

To train the deep learning models, we split our supervised labeled dataset
into three sets: train, validation, and test. The train set contained 70% of
the dataset, while the validation and the test datasets consisted of 15% each.
The test set contained a mix of both TBI and HC subjects and constituted more
than 300 reported instances during the N-back tasks as represented by the
confusion matrix in fig 7. On the other hand, scans from all four subjects in
the BOLD5000 dataset were used to pre-train the model using the self-
supervised approach, as mentioned earlier. For benchmarking, the primary
encoder was initially trained on our collected dataset separately using a
supervised approach. For the supervised method, we used both NIfTI and pre-
processed data to train two different models, as shown in table 2. Finally,
once the encoder was pre-trained on the BOLD5000 dataset using the self-
supervised algorithm, it was finetuned with the NIfTI version of our dataset.
We can analyze the performance of different models from table 2.

Table 2: Performance results for different models on cognitive fatigue classification task. Accuracies are calculated with 3-fold cross-validation Approach | Dataset Used | Data Format | Accuracy
---|---|---|---
HC only |  TBI only |  Overall
Supervised | Ours | NIfTI |  71.72 $\pm$ 0.82 |  78.44 $\pm$ 1.71 |  74.35 $\pm$ 1.27
Supervised | Ours | Pre-processed |  80.19 $\pm$ 0.63 |  84.17 $\pm$ 1.44 |  82.79 $\pm$ 0.73
Self-supervised + Fine-tuning | BOLD5000 + Ours | NIfTI |  85.47 $\pm$ 0.53 |  87.39 $\pm$ 1.26 |  86.84 $\pm$ 1.13

As shown by the results, the model pre-trained on BOLD5000 and later fine-
tuned on our supervised dataset outperformed the other supervised methods.
When testing the models on TBI and HC subjects separately, the models seem to
perform slightly better on the TBI data. This could mean that the enhanced
brain activations in the TBI subjects made it easier for the model to predict
the level of cognitive fatigue compared to the scans from healthy subjects.
However, the difference in the performance is negligent and the model seems to
perform comparatively well on data from both the subjects which makes it
robust for all cases. Also, based on the score distribution of TBI and HC
subjects in fig 2 & 3, TBI subjects seem to induce more fatigue compared to
healthy subjects.

Four NVIDIA GTX 1080 Ti GPUs were used to train the models, whereas, for
testing, only one GPU was used. As shown in table 2, our method beats all
previous approaches including Zadeh et al. (2020b) in classifying the level of
cognitive fatigue from fMRI scans.

##  5 Conclusion

This paper presented a spatio-temporal architecture pre-trained with self-
supervised methods for processing 4D fMRI data that predicts cognitive fatigue
in both TBI and healthy subjects. Unlike previous works that used masks to
focus on a particular brain region, our method learns essential areas of the
brain with maximum neural activations on its own. Motivated by video
classification, we used CNN layers to extract spatial features and an LSTM
network for temporal modeling of the data. Our architecture obtained state-of-
the-art performance for classifying different levels of cognitive fatigue from
fMRI data. Future works include exploring large-scale public datasets with
improved self-supervised algorithms to enhance the overall performance and the
granularity of cognitive fatigue prediction. Additionally, locating different
activation regions in the brain that induce cognitive fatigue can provide
better insights into how the brain works when handling cognitive tasks.

## References

  * Faul et al. [2010] Mark Faul, Marlena M Wald, Likang Xu, and Victor G Coronado.  Traumatic brain injury in the united states; emergency department visits, hospitalizations, and deaths, 2002-2006.  2010\.
  * Christodoulou et al. [2001] C Christodoulou, J DeLuca, JH Ricker, NK Madigan, BM Bly, G Lange, AJ Kalnin, WC Liu, J Steffener, BJ Diamond, et al.  Functional magnetic resonance imaging of working memory impairment after traumatic brain injury.  _Journal of Neurology, Neurosurgery & Psychiatry_, 71(2):161–168, 2001.
  * Kohl et al. [2009] Alexander D Kohl, Glenn R Wylie, HM Genova, Frank Gerard Hillary, and J Deluca.  The neural correlates of cognitive fatigue in traumatic brain injury using functional mri.  _Brain injury_ , 23(5):420–432, 2009.
  * DeLuca et al. [2008] John DeLuca, Helen M Genova, Frank G Hillary, and Glenn Wylie.  Neural correlates of cognitive fatigue in multiple sclerosis using functional mri.  _Journal of the neurological sciences_ , 270(1-2):28–39, 2008.
  * Wylie et al. [2017] GR Wylie, E Dobryakova, J DeLuca, N Chiaravalloti, K Essad, and H Genova.  Cognitive fatigue in individuals with traumatic brain injury is associated with caudate activation.  _Scientific reports_ , 7(1):1–12, 2017.
  * Chaudhuri and Behan [2000] Abhijit Chaudhuri and Peter Behan.  Fatigue and basal ganglia.  _Journal of the neurological sciences_ , 179:34–42, 11 2000\.  doi:10.1016/S0022-510X(00)00411-1.
  * Genova et al. [2013] Helen Genova, Venkateswaran Rajagopalan, John Deluca, Abhijit Das, Allison Binder, Aparna Arjunan, Nancy Chiaravalloti, and Glenn Wylie.  Examination of cognitive fatigue in multiple sclerosis using functional magnetic resonance imaging and diffusion tensor imaging.  _PloS one_ , 8:e78811, 11 2013.  doi:10.1371/journal.pone.0078811.
  * Middleton and Strick [2000] Frank Middleton and Peter Strick.  P.l. basal ganglia and cerebellar loops: motor and cognitive circuits.  _Brain research. Brain research reviews_ , 31:236–50, 04 2000.  doi:10.1016/S0165-0173(99)00040-5.
  * Tricomi and Fiez [2008] Elizabeth Tricomi and Julie Fiez.  Feedback signals in the caudate reflect goal achievement on a declarative memory task.  _NeuroImage_ , 41:1154–67, 08 2008.  doi:10.1016/j.neuroimage.2008.02.066.
  * Dobryakova and Tricomi [2013] Ekaterina Dobryakova and Elizabeth Tricomi.  Basal ganglia engagement during feedback processing after a substantial delay.  _Cognitive, affective and behavioral neuroscience_ , 13, 07 2013.  doi:10.3758/s13415-013-0182-6.
  * Schultz et al. [1993] W. Schultz, P. Apicella, and T. Ljungberg.  Responses of monkey dopamine neurons to reward and conditioned stimuli during successive steps of learning a delayed response task.  In _The Journal of neuroscience : the official journal of the Society for Neuroscience_ , 1993.
  * Delgado et al. [2001] M.R. Delgado, L.E. Nystrom, C. Fissell, Douglas Noll, and Julie Fiez.  Tracking the hemodynamic responses to reward and punishment in the striatum.  _Journal of neurophysiology_ , 84:3072–7, 01 2001.  doi:10.1152/jn.2000.84.6.3072.
  * Delgado et al. [2003] M Delgado, H Locke, VA Stenger, and Julie Fiez.  Dorsal striatum responses to reward and punishment: Effects of valence and magnitude manipulations.  _Cognitive, affective and behavioral neuroscience_ , 3:27–38, 04 2003.  doi:10.3758/CABN.3.1.27.
  * Lewis et al. [2004] S. Lewis, A. Dove, T. Robbins, R. Barker, and A. Owen.  Striatal contributions to working memory: a functional magnetic resonance imaging study in humans.  _European Journal of Neuroscience_ , 19, 2004.
  * Akhlaghpour et al. [2016] H. Akhlaghpour, J. Wiskerke, J. Y. Choi, Joshua P Taliaferro, Jennifer L. Au, and Ilana B. Witten.  Author response: Dissociated sequential activity and stimulus encoding in the dorsomedial striatum during spatial working memory.  _eLife_ , 2016.
  * Dobryakova et al. [2015] Ekaterina Dobryakova, Helen Genova, John Deluca, and Glenn Wylie.  The dopamine imbalance hypothesis of fatigue in multiple sclerosis and other neurological disorders.  _Frontiers in Neurology_ , 6, 03 2015.  doi:10.3389/fneur.2015.00052.
  * Ryan et al. [2016] Nicholas Ryan, Miriam Beauchamp, Richard Beare, Lee Coleman, Michael Ditchfield, Michael Kean, Timothy Silk, Sila Genc, Cathy Catroppa, and Vicki Anderson.  Uncovering cortico-striatal correlates of cognitive fatigue in pediatric acquired brain disorder: Evidence from traumatic brain injury.  _Cortex_ , 83, 08 2016.  doi:10.1016/j.cortex.2016.07.020.
  * Ramesh Babu et al. [2020] Ashwin Ramesh Babu, Mohammad Zaki Zadeh, Ashish Jaiswal, Alexis Lueckenhoff, Maria Kyrarini, and Fillia Makedon.  A multi-modal system to assess cognition in children from their physical movements.  In _Proceedings of the 2020 International Conference on Multimodal Interaction_ , pages 6–14, 2020.
  * Farahanipad et al. [2020] Farnaz Farahanipad, Harish Ram Nambiappan, Ashish Jaiswal, Maria Kyrarini, and Fillia Makedon.  Hand-reha: dynamic hand gesture recognition for game-based wrist rehabilitation.  In _Proceedings of the 13th ACM International Conference on PErvasive Technologies Related to Assistive Environments_ , pages 1–9, 2020.
  * Pereira et al. [2009] Francisco Pereira, Tom Mitchell, and Matthew Botvinick.  Machine learning classifiers and fmri: a tutorial overview.  _Neuroimage_ , 45(1):S199–S209, 2009.
  * Khosla et al. [2019] Meenakshi Khosla, Keith Jamison, Gia H Ngo, Amy Kuceyeski, and Mert R Sabuncu.  Machine learning in resting-state fmri analysis.  _Magnetic resonance imaging_ , 64:101–121, 2019.
  * LeCun et al. [1995] Yann LeCun, Yoshua Bengio, et al.  Convolutional networks for images, speech, and time series.  _The handbook of brain theory and neural networks_ , 3361(10):1995, 1995.
  * Wang et al. [2018] Xiaoxiao Wang, Xiao Liang, Yawen Zhou, Yanming Wang, Jin Cui, Huijuan Wang, Yu Li, Benedictor Alexander Nguchu, and Bensheng Qiu.  Task state decoding and mapping of individual four-dimensional fmri time series using deep neural network.  _arXiv preprint arXiv:1801.09858_ , 2018.
  * Shen et al. [2019] Guohua Shen, Tomoyasu Horikawa, Kei Majima, and Yukiyasu Kamitani.  Deep image reconstruction from human brain activity.  _PLoS computational biology_ , 15(1):e1006633, 2019.
  * Mozafari et al. [2020] Milad Mozafari, Leila Reddy, and Rufin VanRullen.  Reconstructing natural scenes from fmri patterns using bigbigan.  In _2020 International Joint Conference on Neural Networks (IJCNN)_ , pages 1–8. IEEE, 2020.
  * Huang et al. [2017] Heng Huang, Xintao Hu, Yu Zhao, Milad Makkie, Qinglin Dong, Shijie Zhao, Lei Guo, and Tianming Liu.  Modeling task fmri data via deep convolutional autoencoder.  _IEEE transactions on medical imaging_ , 37(7):1551–1561, 2017.
  * Mao et al. [2019] Zhenyu Mao, Yi Su, Guangquan Xu, Xueping Wang, Yu Huang, Weihua Yue, Li Sun, and Naixue Xiong.  Spatio-temporal deep learning method for adhd fmri classification.  _Information Sciences_ , 499:1–11, 2019.
  * Thomas et al. [2018] Armin W Thomas, Hauke R Heekeren, Klaus-Robert Müller, and Wojciech Samek.  Interpretable lstms for whole-brain neuroimaging analyses.  _arXiv preprint arXiv:1810.09945_ , 2018.
  * Cox [1996] Robert W Cox.  Afni: software for analysis and visualization of functional magnetic resonance neuroimages.  _Computers and Biomedical research_ , 29(3):162–173, 1996.
  * Chang et al. [2019] Nadine Chang, John A Pyles, Austin Marcus, Abhinav Gupta, Michael J Tarr, and Elissa M Aminoff.  Bold5000, a public fmri dataset while viewing 5000 visual images.  _Scientific data_ , 6(1):1–18, 2019.
  * Jaiswal et al. [2021] Ashish Jaiswal, Ashwin Ramesh Babu, Mohammad Zaki Zadeh, Debapriya Banerjee, and Fillia Makedon.  A survey on contrastive self-supervised learning.  _Technologies_ , 9(1):2, 2021.
  * Zadeh et al. [2020a] Mohammad Zaki Zadeh, Ashwin Ramesh Babu, Ashish Jaiswal, and Fillia Makedon.  Self-supervised human activity recognition by augmenting generative adversarial networks.  _arXiv preprint arXiv:2008.11755_ , 2020a.
  * He et al. [2020] Kaiming He, Haoqi Fan, Yuxin Wu, Saining Xie, and Ross Girshick.  Momentum contrast for unsupervised visual representation learning.  In _Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition_ , pages 9729–9738, 2020.
  * Zadeh et al. [2020b] Mohammad Zaki Zadeh, Ashwin Ramesh Babu, Jason Bernard Lim, Maria Kyrarini, Glenn Wylie, and Fillia Makedon.  Towards cognitive fatigue detection from functional magnetic resonance imaging data.  In _Proceedings of the 13th ACM International Conference on PErvasive Technologies Related to Assistive Environments_ , pages 1–2, 2020b.


# Quarkyonic stars with isospin-flavor asymmetry

Jérôme Margueron  Univ Lyon, Univ Claude Bernard Lyon 1, CNRS, IP2I Lyon /
IN2P3, UMR 5822, F-69622, Villeurbanne, France  Hubert Hansen  Univ Lyon, Univ
Claude Bernard Lyon 1, CNRS, IP2I Lyon / IN2P3, UMR 5822, F-69622,
Villeurbanne, France  Paul Proust  Univ Lyon, Univ Claude Bernard Lyon 1,
CNRS, IP2I Lyon / IN2P3, UMR 5822, F-69622, Villeurbanne, France  Guy Chanfray
Univ Lyon, Univ Claude Bernard Lyon 1, CNRS, IP2I Lyon / IN2P3, UMR 5822,
F-69622, Villeurbanne, France

###### Abstract

We suggest an extension to isospin asymmetric matter of the quarkyonic model
from McLerran and Reddy McLerran:2019 . This extension allows us to construct
the $\beta$-equilibrium between quarks, nucleons and leptons. The concept of
the quarkyonic matter originates from the large number of color limit for
which nucleons are the correct degrees of freedom near the Fermi surface –
reflecting the confining forces – while deep inside the Fermi sea quarks
naturally appear. In isospin asymmetric matter, we suggest that this new
concept can be implemented within a global isoscalar relation between the
shell gaps differentiating the nucleon and the quark sectors. In addition, we
impose the conservation of the isospin-flavor asymmetry in the nucleon and the
quark phases. Within this model, several quarkyonic stars are constructed on
top of the SLy4 model for the nucleon sector, producing a bump in the sound
speed, which implies that quarkyonic stars are systematically bigger and have
a larger maximum mass than the associated neutron stars. They also predict
lower proton fraction at $\beta$-equilibrium, which potentially quenches fast
cooling in massive compact stars.

Recent observations of neutron stars (NS), such as radio and x-ray astronomy,
or the detection of gravitational wave (GW) detection, have provided the
tightest constraints on the dense matter equation of state (EoS) to date
NewCompStar2019 ; HPY:2007 ; Ligo:2017 . These constraints can be classified
into different groups: the first one refers to the highest NS masses ever
observed Antoniadis:2013 ; Ozel:2016 ; Fonseca:2016 ; Linares:2018 ;
Cromartie:2018 ; Arzoumanian:2018 , estimated to be about $2M_{\odot}$, with
some indications that the maximum mass could eventually be larger
Cromartie:2018 ; Arzoumanian:2018 . The second group assembles constraints
from binary NS (BNS) GW detection in the in-spiral phase, from which the tidal
deformability is estimated Hinderer:2008 ; Flanagan:2008  . It includes
GW170817 Ligo:2017  and further detections. The third group still refers to NS
mergers, more specifically to the analysis of the electromagnetic (EM)
counterpart, see Ref. Coughlin:2019  for instance and refs. therein. The
fourth group matters with x-ray observations, such as thermal emission from
qLMXB Guillot:2013 ; Steiner:2018 ; Baillot:2019 , x-ray burst and
photospheric expansion Ozel:2016 ; Carolyn:2016 , as well as x-ray emission
from hot-spots at the surface of some NS (NICER) NICER .

The analysis of these new data requires to set-up generic models for the EoS,
with various levels of agnosticity to model assumptions, such as the isospin
asymmetry for instance, or the interaction prerequisites, inducing possible
spurious constraints among observables. As an illustration, the nuclear meta-
model Margueron:2018a  can be used to explore predictions compatible with the
assumption that matter is uniquely composed of nucleons (and of course
leptons). It can also be employed for the description of phase transitions.
This is the purpose of this paper: exploring the nucleon-quark phase
transition in compact stars where this transition is described with the
quarkyonic model McLerran:2019 .

The quarkyonic model for dense matter proposed in Ref. McLerran:2007  (and
recently applied to neutron stars in Ref. McLerran:2019 ) is one of them. It
is an interesting candidate to bridge the gap in describing quark matter and
nuclear matter at the phase transition Fukushima:2016 . The quarkyonic model
is not properly a microscopic model since it is not based on the QCD
Lagrangian or an effective version of it, but it implements some features from
in the large number of color ($N_{c}$) limit of QCD McLerran:2007  with few
parameters. New configurations at high and low temperature limit of the
holographic Witten-Sakai-Sugimoto model has recently been interpreted as
holographic realizations of quarkyonic matter in isospin symmetric matter,
based on a quark Fermi sea enclosed by a baryonic layer on momentum space
Kovensky2020 . In the real-world where $N_{c}=3$, it possibly approximates the
actual ground state of dense matter and the confrontation of its predictions
to the data can be used to determine the model parameters. These parameters
happen to be physical and thus meaningful: they are the quarkyonic scale
$\Lambda_{\mathrm{Qyc}}\approx 250-300$ MeV, which is comparable to the QCD
scale $\Lambda_{\mathrm{QCD}}$, and $\kappa_{\mathrm{Qyc}}$ which controls the
saturation of the nucleonic shell McLerran:2007 .

The interesting feature of the quarkyonic model is that it suggests a cross-
over between the hadron phase and the quark one, at variance with other
approaches such as the ones based on Maxwell/Gibbs construction. Note that
cross-over are also suggested from Cooper pairing in the hadron and in the
quark phases. The quarkyonic cross-over is thus another example of such
features. It however disregards one of the essential prediction of QCD, namely
the restoration – as the density increases – of chiral symmetry which is
spontaneously broken in the QCD vacuum. Note that in the holographic approach
from Ref. Kovensky2020 , chirally restored and chirally broken quarkyonic
matter are constructed, and it was found that only chirally restored matter is
energetically preferred. In widely used quark models implementing the chiral
symmetry breaking, e.g. Nambu–Jona-Lasinio approach, the transition between
the broken phase assimilated to the hadronic phase and the restored one
assimilated to the free quark phase is however generally first order (for most
parameter sets). In the future, it would be interesting to combine together in
isospin asymmetric matter the phenomenology of the color gauge symmetry
realized at large $N_{c}$ and the chiral symmetry dynamics, both rooting into
the QCD Lagrangian.

In the cross-over region, the quarkyonic model suggests that the pressure
first increases at the onset of the first quarks, while first order phase
transition models usually suggest a softening of the pressure due to the
increase of the degrees of freedom Zdunik:2013 ; Alford:2013 . The consequence
is the large increase of the energy density, as well as a peaked sound speed
at a density of about 2-3 times the saturation density of nuclear matter,
$n_{\mathrm{sat}}\approx 0.16\,\text{fm}^{-3}$, as expected by some authors
Steiner:2015 ; Tews:2018 .

These features are characteristic to the quarkyonic model – however not
specific – and have motivated further investigations and applications to
neutron star physics. In the original paper by McLerran and Reddy
McLerran:2019 , the quarkyonic matter was studied in the case of isospin
symmetry (symmetric matter, SM) as well as in the case of neutron matter (NM).
The application to neutron star has been performed assuming it is composed
only of neutrons, and $u$ and $d$ quarks within the ratio satisfying local
charge neutrality, $k_{Fd}=2^{1/3}k_{Fu}$ where $k_{Fd}$ ($k_{Fu}$) is the $d$
($u$) quark Fermi momentum. A version of the quarkyonic model for isospin
asymmetric matter, where the isospin asymmetry is controlled by the chemical
equilibrium, was then suggested by Zhao and Lattimer Zhao:2020 . In their
model, Zhao and Lattimer have treated nucleons and quarks as independent
particles for which the energy minimization imposes the equilibrium between
their respective chemical potential. They found that their quarkyonic stellar
model is able to satisfy observed mass and radius constraints with a wide
range of model parameters. They also predict that quarkyonic matter tends to
reduce the proton fraction, compared to the nucleonic case. This reduces the
domain of parameter allowing the direct URCA process Lattimer:1991 . It was
also suggested by Jeong, McLerran and Sen Jeong:2020  that the hard core in
the nucleon interaction could be represented by an excluded volume, which in
turns can be related to the shell gap controlling the cross-over properties.
In such model, the shell gap is directly controlled by the size of the hard
core. Similarly to the original paper by McLerran and Reddy McLerran:2019 ,
this model also predicts the presence of a peak in the sound speed at
2-3$n_{\mathrm{sat}}$. It was also extended to describe three-flavor baryon-
quark mixtures, allowing the onset of strange particles Duarte:2020a ;
Duarte:2020b .

In this paper, we suggest another version of the quarkyonic model for isospin
asymmetric matter (AM) where we investigate the analogy between the quarkyonic
model and the Cooper pair formation around the Fermi energy McLerran:2019 .
While this analogy may appear as rather simplistic, it suggests that quarks
and nucleons may be viewed as two representations of the same quasi-particle
excitation, in the same way as Cooper pairs and single particles coexist in
superfluids or superconductors, but in different regions of the nuclear
spectrum deGennes . In AM, the neutron/proton ratio in the nucleon sector and
the flavor asymmetry in the quark sector are fixed by the compound nature of
the nucleons, since $n:udd$ and $p:uud$. In this spirit, quarks and nucleons
are not distinguished as two independent particles for which the energy
minimization imposes an equilibrium relation, as suggested in the quarkyonic
model of Zhao and Lattimer Zhao:2020 . In addition, the $\beta$-equilibrium
does not involve quark chemical potentials since only nucleons are occupying
the Fermi levels. Our picture requires a new approach for the thermodynamical
construction of the phase equilibrium. The cross-over, as described in the
original quarkyonic model McLerran:2019 , is depicted by an isoscalar
condition connecting the momenta of the quarks and of the nucleons, while the
isospin/flavor asymmetry in the quark and the nucleon sectors is fixed. Under
these two assumptions, the model we propose describes the phase transition
from symmetric to neutron matter.

In our picture, there is no direct contribution of the quarks to the
$\beta$-equilibrium since they do not occupy Fermi levels. The presence of
quarks however influences the $\beta$-equilibrium, through their contribution
to the nucleon chemical potentials. This picture breaks down in the pure quark
phase, which does not occurs in the quarkyonic model since there is always a
small but finite contribution of nucleons at high density. In addition the
chiral symmetry generating the constituent quark mass is assumed to remain at
all density, even in the dense phases where quarks become the dominant
species. This is also an interesting suggestion from the quarkyonic model
which goes against the usual picture of the hadron-quark phase transition
based on chiral symmetry restoration. As we stated earlier, on the one hand
the quarkyonic transition is driven by features of QCD relying on its gauge
theory nature at large $N_{c}$, where only planar graphs survive, whereas on
the other hand the transition follows the chiral symmetry restoration
(property of the quark sector only) which induces a large change of the
constituent masses and also of the baryon properties. In the future, a model
unifying both mechanisms in isospin asymmetric matter would be an interesting
theoretical development.

In this paper, we suggest an extension of the original quarkyonic model
McLerran:2019  for AM in Sec. I. The cold catalyzed NS EoS – at
$\beta$-equilibrium – is derived in Sec. II and finally we calculate NS
properties in Sec. III. We finally conclude and suggest outlook in Sec. IV.

##  I Quarkyonic model in asymmetric matter

The concept of quarkyonic matter has emerged in the large number of color,
$N_{c}$, limit of QCD McLerran:2007 . In this limit and when the nucleon
density $n_{N}$ is much larger than the QCD scale,
$n_{N}\gg\Lambda_{\mathrm{QCD}}^{3}$, the confining potential of QCD is
dominant even though the nucleonic Fermi momentum is large,
$k_{F_{N}}\gg\Lambda_{\mathrm{QCD}}$. The concept of quarkyonic matter has
been introduced in order to resolve this apparent paradox: the ground state of
dense matter is composed of dressed quarks (with mass $\simeq M_{N}/3$) that
are freely moving deep inside the Fermi sea, and of a shell of baryons
generated by the strong confining force, which lies close to the Fermi level
McLerran:2007 . Baryons occupies a momentum shell whose width is taken to be
$\Delta_{\mathrm{Qyc}}\approx\Lambda_{\mathrm{Qyc}}$, while deepest states are
occupied by quarks.

Even if matter is composed of nucleons and quarks, low energy excitations
around the Fermi level involve only quasi-particles of nucleonic type. The
excitation of quarks require the expense of a momentum of the order of the
shell width $\Delta_{\mathrm{Qyc}}$. In this prospect, nucleons and quarks are
not independent particles, but there are realization of matter excitations in
well separated energy regions. This picture is clearly illustrated on the left
panel of Fig. 1 for SM.

In SM, this picture is simple since it is sufficient to fix a condition
between the nucleon and quark single particle energies $\epsilon_{N}(k)$ and
$\epsilon_{Q}(k)$ to separate the deep quark states from the nucleons ones
located around the Fermi level. This condition imposes that the last quark
occupied state coincides with the first nucleon one McLerran:2019 ,

$\epsilon_{N}^{\mathrm{NI}}(k_{F_{N}}-\Delta_{\mathrm{Qyc}})=N_{c}\epsilon_{Q}(k_{F_{Q}})$
(1)

where $k_{F_{N}}$ and $k_{F_{Q}}$ are isoscalar nucleon and quark Fermi
momenta. The factor $N_{c}$ in Eq. (1) stands for the number of quarks forming
each nucleon in its ground state. With this condition, the nucleon and quark
degrees of freedom are reduced to only one free variable, that we fix to be
$k_{F_{N}}$, the nucleon Fermi momentum for convenience. Note that in Eq. (1)
the nucleon single particle energy is the non-interacting (NI) one,
$\mathcal{E}_{N}^{\mathrm{NI}}(k)=\sqrt{M_{N}^{2}+k^{2}}$. This choice is
performed in order to avoid unnecessary complication of the model. We further
assume that chiral symmetry remains broken to set the quark mass
($M_{u}=M_{d}=M_{Q}$) with $M_{Q}\approx M_{N}/N_{c}$, as in the constituent
quark model.

Finally a prescription for the thickness of Fermi layer where nucleons reside
has to be taken. We adopt the same relation from the original paper
McLerran:2019 ,

$\displaystyle\Delta_{\mathrm{Qyc}}=\frac{\Lambda_{\mathrm{Qyc}}^{3}}{\hbar
c^{3}k_{F_{N}}^{2}}+\kappa_{\mathrm{Qyc}}\frac{\Lambda_{\mathrm{Qyc}}}{\hbar
cN_{c}^{2}}\,,$ (2)

with $\Lambda_{\mathrm{Qyc}}\approx 250-300$ MeV and
$\kappa_{\mathrm{Qyc}}\approx 0.3$. $\Delta_{\mathrm{Qyc}}$ defines the energy
scale below which nucleons reside.

The concept of quarkyonic matter however leads to a fundamentally new way to
represent the nucleon and quark densities $n_{N}$ and $n_{Q}$ and their
associated Fermi seas. As in a superfluid, there is only one chemical
potential that enters into the thermodynamic equilibrium, and which is
associated to the last nucleon occupied state $\mu_{N}=E(N_{B})-E(N_{B}-1)$
where $E$ is the total baryon energy in the ground state and $N_{B}$ stands
for the number of baryons which is a partition of nucleons and quarks states.
By adding or removing a baryon, the entire Fermi sea is reorganized leading to
a new partition between quark and nucleon states, with the condition on the
number densities $n_{B}=n_{N}+n_{Q}$. This new concept is first applied to
symmetric and neutron matter McLerran:2019  and we are now suggesting an
extension for AM.

Figure 1: Schematic views of the different Fermi seas at stake in quarkyonic
matter for SM (left) and AM (right). Quarks occupy deep states inside their
Fermi spheres while nucleons occupy the shells close to their Fermi levels.

###  I.1 Global isoscalar relation between nucleon and quark Fermi seas

By breaking the isospin symmetry in AM, the nucleon and quark states are
replaced by four other states: neutrons, protons, as well as $u$ and $d$
quarks are the natural components in AM, where they are represented by their
associated densities, $n_{n}$, $n_{p}$, $n_{u}$ and $n_{d}$. The four Fermi
seas are schematically represented in the right panel of Fig. 1. The baryon
charge density controlling SM is completed in AM by the isospin asymmetry
$\delta_{N}$. Since there are only two charges and four particles, the concept
of quarkyonic matter in asymmetric matter requires two additional relations.

At variance with SM where nucleons and quark Fermi seas could be defined by
imposing a relation between their single particle energies, see Eq. (1),
nothing similar can be done in AM between the four particles. The concept of
quarkyonic matter suggests however that the relation between nucleons and
quarks in SM, see Eq. (1), remains globally valid in AM. The nucleon and quark
Fermi momenta in AM, $k_{F_{N}}$ and $k_{F_{Q}}$, are represented in Fig. 1 by
the red dashed circles.

Expanding the single particle energies in Eq. (1), e.g.
$\mathcal{E}_{Q}(k)=\sqrt{M_{Q}^{2}+k^{2}}$, one gets the following relation
between the nucleon and quark isoscalar Fermi momenta,

$\displaystyle k_{F_{Q}}$ $\displaystyle=$
$\displaystyle\frac{k_{F_{N}}-\Delta_{\mathrm{Qyc}}}{N_{c}}\Theta(k_{F_{N}}-\Delta_{\mathrm{Qyc}})\,.$
(3)

Note that in Eq. (3) we assumed that the nucleon shell gap
$\Delta_{\mathrm{Qyc}}$ is also an isoscalar quantity. Eq. (3) allows quark to
appear as soon as $k_{F_{N}}-\Delta_{\mathrm{Qyc}}>0$ independently of the
isospin asymmetry. This is supported, as we explain previously, by the idea
that exciting quarks in AM requires an energy of the order of
$\Delta_{\mathrm{Qyc}}$ irrespectively of the partition of matter between
neutrons and protons. Without an actual solution of QCD, this assumption is
the simplest one which can be done.

The limit of NM was also explored in the original paper by McLerran and Reddy
McLerran:2019 . The prescription taken there imposes that one has to choose
which among $u$ and $d$ quarks are connected to the neutron states. In
addition, the value for $\Lambda_{\mathrm{Qyc}}=300$ MeV considered in SM was
changed to $380$ MeV in NM. In our model, $\Lambda_{\mathrm{Qyc}}$ is taken
constant, but the change of the Fermi momentum between SM and NM, see Eqs. (3)
and (17), induces an effective modification of $\Lambda_{\mathrm{Qyc}}$
between SM and NM with the ratio $2^{1/3}$, which is exactly the same ratio
considered by McLerran and Reddy McLerran:2019 . In practice, the two
approaches lead to similar results in NM. The isoscalar relation (3) presents
however the advantage to describe AM and to recover the limit of SM, where the
concept of the quarkyonic matter is simple.

In summary, we remark that the isoscalar Fermi momentum $k_{F_{N}}$ controls
both the isoscalar quark Fermi momentum $k_{F_{Q}}$, see Eq. (3), as well as
the nucleon gap $\eqref{eq:kfq}$ from the prescription (2). The isoscalar
nucleon and quark densities can be safely determined as,

$\displaystyle
n_{N}=\frac{2}{3\pi^{2}}\Big{[}k_{F_{N}}^{3}-(k_{F_{N}}-\Delta_{\mathrm{Qyc}})^{3}\Theta(k_{F_{N}}-\Delta_{\mathrm{Qyc}})\Big{]}\,,$
(4)

and the quark density,

$\displaystyle n_{Q}$ $\displaystyle=$
$\displaystyle\frac{2}{3\pi^{2}}k_{F_{Q}}^{3}\Theta(k_{F_{Q}})\,.$ (5)

The total baryon density is built upon the nucleon and quark contributions as,

$\displaystyle n_{B}=n_{N}+n_{Q}\,.$ (6)

While the isoscalar densities $n_{N}$ and $n_{Q}$ can now be calculated from
$k_{F_{N}}$, the connection to the densities of the four particles $n$, $p$,
$u$ and $d$ are yet unknown. They are related to the isoscalar densities from
the following relation:

$\displaystyle n_{N}$ $\displaystyle=$ $\displaystyle n_{n}+n_{p}\,$ (7)
$\displaystyle n_{Q}$ $\displaystyle=$ $\displaystyle(n_{d}+n_{u})/N_{c}\,.$
(8)

In the following, we suggest that the densities $n_{n}$, $n_{p}$, $n_{u}$ and
$n_{d}$ can be obtained in AM by imposing the isospin/flavor asymmetry in the
nucleon and quark phases.

###  I.2 Isospin/flavor asymmetry

We will now determine the particle densities. In the quarkyonic model, there
is a partition between nucleons and quarks which changes as the density
evolves. It reveals the dynamical process converting quarks into nucleons or
breaking the nucleons into their constituents, as it shall be in compound
systems. A simple way to translate this symmetry into the nucleon and quark
phases is to impose that the two phases conserve the isospin/flavor asymmetry.
Since $n:(udd)$ and $p:(uud)$, we obtain the following relations for the quark
density in the nucleon phase, $n_{u}^{\mathrm{nuc}}=n_{n}+2n_{p}$ and
$n_{d}^{\mathrm{nuc}}=2n_{n}+n_{p}$, which leads to the following simple
connection between the isospin asymmetry parameter
$\delta_{N}=(n_{n}-n_{p})/n_{N}$ in the nucleon phase and the flavor asymmetry
parameter $\delta_{Q}=(n_{d}-n_{u})/(n_{d}+n_{u})$ in the quark phase:

$\delta_{N}=N_{c}\delta_{Q}\,.$ (9)

From SM to NM, $\delta_{N}$ goes from $0$ to $1$, while $\delta_{Q}$ goes from
$0$ to $1/N_{c}$. The dynamics of the phase transition thus imposes that the
$u$ and $d$ flavor ratio reflects the isospin asymmetry of the nucleon phase.

As a side note, we are aware that the isospin/flavor asymmetry relation (9)
can possibly be violated by the two phases if the energy minimization is
injected. Such a refinement of the quarkyonic model is indeed very interesting
but it complexifies the quarkyonic approach, which nice feature remains in its
simplicity. Further extensions of the present model will be explored in the
future, especially to analyze their role in the predictions presented here.

Knowing $\delta_{N}$ and $k_{F_{N}}$ – which fixes $n_{N}$ and $n_{Q}$ – one
can deduce all particle densities as

$\displaystyle n_{n}$ $\displaystyle=$
$\displaystyle\frac{1+\delta_{N}}{2}n_{N}\equiv x_{n}n_{N}\,,$ (10)
$\displaystyle n_{p}$ $\displaystyle=$
$\displaystyle\frac{1-\delta_{N}}{2}n_{N}\equiv x_{p}n_{N}\,,$ (11)
$\displaystyle n_{d}$ $\displaystyle=$
$\displaystyle\frac{1+\delta_{Q}}{2}N_{c}n_{Q}\equiv x_{d}N_{c}n_{Q}\,,$ (12)
$\displaystyle n_{u}$ $\displaystyle=$
$\displaystyle\frac{1-\delta_{Q}}{2}N_{c}n_{Q}\equiv x_{u}N_{c}n_{Q}\,.$ (13)

The $u$ and $d$ quark Fermi momenta are simply related to their densities as

$\displaystyle k_{F_{u}}^{3}$ $\displaystyle=$
$\displaystyle\frac{3\pi^{2}}{N_{c}}n_{u}=(1-\delta_{Q})k_{F_{Q}}^{3}\,,$ (14)
$\displaystyle k_{F_{d}}^{3}$ $\displaystyle=$
$\displaystyle\frac{3\pi^{2}}{N_{c}}n_{d}=(1+\delta_{Q})k_{F_{Q}}^{3}\,,$ (15)

since $d$ and $u$ quarks occupy their Fermi sphere, see Fig. 1.

The neutron and proton Fermi layers can be calculated from the difference of
two Fermi spheres with different radii defined as,

$n_{n}=\frac{1}{3\pi^{2}}\left(k_{F_{n}}^{3}-k_{F_{n}^{\min}}^{3}\right)\,,\,\,n_{p}=\frac{1}{3\pi^{2}}\left(k_{F_{p}}^{3}-k_{F_{p}^{\min}}^{3}\right)\,,$
(16)

where $k_{F_{n}^{\min}}$ and $k_{F_{p}^{\min}}$ are the lower bound of the
nucleon shell, see Fig. 1. Injecting Eq. (4) into $n_{i}=x_{i}n_{N}$ ($i=n$,
$p$) and identifying with Eqs. (16), we obtain

$k_{F_{n}}^{3}=(1+\delta_{N})k_{F_{N}}^{3}\,,\hskip
14.22636ptk_{F_{p}}^{3}=(1-\delta_{N})k_{F_{N}}^{3}\,,$ (17)

as well as

$\displaystyle k_{F_{n}^{\min}}^{3}$ $\displaystyle=$
$\displaystyle(1+\delta_{N})\left(N_{c}k_{F_{Q}}\right)^{3}\,,$ (18)
$\displaystyle k_{F_{p}^{\min}}^{3}$ $\displaystyle=$
$\displaystyle(1-\delta_{N})\left(N_{c}k_{F_{Q}}\right)^{3}\,.$ (19)

Note that $k_{F_{n}^{\min}}$ and $k_{F_{p}^{\min}}$ can be re-expressed as
$k_{F_{i}}^{\min}=(2x_{i})^{1/3}N_{c}k_{F_{Q}}=N_{c}(3\pi^{2}x_{i}n_{Q})^{1/3}$,
for $i=n$, $p$.

At low densities, in the absence of quarks $k_{F_{Q}}=0$, so neutrons and
protons occupy entirely their Fermi spheres with radii given by Eqs. (17). It
is interesting to note that Eqs. (17) are identical in the presence or the
absence of quarks.

Now that the Fermi spheres and the Fermi shells for $n$, $p$, $u$ and $d$
particles are well defined, all thermodynamical quantities can be determined,
e.g. the energy of the ground state, the pressure, the chemical potentials, or
the sound speed.

###  I.3 Energy density and derivatives

Model | $E_{\mathrm{sat}}$ | $E_{\mathrm{sym}}$ | $n_{\mathrm{sat}}$ | $L_{\mathrm{sym}}$ | $K_{\mathrm{sat}}$ | $K_{\mathrm{sym}}$ | $Q_{\mathrm{sat}}$ | $Q_{\mathrm{sym}}$ | $Z_{\mathrm{sat}}$ | $Z_{\mathrm{sym}}$ | $m^{*}/m$ | $\Delta m^{*}/m$ | $b_{\mathrm{sat}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---
| MeV | MeV | $\,\text{fm}^{-3}$ | MeV | MeV | MeV | MeV | MeV | MeV | MeV |  |  |
SLy4${}_{\mathrm{MM}}^{\mathrm{NR}}$ | -15.97 | 32.01 | 0.1595 | 46 | 230 | -120 | -225 | 400 | -443 | -690 | 0.69 | -0.19 | 6.90
SLy4${}_{\mathrm{MM}}^{\mathrm{RL}}$ | -15.97 | 32.01 | 0.1595 | 46 | 230 | -120 | -225 | 400 | -443 | -690 | 1.0 | 0.0 | 6.90
Table 1: Parameters of the SLy4 meta-model used in the non-relativistic (NR)
case for the description of nuclear matter and in the relativistic case (RL),
where only relativistic kinematic is considered, for the quarkyonic matter.

The energy density of quarkyonic matter is given by

$\displaystyle\rho_{B}$ $\displaystyle=$ $\displaystyle\rho_{N}+\rho_{Q}\,,$
(20)

where the nucleon and quark terms are given by:

$\displaystyle\rho_{N}$ $\displaystyle=$ $\displaystyle
2\\!\\!\sum_{i=n,p}\\!\int_{k_{F_{i}}^{\min}}^{k_{F_{i}}}\\!\\!\frac{d^{3}k}{(2\pi)^{3}}\sqrt{k^{2}+M_{N}^{2}}+V_{N}(n_{n},n_{p})\,,$
(21) $\displaystyle\rho_{Q}$ $\displaystyle=$ $\displaystyle
2\\!\\!\sum_{q=u,d}\\!\\!N_{c}\\!\\!\int_{0}^{k_{F_{q}}}\\!\\!\\!\frac{d^{3}k}{(2\pi)^{3}}\sqrt{k^{2}+M_{Q}^{2}}\,,$
(22)

with the nuclear residual interaction given by the meta-model Margueron:2018a
,

$\displaystyle
V_{N}(n_{n},n_{p})=\sum_{\alpha}\frac{1}{\alpha!}c_{\alpha}(\delta_{N})x^{\alpha}u_{\alpha}(x)\,,$
(23)

where $x=(n_{N}-n_{\mathrm{sat}})/(3n_{\mathrm{sat}})$,
$c_{\alpha}(\delta_{N})=c_{\mathrm{sat},\alpha}+c_{\mathrm{sym},\alpha}\delta_{N}^{2}$
and $u_{\alpha}(x)=1-(-3x)^{N+1-\alpha}\exp(-b\,n_{N}/n_{\mathrm{sat}})$. The
coefficients $c_{\mathrm{sat},\alpha}$ and $c_{\mathrm{sym},\alpha}$ are
related to the empirical parameters, e.g. $E_{\mathrm{sat}}\approx-16$ MeV,
$n_{\mathrm{sat}}\approx 0.16$ fm-3 and $K_{\mathrm{sat}}\approx 230$ MeV in
nuclear matter, considering the relativistic extension of the meta-model (for
the kinetic term only) RMM . The present calculation is based on the nucleon
Skyrme interaction SLy4, on which the meta-model is adjusted. We consider the
parameters given in table 1.

The binding energy density $\epsilon_{B}$ is defined as

$\displaystyle\epsilon_{B}(k_{F_{N}},\delta_{N})=\rho_{B}(k_{F_{N}},\delta_{N})-M_{N}n_{B}\,,$
(24)

and the binding energy per baryon number is

$\displaystyle
e_{B}(k_{F_{N}},\delta_{N})=\epsilon_{B}(k_{F_{N}},\delta_{N})/n_{B}\,.$ (25)

The other quantities as the chemical potentials, pressure and sound velocities
are computed using the usual definitions, see Ref. Margueron:2018a ;
Margueron:2018b  for more details.

###  I.4 Results

Figure 2: Baryon chemical potential $\mu_{B}$, energy per particle $E_{B}/A$,
pressure $P_{B}$ and sound speed $(v_{s,B}/c)^{2}$ in symmetric (SM,
$\delta_{N}=0.0$), asymmetric (AM, $\delta_{N}=0.5$), and neutron matter (NM,
$\delta_{N}=1.0$) for $\Lambda_{\mathrm{Qyc}}=250$ MeV and
$\kappa_{\mathrm{Qyc}}=0.3$.

In this section we compare a pure nucleon model for matter properties against
quarkyonic models constructed on top of the same nucleon model. The choice of
SLy4 is influenced by the fact that this pure nucleon model reproduces most of
the recent observational data, such as the tidal deformability from GW170817
or the NICER x-ray observation of PSR J0030+0451. We then explore the
influence of the quarkyonic model parameters $\Lambda_{\mathrm{Qyc}}$ and
$\kappa_{\mathrm{Qyc}}$.

The neutron chemical potential $\mu_{n}$, the energy per particle $E_{B}/A$,
the baryon pressure $P_{B}$ and baryon sound speed $(v_{s,B}/c)^{2}$ are shown
in Fig. 2 for SM (solid lines, $\delta_{N}=0$), AM (dashed lines,
$\delta_{N}=0.5$) and NM (dotted lines, $\delta_{N}=1$). The quarkyonic model
parameters are fixed to be $\Lambda_{\mathrm{Qyc}}=250$ MeV and
$\kappa_{\mathrm{Qyc}}=0.3$. The predictions for the quarkyonic phase (green
lines) are confronted to the ones for the pure nucleon phase (magenta lines).
The model predictions are stopped when causality is violated. The sound
velocity in quarkyonic matter has a peak at around $n_{B}\approx 0.4$ fm-3, as
shown in Ref. McLerran:2019  for SM and NM, and confirmed here for AM. The
position of the peak is almost independent of the isospin asymmetry, but the
peak is a bit more pronounced in NM compared to SM. For the chosen parameters,
the sound speed predicted by the quarkyonic model at high density reaches a
value close to the conformal limit, i.e. $1/3$. The bump in the sound speed
density profile present in quarkyonic matter impacts the pressure, the
chemical potential and the binding energy. These thermodynamical quantities
are strongly increased for densities where the sound speed is maximal, and
they are softer at higher densities. The softening is such that the pressure
of quarkyonic matter crosses the pure nucleon one at high density, see Fig. 2.
The softening of the EoS is also predicted by usual construction of first
order phase transitions from nucleon to quark matter, and the interesting
feature of the quarkyonic model is the stiffening of the EoS at low densities,
where it really matters for NS, before the softening at high density. The
region of importance for NS properties coincides mostly with the densities
where the pressure is stiff. This is the reason why this model is of
particular interest for the phenomenology of compact stars.

The increase of the chemical potential in the cross-over region can also be
explained from the behavior of the nucleon Fermi momentum, which can be traced
down from Eq. (4) and re-written as,

$k_{F_{N}}^{3}=\frac{3\pi^{2}}{2}\Big{[}n_{N}+N_{c}^{3}n_{Q}\Big{]}\,,$ (26)

showing that the quark contribution to the nucleon Fermi momentum is strongly
enhanced by the factor $N_{c}^{3}$. As the density increases however, the
quark component of matter becomes more and more dominant and the softening
actually occurs.

The conclusion of this first set of results is that the generic features of
the quarkyonic model predicted by McLerran and Reddy McLerran:2019  are
preserved in our extension of the quarkyonic model for AM, and we can predict
similar features in NM with the same parameters as the one in SM. Our results
are also in qualitative agreement with the ones found by Zhao and Lattimer
Zhao:2020  as well as Jeong McLerran and Sen Jeong:2020  where different
nuclear potential were used.

##  II Quarkyonic model at $\beta$-equilibrium

We are now constructing the beta-equilibrium solutions, which describe the
ground state of dense matter existing in the core of compact stars. In cold
catalyzed NS, matter components are determined from the $\beta$-equilibrium
equations,

$\displaystyle\mu_{n}-\mu_{p}$ $\displaystyle=$ $\displaystyle\mu_{e}\,,$ (27)
$\displaystyle\mu_{e}$ $\displaystyle=$ $\displaystyle\mu_{\mu}\,,$ (28)

and charge neutrality,

$\displaystyle n_{e}+n_{\mu}+\frac{1}{3}n_{d}=n_{p}+\frac{2}{3}n_{u}\,.$ (29)

At fixed total density, these three equations allows the determination of
three variables: the isospin asymmetry $\delta_{N}$ and the electron and muon
densities, $n_{e}$ and $n_{\mu}$.

Note that the charge neutrality condition (29) in NM becomes $n_{d}=2n_{u}$,
which coincides with the relation between $u$ and $d$ quark Fermi momenta,
$k_{Fd}=2^{1/3}k_{Fu}$, imposed in NM in McLerran:2019 .

Figure 3: Particle fractions for the SLy4 nucleon model and the quarkyonic
model taking $\Lambda_{\mathrm{Qyc}}=250$ MeV and $\kappa_{\mathrm{Qyc}}=0.3$.

The particle fractions in dense matter are shown in Fig. 3 for the SLy4
nucleon model (magenta) and the quarkyonic model taking
$\Lambda_{\mathrm{Qyc}}=250$ MeV and $\kappa_{\mathrm{Qyc}}=0.3$ (green). On
the left panel of Fig. 3 are shown only the baryon contributions, $n$, $p$,
$d$ and $u$, while on the right panel, the contributions of the nucleons,
quarks, electrons and muons are represented. In the crossover region, neutron
and proton densities are reduced compare to the original nucleon model while
the amount of quarks increases. In particular, we observe that the fraction of
protons is reduced in the quarkyonic model such that it remains below the
dURCA threshold ($\gtrsim 1/9$% in the presence of muons Lattimer:1991  ).

Figure 4: Comparison of the proton fraction in the pure nucleon model with the
prediction of the quarkyonic model for various sets of the parameters
($\Lambda_{\mathrm{Qyc}}$, $\kappa_{\mathrm{Qyc}}$): ($250,0.3$), ($270,0.3$)
and ($250,0.2$). Left: function of the particle density in $\,\text{fm}^{-3}$,
right: function of the energy density $\rho_{B}c^{2}$ in MeV
$\,\text{fm}^{-3}$. Figure 5: Neutron chemical potential $\mu_{n}$, total
energy per particle $E_{tot}/A$, total pressure $P_{tot}$ and total sound
speed $(v_{s,tot}/c)^{2}$ at $\beta$-equilibrium for the parameters
($\Lambda_{\mathrm{Qyc}}$, $\kappa_{\mathrm{Qyc}}$), same choice as in Fig. 4.

In order to investigate the role of the parameters of the quarkyonic model on
the proton fraction, we show in Fig. 4 a comparison for different choices for
the parameters ($\Lambda_{\mathrm{Qyc}}$, $\kappa_{\mathrm{Qyc}}$):
($250,0.3$), ($270,0.3$) and ($250,0.2$). Increasing $\Lambda_{\mathrm{Qyc}}$
from 250 MeV to 270 MeV induces an increase of the proton fraction at high
density, which passes the dURCA threshold, but at higher density. In terms of
energy density $\rho_{B}c^{2}$, right panel in Fig. 4, the dURCA threshold is
pushed even higher, since the quarkyonic model predicts larger energy
densities than the pure nucleon one.

Finally, the thermodynamical properties of quarkyonic matter at
$\beta$-equilibrium are shown in Fig. 5: the neutron chemical potential
$\mu_{n}$, the total energy per particle $E_{tot}/A$, the equation of state
$P(\rho_{B}c^{2})$ and the total sound speed. The effect of varying the
parameters of the quarkyonic matter is also shown. The largest impact is
observed when the parameter $\Lambda_{\mathrm{Qyc}}$ in increased from 250 MeV
to 270 MeV. Increasing $\Lambda_{\mathrm{Qyc}}$ has the effect of pushing the
onset of the first quarks at higher density, as it can also be observed in
Fig. 5. As a consequence, increasing $\Lambda_{\mathrm{Qyc}}$ makes quarkyonic
matter more and more repulsive, except at low density where the larger
$\Lambda_{\mathrm{Qyc}}$ the softer the EoS, as discussed previously. The
effect of changing $\kappa_{\mathrm{Qyc}}$ is smaller. It was tuned in the
original paper by McLerran and Reddy McLerran:2019  to the conformal limit for
the sound speed.

##  III Quarkyonic stars

The structure of non-rotating neutron stars is provided by the solution of the
spherical hydrostatic equations in general relativity, also named the Tolmann-
Oppenheimer-Volkof equations TOV ,

$\displaystyle\frac{dm(r)}{dr}$ $\displaystyle=$ $\displaystyle 4\pi
r^{2}\rho(r),$ (30) $\displaystyle\frac{dP(r)}{dr}$ $\displaystyle=$
$\displaystyle-\rho(r)c^{2}\Bigg{(}1+\frac{P(r)}{\rho(r)c^{2}}\Bigg{)}\frac{d\Phi(r)}{dr},$
$\displaystyle\frac{d\Phi(r)}{dr}$ $\displaystyle=$
$\displaystyle\frac{Gm(r)}{c^{2}r^{2}}\Bigg{(}1+\frac{4\pi
P(r)r^{3}}{m(r)c^{2}}\Bigg{)}\Bigg{(}1-\frac{2Gm(r)}{rc^{2}}\Bigg{)}^{-1},$

where $G$ is the gravitational constant, $c$ the speed of light, $P(r)$ the
total pressure, $m(r)$ the enclosed mass, $\rho(r)=\rho_{B}(r)$ is the total
mass-energy density and $\Phi(r)$ the gravitational field. $P$ and $\rho$ have
contributions from both the baryons ($P_{B}$, $\rho_{B}$) and the leptons
($P_{L}$, $\rho_{L}$).

The four variables ($m$, $\rho$, $P$, $\Phi$) are obtained from the solution
of the three TOV equations (30) and the EoS for the quarkyonic matter. In the
present calculation, a crust EoS is smoothly connected to the core EoS as
described in Ref. Margueron:2018b . The tidal deformability $\Lambda_{GW}$
induced by an external quadrupole field is expressed in terms of the Love
number $k_{2}$ as $\Lambda_{GW}=2k_{2}/(3C^{5})$, where the compactness is
$C=GM/(Rc^{2})$, and $k_{2}$ is calculated from the pulsation equation at the
surface of NS Hinderer:2008 ; Flanagan:2008 ,

$\displaystyle k_{2}$ $\displaystyle=$
$\displaystyle\frac{8C^{5}}{5}\left(1-2C\right)^{2}\left(2-y_{R}+2C(y_{R}-1)\right)$
(31) $\displaystyle\times\Big{(}2C(6-3y_{R}+3C(5y_{R}-8))$
$\displaystyle+4C^{3}\left(13-11y_{R}+C(3y_{R}-2)+2C^{2}(1+y_{R})\right)$
$\displaystyle+3(1-2C)^{2}\left(2-y_{R}+2C(y_{R}-1)\right)\ln(1-2C)\Big{)}^{-1}\\!\\!\\!\\!,$

where $y_{R}$ is the value of the $y$ function at radius $R$, $y_{R}=y(r=R)$,
and $y(r)$ is the solution of the following differential equation,

$r\frac{dy}{dr}+y^{2}+yF(r)+r^{2}Q(r)=0\,,$ (32)

with the boundary condition $y(0)=2$ and the functions $F(r)$ and $Q(r)$
defined as,

$\displaystyle F(r)$ $\displaystyle=$ $\displaystyle\frac{1-4\pi
r^{2}G[\rho(r)-P(r)]/c^{4}}{1-2M(r)G/(rc^{2})}\,,$ (33) $\displaystyle
r^{2}Q(r)$ $\displaystyle=$ $\displaystyle\frac{4\pi
r^{2}G}{c^{4}}\left(5\rho(r)+9P(r)+\frac{\partial\rho(r)}{\partial
P(r)}[\rho(r)+P(r)]\right)$
$\displaystyle\times\left(1-2M(r)G/(rc^{2})\right)^{-1}-6\left(1-2M(r)G/(rc^{2})\right)^{-1}$
$\displaystyle-\frac{4G^{2}}{r^{2}c^{8}}\left(M(r)c^{2}+4\pi
r^{3}P(r)\right)^{2}\left(1-2M(r)G/(rc^{2})\right)^{-2}\,.$

The NS moment of inertia is obtained from the slow rotation approximation
Hartle:1967 ; Morrison:2004  as

$I=\frac{8\pi}{3}\int_{0}^{R}drr^{4}\rho(r)\left(1+\frac{P}{\rho(c)c^{2}}\right)\frac{\bar{\omega}}{\Omega}e^{\lambda-\Phi}\,,$
(35)

where $\bar{\omega}$ is the local spin frequency, which represents the general
relativistic correction to the asymptotic angular momentum $\Omega$ and
$\lambda$ is defined as $\exp(-2\lambda)=1-Gm/(rc^{2})$.

As usual, for a given EoS the family of solutions is parameterized by the
central density or pressure or enthalpy. The EoS are characterized by their
evolution in the mass-radius diagram, where both masses and radii of compact
stars could in principle be measured, as discussed in our introduction, see
also Ref. Watts:2015 .

Figure 6: Neutron star properties (mass $M$, radius $R$, tidal deformability
$\Lambda_{GW}$, central density $n_{c}$, binding energy $E_{bind}=M_{b}-M$,
where $M_{b}$ is the baryon mass, and the moment of inertia $I$) for
quarkyonic matter with $\Lambda_{\mathrm{Qyc}}$ and $\kappa_{\mathrm{Qyc}}$
fixed as in Fig. 4 (green curves) and nucleon matter (solid magenta curve).

We show in Fig. 6 the predictions for the mass $M$, radius $R$, tidal
deformability $\Lambda_{GW}$, central density $n_{c}$, binding energy
$E_{bind}$ and moment of inertia $I$ associated to various quarkyonic EoS with
($\Lambda_{\mathrm{Qyc}}$, $\kappa_{\mathrm{Qyc}}$): ($250,0.3$), ($270,0.3$)
and ($250,0.2$) (green lines). These predictions are confronted to the ones
for a nucleon EoS (solid magenta line).

The impact of quarkyonic matter on the mass-radius relation is huge, as
already noticed in Refs. McLerran:2019 ; Jeong:2020 . While the maximum mass
for SLy4 is reached at 2.03 M⊙, the quarkyonic stars almost reach 3 M⊙. These
is also a large impact on the radius: the 1.4 M⊙ radius $R_{1.4}$ of the pure
nucleon model is about 11.5 km, while it is pushed up to 13-14 km in the
quarkyonic model. Quarkyonic stars can therefore be much more massive than
pure nucleon ones, and are also bigger in size. Quarkyonic stars have also
different tidal deformabilities $\Lambda_{GW}$ compared to the pure nucleon
case. For the same mass, the quarkyonic stars have larger $\Lambda_{GW}$, and
at fixed $\Lambda_{GW}$, they have larger radii. At a fixed central density
(in $\,\text{fm}^{-3}$), quarkyonic stars are much more massive than the pure
nucleon model we considered. This is an effect of the repulsion observed for
the pressure in Fig. 5. For the same mass, quarkyonic stars have a slightly
lower $E_{bind}$ compared to the associated nucleonic star, they however have
a larger moment of inertia.

One can estimate the influence of the parameters $\Lambda_{\mathrm{Qyc}}$ and
$\kappa_{\mathrm{Qyc}}$. As observed in previous figures, the parameter
$\kappa_{\mathrm{Qyc}}$ is almost not influential at all, while
$\Lambda_{\mathrm{Qyc}}$ is much more critical. By increasing
$\Lambda_{\mathrm{Qyc}}$, the realization of quarkyonic matter is pushed up
and the EoS get closer to the pure nucleon case. The quarkyonic star therefore
also get a bit closer to the neutron star as $\Lambda_{\mathrm{Qyc}}$
increases. The journey in the mass-radius diagram is therefore very much
controlled by the parameter $\Lambda_{\mathrm{Qyc}}$. By decreasing
$\Lambda_{\mathrm{Qyc}}$ the quarkyonic star gets bigger and bigger compared
to the NS, as well as it gets more and more massive. In the future, the
observation of several points in mass and radius, e.g. from NICER
observations, will thus be very useful for the determination of the parameter
$\Lambda_{\mathrm{Qyc}}$.

Finally we discuss the dURCA threshold, represented in the curves by the solid
circle. Only the pure nucleon case reaches the dURCA threshold, and it happens
close to 2 M⊙. Even the quarkyonic model with ($\Lambda_{\mathrm{Qyc}}$,
$\kappa_{\mathrm{Qyc}}$)=($270,0.3$), which satisfies the dURCA condition at
high density, see Fig. 4, reaches the unstable branch before it gets to the
dURCA density. It is thus more difficult to reach the dURCA condition with
quarkyonic stars. The same conclusion was also obtained by Zhao and Lattimer
Zhao:2020  with their version of quarkyonic matter.

Figure 7: Compactness $(M/M_{\odot})/(R/\hbox{km})$ as function of the mass
$M/M_{\odot}$ (left) and radius (right) for various sets of the parameters
$\Lambda_{\mathrm{Qyc}}$ and $\kappa_{\mathrm{Qyc}}$ fixed to be as in Fig. 4.

We show in Fig. 7 the NS compactness defined as $(M/M_{\odot})/(R/\hbox{km})$,
where $R$ is expressed in km, function of the mass (left panel) and of the
radius (right panel). The compactness of the isolated NS RX J0720.4-3135 has
been extracted from observations and estimated to be 0.105$\pm$0.002
Hambaryan:2017 . Reporting this value in the left panel of Fig. 7 we observe
that the nucleonic EoS SLy4 suggests that the mass of the pulsar is 1.25
$M_{\odot}$, compatible with observed masses but close to their lower limit
Ozel:2016b , while the quarkyonic model suggest higher masses, up to 1.43
$M_{\odot}$, which is compatible with the canonical NS mass. Note that a
Bayesian exploration of nucleonic models has predicted a centroid of about
$1.33~{}M_{\odot}$ Margueron:2018b . This value is slightly larger than the
one suggested by SLy4 EoS, but is still lower that the canonical mass. On the
right panel of Fig. 7 the radii associated to the observed compactness are
also reported. While the SLy4 EoS favors 11.9 km, the quarkyonic stars point
towards larger radii, up to 13.6 km in the upper case. So we conclude that for
a fixed value of the compactness, Fig. 7 shows that quarkyonic stars have
larger masses and radii than NS.

Figure 8: Gravitational redshift $z_{\mathrm{surf}}$ as function of the mass
$M/M_{\odot}$ (left) and radius (right) for various sets of the parameters
$\Lambda_{\mathrm{Qyc}}$ and $\kappa_{\mathrm{Qyc}}$ fixed to be as in Fig. 4.

The gravitational redshift $z_{\mathrm{surf}}$ associated to the radial
emission of photons from the surface, which is detected by a distant observer
is defined as $z_{\mathrm{surf}}=\left(1-2GM/(Rc^{2})\right)^{-1/2}-1$. In
Fig. 8 we show $z_{\mathrm{surf}}$ versus the stellar mass (left panel) and
versus the stellar radius (right panel). The emission line feature of the
gamma ray burst GB790305, assumed to originate from the $e^{-}e^{+}$
annihilation (observed peak at 430 keV, line width 150$\pm$20 keV), and
assuming thermal nature of line broadening and taking due account of the
thermal blueshift, leads to the observational constraint
$z_{\mathrm{surf}}^{\mathrm{GB790305}}=0.22$ Mazets:1982 ; Higdon:1990 , which
is reported in Fig. 8. We also deduce from this observational data that SLy4
would favor typical masses of the order of 1.32 $M_{\odot}$, low masses, while
the quarkyonic star built on the same nucleonic star would point towards 1.53
$M_{\odot}$, closer to the canonical NS mass. The radius of extracted from
SLy4 would be 11.9 km, while quarkyonic star would point towards larger radii,
up to about 13.8 km.

Figure 9: Mass-radius relation for SLy4-soft nuclear interaction and for the
quarkyonic model with various sets of the parameters $\Lambda_{\mathrm{Qyc}}$
and $\kappa_{\mathrm{Qyc}}$ fixed to be as in Fig. 4.

Finally, we construct the quarkyonic model on top of a nucleonic model which
does not reach the observation constraint of about $2M_{\odot}$. To do so, we
reduce the value of $Z_{\mathrm{sym}}$ from the SLy4 nucleonic model by 300
MeV, see table 1. The nucleonic model is shown in Fig. 9 under the label
SLy4-soft (solid magenta line) while the quarkyonic models are shown for the
same three cases as before. With this example, we show that the crossover
transition to quark matter, as described by the quarkyonic approach, can bring
enough repulsion to reach large maximum mass, even if the model for the
nucleonic part cannot satisfy the observed requirement that the maximal mass
of NS should be above about $2M_{\odot}$.

##  IV Conclusions

We have proposed an extension of the original quarkyonic model from Ref.
McLerran:2019  to AM, where the original quarkyonic model for SM is recovered
as a limit. Our extension assumes (i) that the description of the quark Fermi
see and nucleon shell is globally isoscalar and (ii) that the isospin-flavor
asymmetry in the quark and nucleon phases is fixed. These assumptions root
into the concept of the quarkyonic model where nucleons result from the strong
confining force, which strength is large close to the Fermi level. The
assumption (i) allows us to smoothly connect to the quarkyonic model in SM,
and suggests a description of NM quite comparable – at least qualitatively –
to the original one suggested by McLerran and Reddy McLerran:2019 . By fixing
the isospin/flavor asymmetry in the nucleon and quark phases, assumption (ii),
the properties of isospin asymmetric quarkyonic matter can be entirely
determined from the nucleon Fermi momentum $k_{F_{N}}$ and the isospin
asymmetry $\delta_{N}$.

NS matter at $\beta$-equilibrium is then been calculated and it provides
qualitatively similar results as in the original model McLerran:2019 . It is
also in agreement with other extensions in asymmetric matter Zhao:2020 ;
Jeong:2020 , while being based on different assumptions. In our model,
quarkyonic stars are larger and heavier than the associated NS, and the
parameter $\Lambda_{\mathrm{Qyc}}$ is playing a dominant role in changing the
radius and the mass of quarkyonic stars. This result is valid even if the
nucleonic component is soft, e.g. too soft to reach $2M_{\odot}$. The proton
fraction at $\beta$-equilibrium is found to be reduced in the quarkyonic
matter, compared to the related pure nucleonic phase, which potentially quench
fast cooling – based on dURCA process – in massive compact stars. The
confrontation to a set of masses and radii, potentially obtained in future
observations like NICER or gravitational wave detections of in-spiral binary
NS, will potentially constrain $\Lambda_{\mathrm{Qyc}}$, as well as cooling
scenarios.

In the future, we aim at incorporating the quarkyonic model in systematical
comparisons to observational data in order to better understand the properties
of dense matter. Extension of the present model to finite temperature is also
on our map for the near future, as well as improving the isospin/flavor
asymmetry relation. Also, adding chiral symmetry consideration in the
quarkyonic model, taking into account the two most striking features of QCD,
will certainly be an interesting extension to study.

###### Acknowledgements.

We thank K. S. Jeong, L. McLerran, S. Reddy, and A. Schmitt for fruitful
discussions, as well as R. Stiele, R. Somasundaram, and A. Pfaff for careful
checks of our formalism. J.M., H.H. and G.C. are supported by the CNRS/IN2P3
NewMAC project, and benefit from PHAROS COST Action MP16214. All authors are
grateful to the LABEX Lyon Institute of Origins (ANR-10-LABX-0066) of the
Université de Lyon for its financial support within the program
Investissements d’Avenir (ANR-11-IDEX-0007) of the French government operated
by the National Research Agency (ANR).

## References

  * (1) L. McLerran and S. Reddy, Phys. Rev. Lett. 122, 122701 (2019).
  * (2) The Physics and Astrophysics of Neutron Stars, L. Rezzolla et al. Editors, Springer International Publishing (2018).
  * (3) P. Haensel, A. Y. Potekhin, and D. G. Yakovlev, Neutron Stars I (Springer, Berlin, 2007).
  * (4) B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017).
  * (5) J. Antoniadis, P.C.C. Freire, N. Wex et al., Science 340, 448 (2013).
  * (6) F. Özel and P. Freire, Annu. Rev. Astron. Astrophys. 54, 401 (2016).
  * (7) E. Fonseca, T.T. Pennucci, J.A. Ellis, et al., Astro. J 832, 167 (2016).
  * (8) M. Linares, T. Shahbaz, and J. Casares, Astrophys. J. 859, 54 (2018).
  * (9) H. T. Cromartie et al., Nat. Astron. 4, 72 (2020).
  * (10) Z. Arzoumanian, A. Brazier, S. Burke-Spolaor et al., Astro. J. SS 235, 37 (2018).
  * (11) T. Hinderer, Astrophys. J. 677, 1216 (2008).
  * (12) É.É. Flanagan and T. Hinderer, Phys. Rev. D 77, 021502 (2008).
  * (13) M.W. Coughlin, T. Dietrich, B. Margalit, and B.D. Metzger, MNRAS 489, L91 (2019).
  * (14) S. Guillot, M. Servillat, N.A. Webb, & R.E. Rutledge, Astrophys. J 772, 7 (2013).
  * (15) A.W. Steiner, C.O. Heinke, S. Bogdanov, et al., MNRAS 476, 421 (2018).
  * (16) N. Baillot d’Etivaux, S. Guillot, J. Margueron, N. A. Webb, M. Catelan, A. Reisenegger, Astrophys. J. 887, 48 (2019).
  * (17) C. A. Raithel, F. Özel and D. Psaltis, Astrophys. J. 831, 44 (2016).
  * (18) K. Gendreau and Z. Arzoumanian, Nat. Astron. 1, 895 (2017).
  * (19) J. Margueron, R. Casali Hoffmann, F. Gulminelli, Phys. Rev. C 97, 025805 (2018).
  * (20) L. McLerran and R.D. Pisarski, Nucl. Phys. A 796, 83 (2007).
  * (21) K. Fukushima, and T. Kojo, Astrophys. Journ. 817, 180 (2016).
  * (22) N. Kovenky and A. Schmitt, Journ. High En. Phys. 09, 112 (2020).
  * (23) J. Zdunik and P. Haensel, Astron. & Astrophys.551, A61 (2013).
  * (24) M.G. Alford, S. Han, Phys. Rev. D 88, 083013 (2013); M.G. Alford, G.F. Burgio, S. Han, G. Taranto, and D. Zappalà, Phys. Rev. D 92, 083002 (2015).
  * (25) P. Bedaque and A.W. Steiner, Phys. Rev. Lett. 114, 031103 (2015).
  * (26) I. Tews, J. Carlson, S. Gandolfi, and S. Reddy, Astrophys. J. 860, 149 (2018).
  * (27) T. Zhao and J. L. Lattimer, Phys. Rev. D 102, 023021 (2020).
  * (28) J. L. Lattimer, C. Pethick, M. Prakash, and P. Haensel, Phys. Rev. Lett. 66, 2701 (1991).
  * (29) K. S. Jeong, L. McLerran, and S. Sen, Phys. Rev. C 101, 035201 (2020).
  * (30) D. C. Duarte, S. Hernandez-Ortiz, and K. S. Jeong, Phys. Rev. C 102, 025203 (2020).
  * (31) D. C. Duarte, S. Hernandez-Ortiz, and K. S. Jeong, Phys. Rev. C 102, 065202 (2020).
  * (32) P. G. De Gennes, Superconductivity Of Metals And Alloys. Advanced Book Classics. Advanced Book Program, (Perseus Books, 1999).
  * (33) J. Margueron, unpublished.
  * (34) J. Margueron, R. Casali Hoffmann, F. Gulminelli, Phys. Rev. C 97, 025806 (2018).
  * (35) R. C. Tolman, Phys. Rev. 55, 364 (1939); J. R. Oppenheimer and G. M. Volkoff, Phys. Rev. 55, 374 (1939).
  * (36) J. B. Hartle and D. H. Sharp, Astrophys. J. 147, 317 (1967).
  * (37) I. A. Morrison, T.W. Baumgarte, S. L. Shapiro, and V. R. Pandharipande, Astrophys. J. 617, L135 (2004).
  * (38) A. Watts et al., PoS 215 (AASKA14), 043 (2015).
  * (39) V. Hambaryan, V. Suleimanov, F. Haberl et al., Astron. Astrophys. 601, A108 (2017).
  * (40) F. Özel and P. Freire, Ann. Rev. Astron. Astrophys. 54, 401 (2016).
  * (41) E. P. Mazets, S. V. Golenetskii, Yu. A. Guryan, and V. N. Ilyinskii, Astrophy. Space Science 84, 173 (1982).
  * (42) J. C. Higdon, and R. E. Lingenfelter, Ann. Rev. Astron. Astrophys. 28, 401 (1990).