11institutetext: Department of Computer Science, University of Georgia
11email<EMAIL_ADDRESS>22institutetext: Department of Computer
Science and Engineering, Shanghai Jiao Tong University
22email<EMAIL_ADDRESS>
# DDR-Net: Dividing and Downsampling Mixed Network for Diffeomorphic Image
Registration

Ankita Joshi 11 Yi Hong 22

###### Abstract

Deep diffeomorphic registration faces significant challenges for high-
dimensional images, especially in terms of memory limits. Existing approaches
either downsample original images, or approximate underlying transformations,
or reduce model size. The information loss during the approximation or
insufficient model capacity is a hindrance to the registration accuracy for
high-dimensional images, e.g., 3D medical volumes. In this paper, we propose a
Dividing and Downsampling mixed Registration network (DDR-Net), a general
architecture that preserves most of the image information at multiple scales.
DDR-Net leverages the global context via downsampling the input and utilizes
the local details from divided chunks of the input images. This design reduces
the network input size and its memory cost; meanwhile, by fusing global and
local information, DDR-Net obtains both coarse-level and fine-level alignments
in the final deformation fields. We evaluate DDR-Net on three public datasets,
i.e., OASIS, IBSR18, and 3DIRCADB-01, and the experimental results demonstrate
our approach outperforms existing approaches. Codes are available –here–.

###### Keywords:

Diffeomorphic image registration Dividing and downsampling Multi-scale
registration.

##  1 Introduction

Deformable image registration establishes pixel- or voxel-level dense
correspondences for 2D or 3D image pairs, which form a deformation that
transforms images into a common space for comparison and analysis. Such a
deformation desires a good property of diffeomorphism, a smooth transformation
with a smooth inverse, to ensure the preservation of topology when warping
images. Classical image registration models, e.g., LDDMM [6], Stationary
Velocity Fields (SVF) [2], successfully estimate diffeomorphic deformations
for building dense correspondences between image pairs. However, these models
face challenges for practical applications, i.e., providing both fast and
accurate solutions. Therefore, researchers have been working on improving the
efficiency of diffeomorphic image registration methods [3, 17].

Recently, deep learning based approaches open an alternative to address the
above challenges, which motivates our work in this paper. Existing deep
registration models focus on tackling the efficiency challenge using
supervised [16] or unsupervised techniques [7, 12]. Supervised approaches [16]
maintain the diffeomorphic property, which is inherited from the classical
diffeomorphic model LDDMM, but it requires extra effort to obtain the ground-
truth deformations. Meanwhile, its registration accuracy is limited by that of
the obtained deformations. The unsupervised approaches [7, 12] have shown
promising diffeomorphic and efficient registration results by introducing an
integration layer into the network design, based on the scaling and squaring
method [9]. Due to the flexibility in selecting the network architecture and
the loss function, the unsupervised framework has the potential to further
improve the registration accuracy. However, because the integration step is
computationally expensive and the network faces the memory challenge for high-
dimensional images, the unsupervised approaches often work on downsampled
images or deformations, which does not fully leverage available information
and limits the registration accuracy.

Figure 1: Architecture of our proposed DDR-Net. Given an image pair $I_{0}$
and $I_{1}$, the network estimates deformations at local scale, i.e., chunk
branch, the global scale, i.e. the downsampled branch, and the original scale.
The deformation $\phi$ at each level is driven by the corresponding velocity
field $v$, which is sampled from the U-Net outputs, the mean $\mu$ and the
variance $\Sigma$. The deformation at the original resolution is obtained by
merging velocity fields generated by the global and local branches and then
used to warp the original source images to the corresponding target images by
using an interpolation layer. The number of the filters used are under the
blocks.

In this paper, we aim to improve the accuracy of unsupervised image
registration, using a multi-scale design to integrate image deformations at
global, local, and original scales. The difference in our work from current
unsupervised registration models [7, 12] is that our model works on
downsampled images and chopped chunks, which reduce the computational and
memory cost compared to working on the original image size directly.
Meanwhile, the integration of these two-scale deformations back to the
original scale improves the registration accuracy because of the information
fusion at different levels. Previous work [8, 11, 16] either use multi-scale
approaches of downsampling the entire image to improve accuracy of the result
or use only patches of the images to reduce the memory cost, but both these
techniques do not adequately leverage the data. Therefore, to gain an accuracy
boost but not run into memory issues, we propose the Dividing and Downsampling
mixed Registration Network (DDR-Net).

Our contributions in this paper are summarized as follows:

  * •

We propose a novel architecture DDR-Net, which can effectively use both global
and local features yielding high quality registration performance. Multi-scale
information benefits the task of deep image registration.

  * •

We demonstrate an effective way to obtain a trade-off between fully leveraging
the available data under limited computing resources and gaining an improved
accuracy of diffeomorphic image registration at the same time.

  * •

We conduct extensive experiments on both 2D and 3D datasets with different
image types, including brain MRIs and liver CT scans. The experimental results
demonstrate that our framework has better registration performance compared to
deep-learning-based method VoxelMorph [7] and the classical registration
method ANTs SyN [5], in terms of image matching, deformation smoothness, and
multi-structure segmentation.

##  2 Dividing and Downsampling mixed Registration Network (DDR-Net)

Architecture Overview. As shown in Figure 1, our proposed DDR-Net includes
three main components: a global branch that handles the registration for
downsampled images, a local branch that handles the registration for the
cropped local chunks of original images, and an original branch that merges
estimated velocity fields from the global and local branches to register
original images. Each branch outputs a deformation field $\phi$ to register
image pairs at its corresponding level, and they share some network designs as
discussed below.

Backbone Registration. At each scale, we have a diffeomorphic image
registration problem. Given an image pair, a source image $I_{0}$ and a target
image $I_{1}$, each of size $n_{x}\times n_{y}\times n_{z}$, the goal of
diffeormorphic image registration is to estimate a smooth deformation field
$\phi:\mathbb{R}^{n_{x}\times n_{y}\times
n_{z}}\rightarrow\mathbb{R}^{n_{x}\times n_{y}\times n_{z}}$ with a smooth
$\phi^{-1}$, such that the image deformed from the source, i.e. $\phi\cdot
I_{0}$, is similar to the target image $I_{1}$. Such a diffeomorphic
deformation field is driven by a smooth velocity field $v_{t},t\in[0,1]$, via
the following differential equation:

$\frac{d}{dt}\phi=v_{t}\circ\phi_{t},\quad\phi_{0}=id.$ (1)

Here, $id$ is an identity deformation. This formulation estimates an optimal
velocity field $v$ that drives a deformation field $\phi$ to match an image
pair. So, the registration network has three sub-tasks, i.e., estimating the
velocity field, solving Eq. (1) for deformations, and deforming an image with
interpolation.

Velocity Field Estimation. The global and local branches follow the same UNet
[14] architecture as shown in Fig. 1. The UNet takes in image pairs and
outputs the mean $\mu$ and the variance $\Sigma$ for sampling a corresponding
stationary velocity field $v$. Here, the stationary velocity field assumption
simplifies the solution of Eq. (1). Given a collection of image pairs
$\\{(I_{0},I_{1})\\}$, where $I_{0},I_{1}\in\mathbb{R}^{n_{x}\times
n_{y}\times n_{z}}$, the downsampling branch takes the low-resolution image
pairs $\\{(I_{0D},I_{1D})\\}$, downsampled by half, i.e.,
$I_{0D},I_{1D}\in\mathbb{R}^{\frac{n_{x}}{2}\times\frac{n_{y}}{2}\times\frac{n_{z}}{2}}$.
The chunk branch receives each input as $k$ divided patches with the same
resolution as $I_{0D}$ and $I_{1D}$, i.e., $k\times\\{(I_{0C},I_{1C})\\}$ and
$I_{0C},I_{1C}\in\mathbb{R}^{\frac{n_{x}}{2}\times\frac{n_{y}}{2}\times\frac{n_{z}}{2}}$.
These cropped chunks have small overlaps on their boundaries to mitigate the
discontinuity of the generated velocity fields at the chunk boundaries, when
merging back to form the original high resolution velocity fields. The
detailed network architecture in terms of the number of convolution layers and
the kernel sizes are shown in Fig. 1.

Deformation Integration. Under that assumption of the stationary velocity
field, Equation (1) is simplified with a constant velocity field, and its
solution is $\phi=e^{v}$. Similar to VoxelMorph, we adopt the scaling and
squaring algorithm [9] to approximate this solution, which is implemented as a
differentiable layer in the network. The downsampling and chunk branches
integrate deformations separately, which are driven by their respective
velocity fields. In the original branch, we use an averaged global and local
velocity field to integrate the deformation. In particular, we upsample the
global velocity field to the original resolution using an upsampling layer and
merge the $k$ chunks to the original volume using a chunk-merge layer. As a
result, we obtain a velocity field at the original resolution which is
integrated to get the deformation at the original level.

Image Interpolation. We use an interpolation layer to deform the source image
at each branch. For each voxel $p$ in the target image, we compute its
location $\phi(p)$ in the source image and compute its intensity value using
linear interpolation. This differentiable operation allows the backpropagation
of the network errors.

Loss Functions. All the velocity fields generated in the network have the
constraint to ensure the smoothness of the deformations. Each branch also
outputs a deformed source image to match the corresponding target image at
each level, i.e., the downsampled, chunk, and original scales. We use the Mean
Square Error (MSE) to measure the matching results and a K-L divergence loss
to encourage the smoothness of the velocity field as done in [7, 11].

Image Pair | $\phi\cdot I_{0}$ | $\phi\cdot I_{0}-I_{1}$ | $\phi$ | $det(J_{\phi})$ |
---|---|---|---|---|---
| $I_{0}$
---
$I_{1}$
|
---
|
---
|
---
|
---
|  SyN
---
VM
Ours
Figure 2: Registration comparison among SyN, VM (VoxelMorph), and our DDR-Net. Left to right: the image pair median slice from OASIS dataset, the source image $I_{0}$ and the target image $I_{1}$; the warped image; the image difference between warped image and the target image; the deformation $\phi$; the determinant of the deformation Jacobian. Image Pair | $\phi\cdot I_{0}$ | $\phi\cdot I_{0}-I_{1}$ | $\phi$ | $det(J_{\phi})$ |
---|---|---|---|---|---
| $I_{0}$
---
$I_{1}$
|
---
|
---
|
---
|
---
|  SyN
---
VM
Ours
Figure 3: Registration comparison among SyN, VM (VoxelMorph), and our DDR-Net for the IBSR dataset. Left to right: the same with Fig. 2. Image Pair | $\phi\cdot I_{0}$ | $\phi\cdot I_{0}-I_{1}$ | $\phi$ | $det(J_{\phi})$ |
---|---|---|---|---|---
| $I_{0}$
---
$I_{1}$
|
---
|
---
|
---
|
---
|  SyN
---
VM
Ours
Figure 4: Registration comparison among SyN, VM (VoxelMorph), and our DDR-Net for the 3DIRCADB-01 dataset. Left to right: the same with Fig. 2. $I_{0}$ | $I_{1}$ |  $\phi\cdot I_{0}$ (SyN) |  $\phi\cdot I_{0}$ (VM) |  $\phi\cdot I_{0}$ (Ours)
---|---|---|---|---
|  |  |  |
|  |  |  |
Figure 5:  IBSR18 sampled source, target, and warped images for SyN, VoxelMorph (VM), and our DDR-Net. Blue: ventricles, green: thalami, red: hippocampi. Dataset | Method | RMSE ($e^{-3}$) | Foldings
---|---|---|---
$|\sum det(J_{\phi})<0|$ |  $\frac{1}{N}\sum\delta(det(J_{\phi})<0)$
$(\mathbf{e^{-5}})$ |  (Ratio: $\mathbb{\permil}$)
OASIS |  SyN [5] |  4.79 $\pm$ 0.001 |  495.37 $\pm$ 0.03 |  0.1 $\pm$ 0.768
VM [7] |  1.25 $\pm$ 0.001 |  1.46 $\pm$ 0.005 |  0.051 $\pm$ 0.014
DDR-Net | 1.10 $\pm$ 0.001 | 1.09 $\pm$ 0.000 | 0.003 $\pm$ 0.019
IBSR18 |  SyN [5] |  10.93 $\pm$ 0.26 | 0.00 $\pm$ 0.000 | 0.00 $\pm$ 0.00
VM [7] |  2.06 $\pm$ 0.001 |  20.17 $\pm$ 0.005 |  0.05 $\pm$ 0.014
DDR-Net | 1.67 $\pm$ 0.001 |  17.266 $\pm$ 0.000 |  0.04 $\pm$ 0.117
3DIRCADB-01 |  SyN [5] |  68.04 $\pm$ 0.024 | 476.98 $\pm$ 0.001 |  35.32 $\pm$ 5.227
VM [7] |  16.84 $\pm$ 0.005 |  1221.21 $\pm$ 0.01 |  0.72 $\pm$ 0.274
DDR-Net | 14.11 $\pm$ 0.004 |  723.89 $\pm$ 0.003 | 0.65 $\pm$ 0.247
Table 1: Comparison of the registration performance on all the three datasets. Dataset | Region | SyN [4] | VM [7] | DDR-Net
---|---|---|---|---
IBSR18 | 3rd ventricle | 0.28 $\pm$ 0.14 | 0.42 $\pm$ 0.16 | 0.44 $\pm$ 0.17
4th ventricle | 0.20 $\pm$ 0.13 | 0.39 $\pm$ 0.16 | 0.36 $\pm$ 0.16
amygdala | 0.34 $\pm$ 0.13 | 0.27 $\pm$ 0.23 | 0.31 $\pm$ 0.20
brainstem | 0.62 $\pm$ 0.12 | 0.69 $\pm$ 0.13 | 0.66 $\pm$ 0.16
caudate | 0.48 $\pm$ 0.09 | 0.40 $\pm$ 0.21 | 0.40 $\pm$ 0.20
cerebellum cortex | 0.55 $\pm$ 0.12 | 0.59 $\pm$ 0.20 | 0.63 $\pm$ 0.13
cerebellum white matter | 0.42 $\pm$ 0.14 | 0.34 $\pm$ 0.20 | 0.44 $\pm$ 0.18
cerebral cortex | 0.43 $\pm$ 0.09 | 0.55 $\pm$ 0.18 | 0.58 $\pm$ 0.13
cerebral white matter | 0.53 $\pm$ 0.06 | 0.58 $\pm$ 0.15 | 0.59 $\pm$ 0.11
csf | 0.22 $\pm$ 0.17 | 0.31 $\pm$ 0.16 | 0.31 $\pm$ 0.17
hippocampus | 0.32 $\pm$ 0.11 | 0.17 $\pm$ 0.14 | 0.34 $\pm$ 0.21
lateral ventricle | 0.38 $\pm$ 0.12 | 0.53 $\pm$ 0.18 | 0.45 $\pm$ 0.17
pallidum | 0.24 $\pm$ 0.16 | 0.23 $\pm$ 0.23 | 0.29 $\pm$ 0.24
putamen | 0.29 $\pm$ 0.21 | 0.27 $\pm$ 0.25 | 0.33 $\pm$ 0.27
thalamus | 0.68 $\pm$ 0.08 | 0.61 $\pm$ 0.18 | 0.60 $\pm$ 0.19
ventraldc | 0.49 $\pm$ 0.11 | 0.49 $\pm$ 0.21 | 0.54 $\pm$ 0.18
Avg. Dice | 0.28 $\pm$ 0.14 | 0.43 $\pm$ 0.1 | 0.47 $\pm$ 0.1
3DIRCADB-01 | bone | 0.223 $\pm$ 0.106 | 0.356 $\pm$ 0.149 | 0.369 $\pm$ 0.153
skin | 0.889 $\pm$ 0.065 | 0.982 $\pm$ 0.012 | 0.991 $\pm$ 0.009
liver | 0.690 $\pm$ 0.099 | 0.731 $\pm$ 0.122 | 0.741 $\pm$ 0.131
Table 2: Segmentation comparison using SyN, VoxelMorph (VM) and DDR-Net.

##  3 Experiments

We evaluate our method on three public datasets, which involves two types of
3D medical images, i.e., brain MRI scans and liver CT scans.

OASIS Dataset [13]. The T1-weighted brain scans from the OASIS dataset
underwent pre-processing steps including down-sampling, skull-stripping,
intensity normalization to the range [0, 1], bias field correction, and co-
registration with affine transformations. We have resulting images of
resampled resolution $128\times 128\times 128$ and the voxel size of
$1.25\times 1.25\times 1.25mm^{3}$. We collect 360 volume pairs for the 3D
experiments.

IBSR18 Dataset [15]. The IBSR18 dataset consists of T1-weighted scans for 18
subjects with dimensions of $256\times 128\times 256$. We resampled these
images to $128\times 128\times 128$. Scans underwent preprocessing steps
including skull-stripping, bias field correction, and intensity normalization.
Each scan also comes with 84 manually labelled anatomical structures. We
collect 306 3D image pairs. Segmentation maps including 28 anatomical
structures (see Table 2) were also obtained in this dataset, which followed
the same preprocessing steps of resampling.

3DIRCADB-01 Dataset [1]. The 3DIRCADB-01 dataset contains 20 CT scans with
masks of the segmented structures available for bone, liver, and skin. Scans
underwent preprocessing steps including histogram equalization and intensity
normalization to the range [0, 1]. We collect and pair the central slices on
the axial plan and obtained 380 2D image pairs with a resolution of $512\times
512$.

Experimental Settings. For all the experiments, we use 70% of the collected
data for training, 10% for validation, and 20% for test, after a random
shuffle. We use the Adam [10] optimizer with a learning rate of $e^{-4}$. All
the experiments were trained using NVIDIA GeForce TITAN X GPUs.

We measure the image matching after registration using the intensity root mean
square error (RMSE), and the smoothness of a deformation by counting the
number of its foldings and the absolute sum of negative determinant of its
Jacobian [3]. Also, we perform the segmentation task using image registration
and report the Dice score in terms of volume overlap.

We compare our approach with ANTsPy package using Symmetric Normalization
(SyN) [5]. We use cross-correlation and other default settings, which are
optimal for our task with 201 iterations for registration. Another baseline
algorithm is VoxelMorph [7], and we also use its default settings for
comparison. Initial affine registration is not performed for any of the
experiments.

Experimental Results. Figures 2, 3, 4 show qualitative results for OASIS,
IBSR18, and 3DIRCADB-01 datasets, respectively. Overall, our approach produces
better matching results as compared to the baseline methods, demonstrated by
the image difference plots. VoxelMorph tends to produce unwanted artifacts in
its deformed images, which are not seen in our results. The deformation plot
and the visualization of the determinant of the Jacobian show that our results
have less deformations happening in the background, as expected. Also, our
method achieves significantly smaller deformations throughout the image space
and provides a better matching result to the target at the same time.

We observe that the VoxelMorph produces more irregular deformation fields
compared to the other methods. Table 1 show the quantitative analysis of the
registration performance on the OASIS, IBSR18, and 3DIRCADB-01 datasets. We
achieve consistently lower number of foldings compared to the baseline methods
while maintaining the consistently higher registration accuracy. Table 2
indicates the Dice score of anatomical structures for the IBSR18 and
3DIRACADB-01 datasets. Our method performs better on most of the anatomical
structures for both the IBSR18 and 3DIRCADB-01 datasets. Figure 5 shows the
segmentation maps for a few anatomical structures from the IBSR18 dataset.

Regarding the inference time, ANTs SyN, which is tested on CPU since it does
not have a GPU implementation, takes on an average 20 minutes, whereas
VoxelMorph tested on GPU takes 121 milliseconds, while DDR-Net tested on GPU
takes 425 milliseconds to register an image pair from the OASIS dataset. A
limited ablation study conducted by us revealed that a single scale UNet
architecture working on the full size velocity field integration for 2D input
size of $128\times 128$ for a single input image pair takes 26.46 MB of
memory, whereas an architecture like DDR-Net for the same input and UNet
architecture will take 21.16 MB of memory. The memory requirement was
calculated considering the forward pass and parameters in the neural network.
For 3D cases, our DDR-Net can take in more images in a batch, compared to a
single scale UNet.

##  4 Conclusion and Discussion

In this paper, we have proposed a diffeomorphic image registration model which
estimates smoother deformations and provides a better image matching result.
It leverages both the global context and local fine structures of the data
effectively to enhance the registration result and handle large deformations
as a result of such an architecture. Generally, our approach produced more
regular deformation fields, which are significantly smoother than the baseline
methods. The dice scores indicate that our approach is comparable to the
existing methods even surpassing them in most cases. Moreover, our
architecture shows that simple yet efficient changes in architectures can lead
to better deep learning strategies. We believe that an optimal balance of GPU
memory and accuracy is essential for registering image pairs in 3D in order to
maximize the utilization of 3D information. Our work will push new avenues of
research in considering memory efficient registration models in deep learning,
which fully leverage the data.

## References

  * [1] 3dircadb-01. https://www.ircad.fr/research/3dircadb/
  * [2] Arsigny, V., Commowick, O., Pennec, X., Ayache, N.: A log-euclidean framework for statistics on diffeomorphisms. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 924–931. Springer (2006)
  * [3] Ashburner, J.: A fast diffeomorphic image registration algorithm. Neuroimage 38(1), 95–113 (2007)
  * [4] Avants, B.B., Epstein, C.L., Grossman, M., Gee, J.C.: Symmetric diffeomorphic image registration with cross-correlation: evaluating automated labeling of elderly and neurodegenerative brain. Medical image analysis 12(1), 26–41 (2008)
  * [5] Avants, B.B., Tustison, N.J., Song, G., Cook, P.A., Klein, A., Gee, J.C.: A reproducible evaluation of ants similarity metric performance in brain image registration. Neuroimage 54(3), 2033–2044 (2011)
  * [6] Beg, M.F., Miller, M.I., Trouvé, A., Younes, L.: Computing large deformation metric mappings via geodesic flows of diffeomorphisms. International journal of computer vision 61(2), 139–157 (2005)
  * [7] Dalca, A.V., Balakrishnan, G., Guttag, J., Sabuncu, M.R.: Unsupervised learning for fast probabilistic diffeomorphic registration. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 729–738. Springer (2018)
  * [8] Hering, A., van Ginneken, B., Heldmann, S.: mlvirnet: Multilevel variational image registration network. In: International Conference on Medical Image Computing and Computer-Assisted Intervention. pp. 257–265. Springer (2019)
  * [9] Higham, N.J.: The scaling and squaring method for the matrix exponential revisited. SIAM Journal on Matrix Analysis and Applications 26(4), 1179–1193 (2005)
  * [10] Kingma, D.P., Ba, J.: Adam: A method for stochastic optimization. arXiv preprint arXiv:1412.6980 (2014)
  * [11] Krebs, J., Delingette, H., Mailhé, B., Ayache, N., Mansi, T.: Learning a probabilistic model for diffeomorphic registration. IEEE transactions on medical imaging 38(9), 2165–2176 (2019)
  * [12] Krebs, J., Mansi, T., Mailhé, B., Ayache, N., Delingette, H.: Unsupervised probabilistic deformation modeling for robust diffeomorphic registration. In: Deep Learning in Medical Image Analysis and Multimodal Learning for Clinical Decision Support, pp. 101–109. Springer (2018)
  * [13] Marcus, D.S., Wang, T.H., Parker, J., Csernansky, J.G., Morris, J.C., Buckner, R.L.: Open access series of imaging studies (oasis): cross-sectional mri data in young, middle aged, nondemented, and demented older adults. Journal of cognitive neuroscience 19(9), 1498–1507 (2007)
  * [14] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: International Conference on Medical image computing and computer-assisted intervention. pp. 234–241. Springer (2015)
  * [15] Valverde, S., Oliver, A., Cabezas, M., Roura, E., Lladó, X.: Comparison of 10 brain tissue segmentation methods using revisited ibsr annotations. Journal of Magnetic Resonance Imaging 41(1), 93–101 (2015)
  * [16] Yang, X., Kwitt, R., Styner, M., Niethammer, M.: Quicksilver: Fast predictive image registration–a deep learning approach. NeuroImage 158, 378–396 (2017)
  * [17] Zhang, M., Fletcher, P.T.: Finite-dimensional lie algebras for fast diffeomorphic image registration. In: International conference on information processing in medical imaging. pp. 249–260. Springer (2015)


# Quark formation and phenomenology in binary neutron-star mergers using V-QCD

Samuel Tootle  Christian Ecker  Konrad Topolski  Institut für Theoretische
Physik, Goethe Universität, Max-von-Laue-Str. 1, 60438 Frankfurt am Main,
Germany  Tuna Demircik  Asia Pacific Center for Theoretical Physics, Pohang,
37673, Korea  Matti Järvinen  Asia Pacific Center for Theoretical Physics,
Pohang, 37673, Korea  Department of Physics, Pohang University of Science and
Technology, Pohang, 37673, Korea  Luciano Rezzolla  Institut für Theoretische
Physik, Goethe Universität, Max-von-Laue-Str. 1, 60438 Frankfurt am Main,
Germany  Frankfurt Institute for Advanced Studies, Ruth-Moufang-Str. 1, 60438
Frankfurt, Germany  School of Mathematics, Trinity College, Dublin 2, Ireland

###### Abstract

Using full 3+1 dimensional general-relativistic hydrodynamic simulations of
equal- and unequal-mass neutron-star binaries with properties that are
consistent with those inferred from the inspiral of GW170817, we perform a
detailed study of the quark-formation processes that could take place after
merger. We use three equations of state consistent with current pulsar
observations derived from a novel finite-temperature framework based on V-QCD,
a non-perturbative gauge/gravity model for Quantum Chromodynamics. In this
way, we identify three different post-merger stages at which mixed baryonic
and quark matter, as well as pure quark matter, are generated. A phase
transition triggered collapse already $\lesssim 10\,\rm{ms}$ after the merger
reveals that the softest version of our equations of state is actually
inconsistent with the expected second-long post-merger lifetime of GW170817.
Our results underline the impact that multi-messenger observations of binary
neutron-star mergers can have in constraining the equation of state of nuclear
matter, especially in its most extreme regimes.

neutron star mergers, critical point, gauge/gravity duality

††preprint: APCTP Pre2022 - 007

##  I Introduction

Multi-messenger observations from binary neutron-star mergers are expected to
provide new insights into the properties of dense and hot Quantum
Chromodynamics (QCD). However, due to the lack of first principle techniques,
the precise phase structure of QCD and its equation of state (EOS) are
currently not known at densities and temperatures that occur during and after
the merger of two neutron stars. It is, therefore, not clear if a transition
from dense nuclear to quark matter can happen in such events nor if such a
transition leaves a discernible imprint on the gravitational waveform or the
electromagnetic counterpart. Such an imprint may be best visible in the kHz
gravitational-wave signal emitted during a possible post-merger hypermassive
neutron star (HMNS) stage as it encodes information on the hot and dense part
of the EOS that is inaccessible during the inspiral phase, however, this
parameter space will only be able to be reliably explored by future detectors.
In order to resolve the question of the detectability of quark matter in
binary neutron-star mergers, it is therefore crucial to explore state-of-the-
art models that make predictions for the phase structure of QCD beyond the
saturation density of atomic nuclei $n_{s}=0.16\,\rm fm^{-3}$ and, at the same
time, satisfy known constraints from theory and astrophysical observations.

Numerical simulations of binary neutron-star systems play a critical role to
ascertain the influence of the neutron-star characteristics (mass and spin)
and the EOS on the post-merger remnant and lifetime [1, 2]. For EOSs that
include a transition to quark matter, these simulations also provide insights
into the mechanisms for and the abundance of quark formation in the post-
merger remnant. There exists a number of works that study mergers using models
with a phase transition either in an effective particle approach [3] or in
fully general-relativistic hydrodynamics [4, 5, 6]. The recent study [7] also
analyses the composition of the merger ejecta. In [8] possible signatures of
phase transitions in the gravitational wave signal have been classified using
a polytropic EOS with the phase transition to quark matter introduced via a
Gibbs-like construction and temperature dependence using the standard
$\Gamma_{\rm th}$ component based on an ideal-fluid EOS [9].

In this article we perform the first study to analyse quark matter production
in BNS mergers using a novel framework [10] for describing nuclear matter
which combines the gauge/gravity duality with nuclear theory models. In this
combined framework, the EOS at large densities is given by the gauge/gravity
duality [11, 12], more precisely the V-QCD model [13, 14, 15], which is used
to describe the deconfinement phase transition from dense baryonic to quark
matter. The DD2 version of the Hempel-Schaffner-Bielich statistical model [16,
17] is used at low densities, and a van der Waals model for the temperature
dependence in the dense nuclear matter phase. The advantage of this combined
framework is that it is, by construction, consistent with theoretic
predictions from nuclear theory at low densities and perturbative QCD at high
densities. At the same time, it provides a description for the deconfinement
phase transition within a single model, i.e., V-QCD. The framework also yields
predictions for the location of the QCD critical point which allows us to
investigate the possibility that neutron star matter during and/or after the
merger reach temperatures sufficient to probe this part of the QCD phase
diagram.

Figure 1: Left: Cold, beta equilibrium slices of the V-QCD EOS. Shown are the
Soft (red), Inter (green) and Stiff (blue) versions of the model. Dark blue
and orange bands mark the uncertainties in nuclear theory [18] and
perturbative QCD [19] calculations, respectively. The light blue region marks
the theoretically allowed region spanned by polytropic EOSs, where the pink
region marks the parametric freedom of V-QCD after imposing constraints from
NS observations. Right: Corresponding mass-radius relations where dashed lines
mark unstable quark matter branches. In addition we show various error bands
of direct mass [20, 21] and radius [22, 23, 24, 25, 26] measurements of heavy
pulsars; shown with a black dashed line is a theoretical lower bound for the
radii as computed from considerations on the threshold mass by Koeppel et al.
[27].

When taking into account different individual masses and spins of the stars as
well as variations of the EOS, the parameter space of binary neutron star
simulations grows rapidly. As a first test of the novel V-QCD EOS, we study
configurations motivated by information gained from the inspiral part of the
gravitational-wave event GW170817 [28] in order to determine whether or not
the constructed EOSs agree with the observation. To do so, we perform full 3+1
dimensional general-relativistic hydrodynamics (GRHD) simulations consistent
with GW170817 so as to determine mechanisms for quark production as well as
ascertain the importance of multi-messenger detections on constraining the
model. More specifically, we simulate non-spinning binaries of fixed chirp
mass $\mathcal{M}_{\rm
chirp}:={(M_{1}M_{2})^{3/5}}/{(M_{1}+M_{2})^{1/5}}=1.186\,M_{\odot}$ and two
different mass ratios $q:=M_{2}/M_{1}=1,\,0.7$. We use three different EOSs,
based on a soft, an intermediate and a stiff version of V-QCD, whose cold
beta-equilibrium parts satisfy [29, 30, 31] all currently known observational
constraints from direct mass [32, 20, 21] and radius [25, 26, 22, 23]
measurements of heavy pulsars as well as the constraint on the tidal
deformability obtained from GW170817 [33].

As an important outcome of our simulations we identify three different stages
of quark matter formation during the post-merger HMNS evolution which will be
referenced as hot-quark (HQ), warm-quark (WQ) and cold-quark (CQ) stages. The
basis and definition of these stages are discussed in detail in Section IV.
Furthermore, we provide an analysis of the gravitational waveforms and their
frequency spectra of all our simulations. The clearest imprint on the waveform
is their early termination caused by the phase-transition-triggered collapse
to a black hole in the CQ stage. From their frequency spectra we extract a
number of characteristic post-merger frequencies.

The rest of the article is structured as follows. In Section II we provide
more details about the V-QCD model. In Section III we explain our numeric
setup for the binary initial data and the subsequent time evolution. In
Section IV we present our results for the merger dynamics and show examples
for the three post-merger stages mentioned earlier. In Section V we discuss
our results for the gravitational waves and their spectra. In Section VI we
summarise and conclude. Finally, in Appendix A we show results obtained with
lower resolution as in the main text and study the consistency of our
simulations. Hereafter, we use units where we set $k_{{}_{\rm B}}=G=c=1$
unless otherwise noted.

##  II Matter Model

| $q$ | ${M_{\rm TOV}}$ | ${M_{1}}$ | ${M_{2}}$ | ${R_{1}}$ | ${R_{2}}$ | $\Lambda_{1}$ | $\Lambda_{2}$ | $\tilde{\Lambda}$ | $\Lambda_{1.4}$ | ${f_{\rm mer}}$ | ${f_{3}}$ | ${f^{2,1}}$ | ${f^{2,2}}$ | ${t_{\rm BH}}$
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---
|  | $[M_{\odot}]$ | $[M_{\odot}]$ | $[M_{\odot}]$ | $[{\rm km}]$ | $[{\rm km}]$ |  |  |  |  | ${\rm kHz}$ | ${\rm kHz}$ | ${\rm kHz}$ | ${\rm kHz}$ | $[{\rm ms}]$
Soft_q10 | $1.0$ | $2.02$ | $1.37$ | $1.37$ | $12.37$ | $12.37$ | $537$ | $~{}537$ | $537$ | $475$ | $1.77$ | $4.00$ | $1.44$ | $2.92$ | $9.5$
Soft_q07 | $0.7$ | $2.02$ | $1.64$ | $1.15$ | $12.42$ | $12.24$ | $183$ | $1393$ | $517$ | $475$ | $1.63$ | $3.80$ | $1.55$ | $2.80$ | $5.8$
Soft_q10-NPT | $1.0$ | $2.06$ | $1.37$ | $1.37$ | $12.37$ | $12.37$ | $537$ | $~{}537$ | $537$ | $475$ | $1.76$ | $4.00$ | $1.62$ | $2.85$ | $>37$
Soft_q07-NPT | $1.0$ | $2.06$ | $1.64$ | $1.15$ | $12.42$ | $12.24$ | $183$ | $1393$ | $517$ | $475$ | $1.64$ | $3.85$ | $1.47$ | $2.79$ | $11$
Inter_q10 | $1.0$ | $2.14$ | $1.37$ | $1.37$ | $12.45$ | $12.45$ | $565$ | $~{}565$ | $565$ | $511$ | $1.74$ | $4.00$ | $1.51$ | $2.73$ | $>35$
Inter_q07 | $0.7$ | $2.14$ | $1.64$ | $1.15$ | $12.56$ | $12.30$ | $201$ | $1437$ | $543$ | $511$ | $1.63$ | $3.70$ | $1.42$ | $2.67$ | $>37$
Stiff_q10 | $1.0$ | $2.34$ | $1.37$ | $1.37$ | $12.58$ | $12.58$ | $617$ | $~{}617$ | $617$ | $560$ | $1.74$ | $3.90$ | $1.39$ | $2.49$ | $>37$
Stiff_q07 | $0.7$ | $2.34$ | $1.64$ | $1.15$ | $12.76$ | $12.38$ | $231$ | $1525$ | $591$ | $560$ | $1.59$ | $3.50$ | $1.39$ | $2.48$ | $>38$
Table 1: Properties of cold non-spinning isolated neutron stars used in the
construction of the binary initial data. Specifically, for each model
considered, we list the mass ratio, $q$; the maximum non-rotating mass,
$M_{\rm TOV}$; the gravitational mass of each object at infinite separation,
$M_{<1/2>}$; the proper radius of each object, $R_{<1/2>}$; the tidal
deformability of each object, $\Lambda_{<1/2>}$; the binary tidal
deformability, $\tilde{\Lambda}$; the tidal deformability of a
$1.4~{}M_{\odot}$ NS, $\Lambda_{1.4}$; the instantaneous frequency at the
maximum GW strain amplitude, $f_{\rm mer}$; the characteristic frequencies of
the post-merger phase, $f_{3}\,,f^{2,1}\,,f^{2,2}$; and the collapse time to
black hole formation, $t_{\rm BH}$, for the configurations that collapsed
within the evolution time.

In this section we provide a brief summary of the combined framework used in
our merger simulations. The most important novel ingredient in this framework
is the use of the gauge/gravity duality. The gauge/gravity duality allows one
to translate intractable problems in quantum field theory at strong coupling
and finite density to tractable problems in classical five-dimensional
gravity. This approach has been successfully applied in the context of heavy-
ion collisions [34, 35], condensed matter physics [36, 37] and, more recently,
also astrophysics [38, 39, 40, 41, 42, 43]. The specific gauge/gravity model
which we are using is V-QCD, which we will now briefly review, however, a
detailed introduction to V-QCD and to the EOS at finite temperature and charge
fraction based on V-QCD can be found in [40, 10].

V-QCD includes sectors both for quarks [44, 45] and for gluons [46, 47]. It is
based on a string theory setup in the Veneziano limit [48, 49] in which both
the number of colors $N_{c}$ and flavors $N_{f}$ are large, but their ratio is
$\mathcal{O}(1)$ [13, 50, 51] such that the quark degrees of freedom strongly
back-react to the gluon dynamics. However, the model cannot be strictly
derived from string theory, and one switches to an effective “bottom-up”
approach where the action of the model is tuned to agree with QCD data at
finite $N_{c}$ and $N_{f}$. Therefore the model includes a relatively large
number of parameters that are tuned to reproduce known features of QCD such as
asymptotic freedom, confinement, linear glueball and meson trajectories and
chiral symmetry breaking. The remaining freedom is constrained by lattice data
for large-$N_{c}$ pure Yang-Mills theory [52] and data for $N_{c}=3$ QCD with
$N_{f}=2+1$ flavors at physical quark masses [53, 54] at small baryon-number
density. Such an effective approach turns out to be useful to model dense
nuclear and quark matter in the region where the QCD-coupling is large and
traditional nuclear theory calculations and perturbation theory are not
reliable [54, 55, 56].

In the current version of the model baryons are implemented in a homogeneous
approximation [14]. This approximation is natural at densities well above
$n_{s}$, where the distance between neighbouring nucleons becomes smaller than
their diameter and a description in terms of homogeneous matter is expected to
work well. At densities below $n_{s}$ this homogeneous approximation breaks
down and traditional nuclear matter models become more reliable. In the
combined framework this is taken into account by constructing hybrid EOSs
whose cold low-density part is modelled by nuclear theory and the dense
baryonic and quark matter part via the gauge/gravity duality. For the cold
nuclear theory part we use a combination 111By using the combination of the
HS(DD2) and APR models together with V-QCD, we obtain EoSs that agree well
with all observational constraints. Using only the HS(DD2) model would lead to
a slight tension with the measurement of the tidal deformability from
GW170817. of the Hempel-Schaffner-Bielich (HS) EOS [17] with DD2 relativistic
mean field theory interactions [16] and the Akmal-Pandharipande-Ravenhall
(APR) model [58].

The prediction of V-QCD for the EOS in the nuclear matter phase suffers from a
generic limitation of gauge/gravity duality: the temperature dependence is
trivial. In [10] this issue was solved by using a simple van der Waals model,
i.e., gas of nucleons, electrons and mesons with excluded volume correction
for the nucleons and an effective potential. The van der Waals model was
chosen to match with the V-QCD prediction for the cold EOS of nuclear matter
and then used to extrapolate the result to finite temperature. The model was
further improved to allow deviation from beta-equilibrium by imposing the
charge fraction dependence of the HS(DD2) model [17, 16] for nuclear matter
and a simple model of free electrons for quark matter. In the regime of low
density, both the dependence on the temperature and the charge fraction are
following the HS(DD2) model. Finally, our framework contains a strong first
order phase transition between the nuclear and quark matter phases at low
temperatures which ends on a critical point at finite temperature and density.
Consequently, there is a mixed nuclear-quark matter phase which was obtained
in [10] by carrying out a Gibbs construction depending on all the three
variables, i.e., the density, the temperature, and the charge fraction.
Surface tension of the nuclear to quark matter interface is neglected. This is
expected to be a good approximation at the temperatures that we encounter in
our simulations, because the latent heat at the transition is very high,
comparable to the free energy density of the nuclear matter phase.

We here use three variants of the EOS, which were established in [10] in order
to represent the parameter dependence of the model. These are called the soft
(Soft), intermediate (Inter), and stiff (Stiff) variants; and reflect three
different choices of the parameters in the action of V-QCD [54]. In Fig. 1
(left) we show cold beta-equilibrium slices of these three EOSs with
uncertainty bands from nuclear theory [18] (blue) and perturbative QCD [19]
(orange). Red, green and blue lines are the Soft, Inter and Stiff versions of
V-QCD where the dotted part of these curves mark the first order phase
transition between the baryonic (solid) the quark matter (dashed) phase. In
addition we show markers for the central densities reached in various isolated
non-rotating stars that we use to initialise the binary systems in our
simulations. Light blue and pink regions mark the residual freedom of
polytropic parametrizations of the EOS [59, 60] and of the V-QCD model [29],
respectively, after imposing constraints from neutron star observations. In
Fig. 1 (right) we show the corresponding mass-radius relations of non-rotating
stars together with various error bands of the direct mass measurement (green
area) of the heaviest known pulsar PSR J0740+6620 [20, 21] ($M=2.08\pm
0.07M_{\odot}$) and direct radius measurements of PSR J0740+6620 [25, 26] and
PSR J0030+0451 [22, 23] obtained by the NICER experiment (pink ellipses) as
well as from the measurement of the X-ray binary 4U 1702-429 [24] (cyan area;
see [31] for a more detailed analysis of the impact of the NICER results).
Finally, shown with a black dashed line is a theoretical lower bound for the
radii as computed from considerations on the threshold mass by Koeppel et al.
[27]. We note that we have checked that our EOSs respect the upper bound on
the binary tidal deformability $\tilde{\Lambda}<720$ (low-spin priors)
obtained from GW171817 and are well within the bounds recently proposed on
general parameterization of the sound speed in neutron stars [61] (see Table
1).

##  III Initial Data and Evolution Setup

To generate the binary neutron-star initial data we use the recently developed
Frankfurt University/Kadath (FUKA) [62] solver library. FUKA is based on an
extended version of the KADATH spectral solver library [63], which has been
specifically designed for numerical relativity applications. To generate the
initial data for compact binaries, FUKA uses the eXtended Conformal Thin
Sandwich formulation of Einstein’s field equations [64, 65, 62]. Initially
FUKA solves the initial data using force-balance equations to obtain a binary
in the quasi-circular orbit approximation. In an effort to minimize the
residual eccentricity during the inspiral, an eccentricity reduction step is
performed using estimates for the orbital- and radial-infall- velocities at
$3.5$-th post-Newtonian order. As shown in [62], this procedure leads to a
significant reduction of the eccentricity in asymmetric and spinning compact
object binaries. We find the advanced handling of eccentricities in FUKA
essential to obtain accurate initial data for the unequal mass binaries
studied in this work.

Figure 2: Correlation between the maximal temperature, maximal density and
quark fraction in the three regions of interest, for all three equations of
state and both mass ratios. We mark the time slices used for Fig. 3 with
vertical black lines. The abbreviations HQ, WQ and CQ stand for the hot-quark,
warm-quark and the cold-quark stages respectively.

For the binary evolution we make use of the Einstein Toolkit [66]
infrastructure that includes the fixed-mesh box-in-box refinement framework
Carpet [67]. For the simulations presented in the main text we use six
refinement levels with finest grid-spacing of $\Delta_{\texttt{H}}:=221\,{\rm
m}$. To check the consistency of our results we have performed a series of
simulations for the Soft EOS with low ($\Delta_{\texttt{L}}:=369\,{\rm m}$)
and medium ($\Delta_{\texttt{M}}:=295\,{\rm m}$) resolution which is discussed
in Appendix A. For the spacetime evolution we use Antelope [68], which solves
a constraint damping formulation of the Z4 system [69, 70]. Furthermore, to
evolve the hydrodynamic part we use the Frankfurt/Illinois (FIL) general-
relativistic magnetohydrodynamic code [68], which is based on the
IllinoisGRMHD code [71]. FIL implements fourth order conservative finite-
differencing methods [72], enabling a precise hydrodynamic evolution even at
low resolution. Furthermore, FIL is able to handle tabulated EOSs that are
dependent on temperature and electron-fraction. In principle, FIL also
includes a neutrino leakage scheme that can handle neutrino cooling and weak
interactions. However, in its current form V-QCD does not include a
description for neutrinos which is why we do not make use of the leakage
scheme in this work.

Finally, in Table 1 we summarise the properties of cold non-spinning isolated
neutron stars used in the construction of the binary initial data as well as
various characteristic frequencies extracted from the post-merger waveforms
presented in Sec. V.

Figure 3:  From top to bottom rows we show snapshots of the baryon-number
density, temperature and quark fraction in the orbital plane at three
different times $t-t_{\rm mer}=1.4,\,3.3,\,5.6\,\,{\rm ms}$ (from left to
right) that are representative for HQ, WQ and CQ stages in the Soft_q07
simulation. In addition we mark contours for the quark fraction $Y_{\rm
quark}=0.01,\,0.1,\,0.5$ by dotted, dashed and solid lines, respectively.

##  IV Merger Dynamics and Quark Formation

Figure 4:  Left: Phase diagram of the Soft_q07 post-merger evolution, with
regions of confined, mixed and deconfined phases of matter shown. Points with
a specific colour all occur within time intervals indicated by darker colour
bars. Right: The amount of quarks in the orbital plane, measured using the
same times as in Fig. 3. We normalize the measured amount of quarks by the
maximal amount found in all three time slices presented here.

In this section we present the merger dynamics and discuss the three
mechanisms in which quarks are formed during the post-merger evolution. We
have performed numerical simulations of BNS systems using mass ratios
$q=1\,,0.7$ where the neutron stars have an irrotational fluid profile. For
each configuration, the initial data and evolution have been computed using
the Soft, Inter, and Stiff EOSs in order to compare the influence of stiffness
on the quark production and lifetime of the remnant HMNS whose lifetime is
expected to be on the order of one second [73, 74]. A comparative overview of
the results from these numerical simulations can be seen in Fig. 2 where we
show the evolution of the normalized maximum density, $(n_{\rm b}/n_{\rm
s})_{\rm max}$, the maximum temperature, $T_{\rm max}$, and the quark
fraction, $\langle Y_{\rm quark}\rangle$. To further characterise the
distribution of quark matter within the remnant HMNS, we define three regions
that follow the centre-of-mass of the HMNS in the orbital plane. We will refer
to these regions from now on as the inner core ($R<3\,{\rm km}$), outer core
($3\,{\rm km}$ $<R<$ $6\,{\rm km}$) and inner crust ($6\,{\rm km}$ $<R<$
$9\,{\rm km}$). An illustration along with further details regarding the
motivation behind these choices is discussed in Appendix C. The quark fraction
$Y_{\rm quark}$ is currently computed in the post-analysis of the evolutions.
In order to obtain the associated grid function values we use a multi-linear
interpolation scheme for the given tabulated EOS of a respective model, in
which discrete data points $Y_{\rm quark}=Y_{\rm quark}(n_{\rm b},T,Y_{\rm
e})$ are stored, as well as the 2D slices of the $n_{\rm b}$, $T$ and $Y_{\rm
e}$.

In all three models, comparing $\langle Y_{\rm quark}\rangle$ to the maximal
temperature shows a clear correlation between the appearance of quark matter
and the increase in temperature due to shock-heating during the initial merger
stage which is coloured red in Fig. 2. For the remainder of this work, we will
reference this as hot quark (HQ) production. The amount of HQs produced is
significantly larger for the Soft EOS than for the stiffer models which is not
unexpected as the densities inside a neutron star with fixed mass is larger
when matter is soft and is, therefore, easier to produce strong shocks during
the collision. Furthermore, the amount of HQs produced is significantly larger
in the unequal mass case than in the equal mass case because the initial
density of the heavier companion is higher before merger in unequal binaries
which subsequently results in a more violent collision and stronger shocks.
Another significant difference between equal and unequal mass binaries is that
in the latter HQs are formed not only in the inner core (and a tiny amount in
the outer core), but also in the outer core and inner crust region. This has
to do with the fact that the distribution of matter in the unequal mass case
is highly non-axisymmetric which leads to HQs in the hot off-central region of
the HMNS.

In the case of the Soft EOS, we see during the cooling of the remnant and as
the core becomes more dense an additional quark production channel occurs. We
will designate this production mechanism as warm quarks (WQ). The formation of
WQs appears to be a result of the complicated post-merger dynamics that leads
to a periodic expansion and contraction of the HMNS. During the contraction
the density increases close to the rotation axis and leads to the formation of
quarks. In the subsequent expansion the density decreases and quarks
transition back to baryonic matter. In the Soft model we find that WQs are
produced in the equal mass case $\approx 2\,{\rm ms}$ after the end of the HQ
production, while for unequal masses we find a continuous transition between
the production of HQs and WQs. For the stiffer models we do not find WQs,
however, we expect their formation to be possible with a heavier binary which
will be studied thoroughly in a future work. To capture accurately the
delicate combination of temperature and density required for WQs, the use of
high resolution numerical methods is essential. This is demonstrated in
Appendix A where we compare our results with two lower resolutions for the
Soft model and find that WQs can be absent in our intermediate and low
resolution simulations.

Finally, when the HMNS becomes approximately circular and its compactness
increases, a cold but dense quark matter core can be formed which we denote as
cold quark (CQ) production. We find CQ production to be present only for the
Soft model, which ultimately undergoes a phase-transition-triggered collapse
for the configurations considered. Notice that due to the large latent heat in
the V-QCD model, no stable quark matter cores can be formed and therefore
formation of a quark core inevitably leads to the collapse of the remnant. The
starting time of CQ production depends on the mass ratio and the total mass of
the system. For the unequal-mass binary it starts already $\approx 4\,{\rm
ms}$ after merger-time, while in the equal-mass case the starting time is
approximately two milliseconds later ($\approx 6\,{\rm ms}$). A close
inspection of the maximum number density in the core of the unequal-mass
binary reveals that CQs set in during the fifth bounce of the merger remnant.
In this stage the maximal temperature in the core is lower than in the
exterior, in particular for the equal-mass binary. We find this stage to last
for $\approx 4\,{\rm ms}$ ($\approx 2\,{\rm ms}$) for the equal-(unequal-)mass
binary, before the density and temperature in the core eventually start to
increase rapidly leading to the formation of a black-hole horizon.

To further explore the appearance of HQ, WQ and CQ we show in Fig. 3 snapshots
of the number density (top), temperature (middle) and quark fraction (bottom)
in the orbital plane at three representative subsequent times (left to right)
of the post-merger evolution of the Soft_q07 configuration. The times have
been chosen based on the times of significant HQ, WQ, and CQ production,
respectively, as denoted by vertical dotted lines in Fig. 2. Additionally in
Fig. 3 we indicate by dotted, dashed and solid lines the outer contours of
regions that contain quark fractions larger than $0.01$, $0.1$ and $0.5$,
respectively. The plots for $t-t_{\rm mer}=1.4\,{\rm ms}$ clearly show the
presence of HQs in the hottest regions outside the dense core of the HMNS.
Typical for this early post merger state is the formation of multiple
disconnected hot pockets with temperatures well above $50$ MeV inside the HMNS
which, in this case, result in two disconnected regions where HQs are formed.
The second column shows snapshots at $t-t_{\rm mer}=3.3\,{\rm ms}$
corresponding to the formation of a single patch of WQs. We note that a
precise definition of the WQ stage is difficult since it sensitively depends
on the combination of densities and temperatures outside the hottest and
densest regions of the star such that a transition to a mixed-phase is
possible. These conditions are typically realized in regions outside, but
close to the center of the HMNS where both temperature and density are
significantly below their maximal values. Finally, in the third column we show
snapshots at $t-t_{\rm mer}=5.6\,{\rm ms}$ where a pure quark matter core
($Y_{\rm quark}=1$) has already been formed and ultimately leads to a phase-
transition-triggered collapse of the HMNS. The maximal density of the quark
matter core is well above $9\,n_{s}$ (and rising), but the temperature is very
low further motivating our classification of quark matter in this region as
CQs. Also noteworthy is the appearance of a small disconnected portion of WQs
that forms in a small region outside the dense and cold centre. It is
interesting to compare our cold quark core to [75], who find the quark core to
be hot. The reason for this difference is that [75] models the temperature
dependence by adding a gamma-law to the cold EOS model where the temperature
simply scales with the density and does not take into account the change in
composition as in V-QCD.

To study the composition of matter in terms of the maxima of $\langle Y_{\rm
quark}\rangle$ compared to $T_{\rm max}$ and $n_{b}^{\rm max}$, we show in
Fig. 4 (left) a phase diagram of the Soft V-QCD model in the temperature and
baryon-number density plane where red, green and blue regions indicate
confined, mixed and deconfined phases, respectively; and solid lines mark the
phase boundaries in beta equilibrium which end at a critical point marked by a
black cross. Additionally, we show as a reference dashed and dotted lines
indicating the phase boundaries at approximately the minimal ($Y_{e}=0.05$)
and maximal ($Y_{e}=0.09$) values of the charge fraction inside the HMNS
during the evolution. These values differ from that of the beta equilibrium,
because we determine $Y_{e}$ by demanding beta equilibrium inside the stars
that are used to initialize the binary system, but advect $Y_{e}$ during the
subsequent evolution thereby driving the fluid out of beta equilibrium.
Finally, colored symbols mark representative results of the $q=0.7$ simulation
for the maximal temperature (squares), baryon-number density (stars) and quark
fraction (circles) collected during the early HQ stage ($t-t_{\rm
mer}=0-2\,{\rm ms}$, red), the WQ stage ($t-t_{\rm mer}=2-4\,{\rm ms}$, green)
and the CQ stage ($t-t_{\rm mer}=4-6\,{\rm ms}$, blue). From the red symbols
one clearly sees the maxima of the quark fraction appearing close to the
maxima in temperature at approximately $60$ MeV, while the maxima in the
number density appear at much lower temperature at the bottom of the phase
diagram, illustrating our motivation to classifying these quarks as HQs. The
green symbols represent examples for the WQ stage where $Y_{\rm quark}$ is
maximal at intermediate temperature ($\approx 20$ MeV) and where the maximal
density ($\approx 4\,n_{s}$) occurs. The maximal temperature is significantly
higher ($\approx 60$ MeV) in this case, but appears at much lower densities
($\approx 1\,n_{s}$). Furthermore, the blue symbols mark the CQ production
where the maximal quark fraction coincides with the maximum number density at
very low temperature.

Finally, in Fig. 4 (on the right) we show our measurements of the amount of
quark matter during the three separate time-slices used in Fig. 3. For this we
define $126\times 35$ bins with equal size in linear space, in the
$[n_{s},\,10\,n_{s}]\times[0,140\,\rm{MeV}]$ region of the number density-
temperature $(n_{b},T)$ plane, which captures the whole domain of quark
production. To do the counting we evaluate the quantity $N_{\rm
quark}=\sum_{j}V_{j}Y_{{\rm quark},j}n_{b,j}$ separately for each bin, where
the sum is performed over all the grid-cells $j$ of volume $V_{j}$ on a fixed
time-slice in the orbital plane, and then normalized by the maximal count on
the whole grid - where the maximal $N_{\rm quark}$ is marked with a star. The
first row shows HQs at $t-t_{\rm mer}=1.4\,{\rm ms}$, which span temperatures
$\approx 40-60$ MeV and densities $\approx 1-2\,n_{s}$. The second row shows
WQs at $t-t_{\rm mer}=3.3\,{\rm ms}$ produced at temperatures $\approx 10-30$
MeV and densities $\approx 2-3\,n_{s}$. Finally, the third row for $t-t_{\rm
mer}=5.6\,{\rm ms}$ clearly shows that significant amounts of WQs and CQs can
be present simultaneously. However, from the the colour code of this plot one
can also see that the majority of quarks (indicated in yellow) are cold and
dense.

Figure 5:  Shown above is the $h^{22}_{+}$ of the strain, the power-spectral
density of the effective strain, and the spectrogram of the $h^{22}_{+}$ of
the Soft, Inter, and the Stiff EOSs considered in this work for $q=1\,,0.7$
and $\mathcal{M}_{\rm chirp}=1.188$. Also included are the related spectra for
Soft-NPT which does not include a phase transition to quark matter. In the
right two panels the dashed blue and orange lines correspond to the $f^{2,2}$
and $f^{2,1}$ peaks respectively, the dashed yellow lines corresponds to the
$f_{3}$ peak as measured in the spectrogram, and the star denotes the peak
merger frequency, $f_{\rm mer}$. In the right figure, the white line traces
the maximum in the spectrogram. Finally, the sensitivity curves in the middle
plot are related the current sensitivity of advanced LIGO and the Einstein
Telescope respectively [76, 77].

##  V Gravitational Wave Analysis

In this section we analyse the gravitational waveforms and their spectral
properties. We use the Newman-Penrose formalism [78, 79] to relate the Weyl
curvature scalar $\psi_{4}$ to the second time derivative of the polarization
amplitudes of the gravitational wave strain $h_{+,\times}$ via

$\ddot{h}_{+}+i\ddot{h}_{\times}=\psi_{4}:=\sum_{\ell=2}^{\infty}\sum_{m=-\ell}^{m=\ell}\psi_{4}^{\ell,m}{{}_{-2}Y}_{\ell,m}\,,$
(1)

where ${}_{s}Y_{\ell,m}(\theta,\phi)$ are spin-weighted spherical harmonics of
weight $s=-2$. From our simulations we extract the modes $\psi_{4}^{\ell,m}$
with a sampling rate of $\approx 26\,{\rm kHz}$ from a spherical surface with
$\approx 440\,{\rm km}$ radius centred at the origin of our computational
domain and extrapolate the result to the estimated luminosity distance of
$40\,{\rm Mpc}$ of the GW170817 event [28]. In addition, we fix the angular
dependence of the spherical harmonics part by the viewing angle
$\theta=15^{\circ}$ determined from the jet of GW170817 [80] and set
$\phi=0^{\circ}$ without loss of generality. We restrict our analysis to the
dominant $\ell=m=2$ and sub-dominant $\ell=2,\,m=1$ modes of the expansion
(1). To analyze the spectral features of the waveforms we follow [81] and
compute the power spectral density (PSD) defined as

$\displaystyle\tilde{h}^{\ell,m}(f):=\frac{1}{\sqrt{2}}\biggl{(}$
$\displaystyle\left|\\!\int\\!\mathrm{d}t\,{\rm e}^{-2\pi
ift}h^{\ell,m}_{+}(t)\right|^{2}+$
$\displaystyle\left|\\!\int\\!\mathrm{d}t\,{\rm e}^{-2\pi
ift}h^{\ell,m}_{\times}(t)\right|^{2}\biggr{)}^{1/2}\,,$ (2)

where the time integration is performed over the interval $t-t_{\rm
mer}=0-10\,{\rm ms}$ or up to the collapse to a black hole as is the case for
the Soft binaries. Furthermore, we use spectrograms to study the time
dependence of the gravitational wave frequency distribution. To perform the
Fourier transform for the spectrograms we use time-windows of $3\,{\rm ms}$
centred at every $\approx 0.04\,{\rm ms}$ of our waveform data. We also study
the phase difference between the $h_{+}^{2,2}$ and $h_{\times}^{2,2}$ modes by
computing the instantaneous gravitational wave frequency $f_{\rm GW}$ defined
as

$f_{\rm GW}:=\frac{1}{2\pi}\frac{d\phi}{dt}\,,\qquad\phi:={\rm
arctan}\left(\frac{h_{\times}^{2,2}}{h_{+}^{2,2}}\right)\,.$ (3)

Finally, we can measure from $f_{\rm GW}(t)$ the instantaneous frequency at
merger time defined as $f_{\rm mer}=f_{\rm GW}(t=0)$.

In Fig. 5 we include a summary of our gravitational wave analysis of the
configurations described in Tab. 1. From top to bottom we show results for the
different EOS models and the two mass ratios we have simulated. In the first
column we show the dominant $h_{+}^{2,2}$ gravitational wave strain component
with the time periods where the dominant quark production mechanism, HQ, WQ
and CQ, are highlighted in red, green and blue, respectively.

The most prominent signature of quark matter in the waveform is clearly the
earlier termination of the signal due to a phase-transition-triggered collapse
to a black hole. This happens in our simulations only for the Soft model,
which for both mass ratios forms a CQ stage that leads to a collapse of the
HMNS. We note that models that include hyperons [82, 83] can lead to
signatures similar as those from a first order deconfinement phase transition
[3]. In both cases a softening of the EOS leads to more compact merger
remnants that is closer to collapse. We estimate the lifetime $t_{\rm BH}$ of
our merger remnants by the time at which the strain amplitude starts to decay
exponentially as characteristic for the black-hole ring-down. In cases where
no horizon is formed during simulation time, we provide lower limits on
$t_{\rm BH}$ by the final time of the simulation. The values for $t_{\rm BH}$
we obtain in this way are listed in Table 1. Comparing the waveforms of the
Soft model for $q=0.7$ (second row) and $q=1$ (sixth row) shows that the
unequal-mass case collapses $\approx 4\,{\rm ms}$ earlier than the equal mass
case. This behaviour is consistent with the threshold mass analysis [84, 85,
86], which finds that highly asymmetric mergers tend to collapse at lower
total mass than symmetric ones.

In the second column of Fig. 5 we show the PSDs computed via Eq. (V). Solid
lines represent the dominant $\ell=m=2$ mode, while light dotted lines are the
sub-dominant $\ell=2,\,m=1$ mode. The dotted black line is the contribution to
the PSD that comes from the inspiral part ($t-t_{\rm mer}<0$), which is not
included in the other curves. In addition we show estimated sensitivities of
the advanced LIGO [76] (purple dot-dashed) and the Einstein Telescope [77]
(cyan dot-dashed line). The PSD shows the typical three-peak structure from
which a number of characteristic post merger frequencies can be extracted. The
most prominent peak $f^{2,2}$ (dashed blue) corresponds to the dominant post-
merger frequency of the $\ell=m=2$ mode of the signal. Comparing the three EOS
models for $q=1,0.7$ we find good agreement amongst their values of $f^{2,2}$
and only a small shift to lower frequencies with increasing stiffness of the
EOS. Most notably, the $f^{2,2}$ peaks for the Soft_q10 and Soft_q10-NPT are
equivalent though the peak-width for Soft_q10-NPT is slightly broader (where
the NPT simulations do not include a transition to quark matter as we will
discuss below). In addition, we mark the maximal frequency $f^{2,1}$ of the
sub-dominant $\ell=2,m=1$ mode by a dashed orange line and the merger
frequency $f_{\rm mer}$ by a red star. Similar to the dominant mode, the
$f^{2,1}$ frequencies (dashed orange) decrease only slightly with the increase
in stiffness, but there is an appreciable difference when measuring the
$f^{2,1}$ for the Soft EOS with and without a phase transition to quark matter
in the equal-mass binaries. Overall we find that the sub-dominant mode in the
equal mass case is by at least one order of magnitude suppressed with respect
to the unequal-mass case. This is expected, since unequal-mass binaries result
in less axially symmetric remnants than equal-mass binaries and, therefore,
produce a richer post-merger frequency spectrum. In this sense gravitational
waves from asymmetric binaries potentially provide more information about the
EOS than symmetric binaries. However, in both cases the amplitude of the sub-
dominant mode is orders of magnitude too low to be seen by current detectors,
even when assuming the advantageous luminosity distance and sky localization
of GW170817.

In the third column of Fig. 5 we show spectrograms where the colours from dark
blue to light yellow indicate the distribution of low and high spectral weight
during the evolution. In addition, we show lines for various characteristic
frequencies: dashed blue and orange lines mark the$f^{2,2}$ and $f^{2,1}$
respectively - as extracted from the PSD. The dashed yellow line marks $f_{3}$
which we extract directly from the spectrogram because it is difficult to
obtain unambiguously from the PSD. Finally, the solid red line is $f_{\rm GW}$
as computed from (3) and the solid white line follows the maximum values of
the spectrogram. Note that the values of $f^{2,2}$ and $f^{2,1}$ extracted
from the PSD align exceptionally well with the $f_{1}$ and $f_{2}$ peaks seen
in the spectrogram; however, the $\ell,m=(2,2),(2,1)$ PSDs do not provide
identifying information related to the $f_{3}$ peak which is measured directly
from the spectrogram.

###  V.1 Collapse Behaviour and Lifetime

To isolate the impact of quark matter on the lifetime we also perform
simulations where we exclude the thermodynamically preferred quark matter
phase at high temperature and densities in the Soft V-QCD model (see first and
fifth row in Fig. 5) which we designate as the Soft-NPT model. Unlike the case
with quark matter, the equal-mass case without quark matter does not collapse
during our simulation time, however, we find that in the unequal-mass case, a
binary using Soft-NPT model collapses after $\approx 11\,{\rm ms}$, i.e. ,
roughly $5\,{\rm ms}$ later than with the Soft model. The early collapse of
our simulations with the Soft V-QCD model point to a strong tension with the
expected second-long lifetime of GW170817 [73, 74]. This means our results
disfavour the Soft EOS and, therefore, allow us to further constrain the V-QCD
model with information gained from the post-merger phase. This is an important
result since, as discussed in Sec. II, the cold part of the Soft V-QCD model
is fully consistent with currently known observations and the constraints from
the post-merger lifetime can only be inferred from merger simulations.

Furthermore, obtaining faithful estimates for the lifetime of merger remnants
that survive significantly longer than $\approx 10\,{\rm ms}$ is difficult in
our setup. There are mainly two reasons for this: First of all, resolving the
formation of quark matter requires relatively fine resolution, which makes
long simulation times numerically extremely expensive. We therefore restrict
our high-resolution simulations to $\approx 10\,{\rm ms}$ after merger.
However, we present in Appendix B simulation results with medium resolution up
to $>35\,{\rm ms}$ for the Inter and Stiff model. This allows us to
approximate the lower bound of their lifetime to be $>35\,{\rm ms}$. Second of
all, the competition between turbulent heating in the presence of magnetic
fields and cooling effects by neutrino emission can become important on time-
scales on the order of several $100\,{\rm ms}$ [87]. Since we do not include
either of these effects in our simulations we expect the thermodynamic
properties, in particular the temperature distribution inside the star, not to
be reliable on such long time-scales. However, we expect these thermal effects
to be negligible on the comparably short time-scales the life-time estimates
listed in Table 1 are based on.

##  VI Summary and Conclusion

We have studied fully 3+1 dimensional GRHD simulations of GW170817-like events
using a unified model for QCD that includes the deconfinement phase transition
from baryonic to quark matter. The focus of our work was the formation of
quark matter during and after the merger of the two neutron stars. For this we
studied three variants of the unified model, namely a soft, an intermediate
and a stiff version. We identified three distinct mechanisms in the post-
merger evolution during which mixed baryonic and quark matter and pure quark
matter is formed and which we denote as: i) hot quark (HQ), ii) warm quark
(WQ) and iii) cold quark (CQ) production. We find that some or all of these
stages may be present in the evolution depending on the masses of the stars
and the EOS model employed.

Specifically, the HQ production appears at merger time and we find it to last
for approximately two milliseconds, except for the stiff model where only a
small amount of quarks are produced at merger time. During this time, strong
shocks that are formed during the initial contact of the two stars lead to a
steep rise in temperature. This causes a local formation of quark matter in
the hottest, but not in the dense regions of the merger remnant. As the
remnant begins to cool and lose its differential rotation causing an increase
in density, the production of WQs may be possible depending on the model and
binary configuration. The WQ production is the result of complicated merger
dynamics and is observed in regions that are typically neither the hottest nor
the densest regions of the HMNS. In this work WQs are only observed in the
soft model. It is important to note that the existence and starting time of
the WQ stage depends strongly on the details of the complicated post-merger
evolution during which the HMNS pulsates in a highly non-axisymmetric way and,
as such, turns out to be particularly sensitive to the mass ratio of the
binary system. Additionally, we find that WQs are produced earlier and also in
a larger quantities in non-equal mass binaries than in equal-mass binaries.
Finally, the CQ production appears after the most violent post-merger period
has settled and the density in the centre increases and cools down. The
production of CQ is then formed in the dense and cold centre of the merger
remnant. In this work, the CQ stage is followed by a phase-transition-
triggered collapse to a black hole, using the terminology of [8]. The reason
for this is the large latent heat of the deconfinement phase transition in the
V-QCD model at low temperatures which does not allow stable cold quark matter
cores inside isolated stars or HMNSs.

Finally, we studied the waveforms and their spectral properties and find the
most significant imprint is the early termination of the gravitational wave
signal due to the phase-transition-triggered collapse due to CQ production. By
this we find the short post-merger lifetime of $\lesssim 10\,{\rm ms}$
obtained from the soft version of V-QCD to be in tension with the expected
second long lifetime of GW170817. It is remarkable that, even though the
majority of the post-merger period ($\sim 6ms$) of the unequal-mass case for
the soft model includes a phase transition to a mixed or pure quark phase,
there is essentially no measurable imprint on the amplitude or frequency of
the gravitational wave strain. Moreover, this tension also means that there
are new constraints for the predictions for the QCD equation of state arising
from the merger simulations. As the soft version of the unified model is
basically excluded by the lifetime estimate of the remnant, the soft versions
should be also excluded from the analysis of the cold [29, 31] and hot [10]
equation of state. We do not attempt to make such constraints precise in this
article, but as an example we discuss the location of the critical point of
the nuclear to quark matter transition. It was found in [10] that the
temperature $T_{c}$ of the critical point varies roughly between 110 MeV and
130 MeV, with the soft version predicting the highest temperature. Exclusion
of the soft model therefore means that $T_{c}\lesssim 120$ MeV within the
framework.

There are multiple directions we plan to explore in future work. In the
present work we only considered non-spinning binaries with individual masses
consistent with GW170817, which prefers the low-spin prior. It will be
important to investigate the dependence of various quark matter production
stages also on different individual masses and spins of the binary components.
Another direction worth studying is the impact of the deconfinement phase
transition on the threshold mass, which has been analysed extensively in [84,
85, 86, 88] for purely hadronic EOSs as well as some hybrid models that
include first order phase transitions[85, 86, 88]. Notice also that our
simulations did not take into account the surface tension at the interface
between the nuclear and quark matter. The value of the surface tension has not
yet been computed in the V-QCD framework, and even if the value was known, it
is not straightforward to take it into account: strong first order phase
transitions in the presence of a sizeable surface tension proceed through
bubble nucleation, which is a dynamical process that cannot be easily taken
into account at the level of the equation of state. It is however clear that
the impact of the surface tension would be to somewhat reduce the amount of
mixed phase and suppress the quark production. We plan to study in a future
work how large this effect is. Furthermore, it would be interesting to include
cooling effects from neutrino emission which could have non-trivial impact on
the quark formation, especially during the HQ and WQ stages. Finally, because
of the high computational costs we restricted our high-resolution (medium-
resolution) simulations to $\approx 10\,{\rm ms}$ ($\approx 40\,{\rm ms}$)
after the merger. This allowed us to obtain reliable predictions for the quark
matter production and the lifetime from the simulations that collapse during
this short post-merger period. However, we can only speculate about the
ultimate fate of the cases that do not collapse during our simulation time. We
plan to extend these simulation to the order of $1$ second by imposing axial
symmetry in future long-term evolutions to obtain improved estimates for the
lifetime of our long-lived HMNSs.

###### Acknowledgements.

CE acknowledges support by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) through the CRC-TR 211 ’Strong-interaction matter under
extreme conditions’– project number 315477589 – TRR 211. ST, KT, and LR
acknowledge the support by the State of Hesse within the Research Cluster
ELEMENTS (Project ID 500/10.006). TD and MJ have been supported by an
appointment to the JRG Program at the APCTP through the Science and Technology
Promotion Fund and Lottery Fund of the Korean Government. TD and MJ have also
been supported by the Korean Local Governments – Gyeongsangbuk-do Province and
Pohang City – and by the National Research Foundation of Korea (NRF) funded by
the Korean government (MSIT) (grant number 2021R1A2C1010834). LR acknowledges
funding by the ERC Advanced Grant “JETSET: Launching, propagation and emission
of relativistic jets from binary mergers and across mass scales” (Grant No.
884631). The simulations were performed on HPE Apollo HAWK at the High
Performance Computing Center Stuttgart (HLRS) under the grant BNSMIC. For our
visualisations we made an extensive use of the Kuibit library [89].

Figure 6: Shown are the consistency results of the Soft_q10 dataset. Left:
consistency as measured in the gravitational wave strain $h^{22}_{+}$
polarisation amplitude. Although slight changes in the amplitude and phase are
observed, the collapse time is the same in all resolutions. Right: the
corresponding quark fraction for the three resolutions considered. Most
notable is the appearance of WQ in the high-resolution run further emphasising
the need for high resolutions to capture fine-structure details.

##  Appendix A Convergence Analysis

To ensure a high level of confidence in the discussed results, we have
performed a consistency study with the Soft_q10 initial data to ensure that
the overall dynamics of the merger and the phase transition-induced collapse
are not simply artefacts of the numerical accuracy. To this end, we have
employed three resolutions which are characterised by the finest grid spacing
of $\Delta_{\texttt{L}}:=369\,{\rm m}$, $\Delta_{\texttt{M}}:=295\,{\rm m}$
and $\Delta_{\texttt{H}}:=221\,{\rm m}$ where $\Delta_{\texttt{H}}$ is the
resolution used for the results presented in the main body of the text.
Additionally, we employ a Courant static time step such that $dt=C_{dt}{\rm
max}(dx)$ where $C_{dt}=0.2$ for our simulations. Therefore, an increase in
spatial resolution results in an increase in temporal resolution.

As shown in Fig. 6 (left) the gravitational waveforms show only a slight
variation in the relative phase during the inspiral and a systematic decrease
in the post-merger amplitude, however, all resolutions considered result in
the same collapse time. Furthermore, we see in Fig. 6 (right) that the
existence and identification of the different stages of quark production is
robust under the change of resolution. Notably, the intermediate production of
WQ appears for $\Delta_{\texttt{H}}$ which is not surprisingly absent in the
$\Delta_{\texttt{M}}$ resolution given the sensitivity required in temperature
and density for WQ to appear especially for such a short duration ($\sim
0.3\,{\rm ms}$). Remarkably, the peaks in $Y_{\rm quark}$ follow a similar
spacing in time, where individual peaks tend to form more often with an
increase in resolution.

##  Appendix B Medium Resolution Runs

Given the consistent nature of our evolutions as discussed in A, we conducted
additional evolutions at $\Delta_{\texttt{M}}$ resolution up to $\approx
35\,{\rm ms}$ for the Inter and Stiff models. The results of these runs are
summarised in Fig. 7 where the behaviour of the thermodynamic quantities is
shown to be consistent with the results in Fig. 2 as well as the waveforms
shown in Fig. 5. This allows us to place lower bounds on the collapse time,
namely $\approx 35\,{\rm ms}$. Although these results provide evidence that
the Inter and Stiff models are consistent with the expected second long
lifetime of GW170817, we want to emphasize that our simulations do not include
contributions of magnetic fields and neutrino cooling which could be relevant
for the long-term stability of a HMNS.

Figure 7: The results of the Stiff and Inter medium resolution runs for both
investigated mass ratios. In neither do we observe a collapse throughout the
whole $\sim 35$ ms of the post-merger evolution.  Figure 8: A slice on $(x,y)$
plane from the Inter_q10-MR run, plotted at the latest available time. Shown
are the temperature (left part) and the baryon-number density (right part),
while the black dashed lines represent the $3,6$, and $9\,{\rm km}$ ranges
over which the integrals are computed. Note the good match with the normalised
baryon-number $n_{\rm b}/n_{s}$ at values of $3,2$, and $1$, respectively
(green solid isocontours).

##  Appendix C Extraction Radii of Thermodynamic Quantities

In our analysis of the quark production within the remnant HMNS it seemed
pertinent to not only discuss the amount of quarks produced, but also the
relative location of quark production. To illustrate the extraction radii we
present in Fig. 8 these regions for the Inter_q10 run, for which there is no
collapse to a black hole on the timescale that we cover. Lime contours within
the figure correspond to baryon-number density $n_{\rm b}$ equal to $1n_{\rm
s}$, $2n_{\rm s}$ and $3n_{\rm s}$. The regions enclosed by the black dotted
circles correspond to radii of $3$ km, $6$ km and $9$ km which we designate as
the inner core, outer core and inner crust. The annuli are centred on the
maximal density point for $t<t_{\rm mer}$ and on the centre of mass of the
system for $t\geq t_{\rm mer}$. These extraction domains are then used to
extract the local average and maximal values of thermodynamic quantities as
shown in Fig. 2 and Fig. 7. The averages over any annuli $\Omega$ are obtained
using the following formula:

$\langle f\rangle:=\frac{1}{A}\int_{\Omega}fd^{2}x\,,\qquad
A:=\int_{\Omega}d^{2}x\,,$ (4)

where, for simplicity, we consider the metric determinant of flat spacetime in
the integrand.

## References

  * Baiotti and Rezzolla [2017] L. Baiotti and L. Rezzolla, Rept. Prog. Phys. 80, 096901 (2017), arXiv:1607.03540 [gr-qc] .
  * Paschalidis [2017] V. Paschalidis, Classical and Quantum Gravity 34, 084002 (2017), arXiv:1611.01519 [astro-ph.HE] .
  * Bauswein  _et al._ [2019] A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, Phys. Rev. Lett. 122, 061102 (2019), arXiv:1809.01116 [astro-ph.HE] .
  * Most  _et al._ [2019a] E. R. Most, L. J. Papenfort, V. Dexheimer, M. Hanauske, S. Schramm, H. Stöcker, and L. Rezzolla, Physical Review Letters 122, 061101 (2019a), arXiv:1807.03684 [astro-ph.HE] .
  * Most  _et al._ [2020] E. R. Most, L. Jens Papenfort, V. Dexheimer, M. Hanauske, H. Stoecker, and L. Rezzolla, European Physical Journal A 56, 59 (2020), arXiv:1910.13893 [astro-ph.HE] .
  * Ecker  _et al._ [2020] C. Ecker, M. Järvinen, G. Nijs, and W. van der Schee, Phys. Rev. D 101, 103006 (2020), arXiv:1908.03213 [astro-ph.HE] .
  * Prakash  _et al._ [2021] A. Prakash, D. Radice, D. Logoteta, A. Perego, V. Nedora, I. Bombaci, R. Kashyap, S. Bernuzzi, and A. Endrizzi, Phys. Rev. D 104, 083029 (2021), arXiv:2106.07885 [astro-ph.HE] .
  * Weih  _et al._ [2020] L. R. Weih, M. Hanauske, and L. Rezzolla, Phys. Rev. Lett. 124, 171103 (2020), arXiv:1912.09340 [gr-qc] .
  * Rezzolla and Zanotti [2013] L. Rezzolla and O. Zanotti,  _Relativistic Hydrodynamics_ (Oxford University Press, Oxford, UK, 2013).
  * Demircik, Tuna and Ecker, Christian and Järvinen, Matti [2021] Demircik, Tuna and Ecker, Christian and Järvinen, Matti, arXiv e-prints (2021), arXiv:2112.12157 [hep-ph] .
  * Maldacena [1998] J. M. Maldacena, Adv.Theor.Math.Phys. 2, 231 (1998), arXiv:hep-th/9711200 [hep-th] .
  * Witten [1998] E. Witten, Adv. Theor. Math. Phys. 2, 253 (1998), arXiv:hep-th/9802150 .
  * Jarvinen and Kiritsis [2012] M. Jarvinen and E. Kiritsis, JHEP 03, 002, arXiv:1112.1261 [hep-ph] .
  * Ishii  _et al._ [2019] T. Ishii, M. Järvinen, and G. Nijs, JHEP 07, 003, arXiv:1903.06169 [hep-ph] .
  * Ecker  _et al._ [2020] C. Ecker, M. Järvinen, G. Nijs, and W. van der Schee, Phys. Rev. D 101, 103006 (2020), arXiv:1908.03213 [astro-ph.HE] .
  * Typel  _et al._ [2010] S. Typel, G. Ropke, T. Klahn, D. Blaschke, and H. H. Wolter, Phys. Rev. C 81, 015803 (2010), arXiv:0908.2344 [nucl-th] .
  * Hempel and Schaffner-Bielich [2010] M. Hempel and J. Schaffner-Bielich, Nucl. Phys. A 837, 210 (2010), arXiv:0911.4073 [nucl-th] .
  * Hebeler  _et al._ [2013] K. Hebeler, J. M. Lattimer, C. J. Pethick, and A. Schwenk, Astrophys. J. 773, 11 (2013), arXiv:1303.4662 [astro-ph.SR] .
  * Fraga  _et al._ [2014] E. S. Fraga, A. Kurkela, and A. Vuorinen, Astrophys. J. Lett. 781, L25 (2014), arXiv:1311.5154 [nucl-th] .
  * Cromartie  _et al._ [2019] H. T. Cromartie _et al._ (NANOGrav), Nature Astron. 4, 72 (2019), arXiv:1904.06759 [astro-ph.HE] .
  * Fonseca  _et al._ [2021] E. Fonseca _et al._ , Astrophys. J. Lett. 915, L12 (2021), arXiv:2104.00880 [astro-ph.HE] .
  * Riley  _et al._ [2019] T. E. Riley _et al._ , Astrophys. J. Lett. 887, L21 (2019), arXiv:1912.05702 [astro-ph.HE] .
  * Miller  _et al._ [2019] M. C. Miller _et al._ , Astrophys. J. Lett. 887, L24 (2019), arXiv:1912.05705 [astro-ph.HE] .
  * Nättilä  _et al._ [2017] J. Nättilä, M. C. Miller, A. W. Steiner, J. J. E. Kajava, V. F. Suleimanov, and J. Poutanen, Astron. Astrophys. 608, A31 (2017), arXiv:1709.09120 [astro-ph.HE] .
  * Miller  _et al._ [2021] M. C. Miller _et al._ , Astrophys. J. Lett. 918, L28 (2021), arXiv:2105.06979 [astro-ph.HE] .
  * Riley  _et al._ [2021] T. E. Riley _et al._ , Astrophys. J. Lett. 918, L27 (2021), arXiv:2105.06980 [astro-ph.HE] .
  * Koeppel  _et al._ [2019] S. Koeppel, L. Bovard, and L. Rezzolla, Astrophys. J. Lett. 872, L16 (2019), arXiv:1901.09977 [gr-qc] .
  * Abbott  _et al._ [2017] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 119, 161101 (2017), arXiv:1710.05832 [gr-qc] .
  * Jokela  _et al._ [2021] N. Jokela, M. Järvinen, G. Nijs, and J. Remes, Phys. Rev. D 103, 086004 (2021), arXiv:2006.01141 [hep-ph] .
  * Demircik  _et al._ [2021] T. Demircik, C. Ecker, and M. Järvinen, Astrophys. J. Lett. 907, L37 (2021), arXiv:2009.10731 [astro-ph.HE] .
  * Jokela  _et al._ [2022] N. Jokela, M. Järvinen, and J. Remes, Phys. Rev. D 105, 086005 (2022), arXiv:2111.12101 [hep-ph] .
  * Antoniadis  _et al._ [2013] J. Antoniadis _et al._ , Science 340, 6131 (2013), arXiv:1304.6875 [astro-ph.HE] .
  * Abbott  _et al._ [2018] B. P. Abbott _et al._ (LIGO Scientific, Virgo), Phys. Rev. Lett. 121, 161101 (2018), arXiv:1805.11581 [gr-qc] .
  * Casalderrey-Solana  _et al._ [2014] J. Casalderrey-Solana, H. Liu, D. Mateos, K. Rajagopal, and U. A. Wiedemann,  _Gauge/String Duality, Hot QCD and Heavy Ion Collisions_ (Cambridge University Press, 2014) arXiv:1101.0618 [hep-th] .
  * Busza  _et al._ [2018] W. Busza, K. Rajagopal, and W. van der Schee, Ann. Rev. Nucl. Part. Sci. 68, 339 (2018), arXiv:1802.04801 [hep-ph] .
  * Zaanen  _et al._ [2015] J. Zaanen, Y.-W. Sun, Y. Liu, and K. Schalm,  _Holographic Duality in Condensed Matter Physics_ (Cambridge Univ. Press, 2015).
  * Hartnoll  _et al._ [2016] S. A. Hartnoll, A. Lucas, and S. Sachdev, (2016), arXiv:1612.07324 [hep-th] .
  * Hoyos  _et al._ [2016] C. Hoyos, D. Rodríguez Fernández, N. Jokela, and A. Vuorinen, Phys. Rev. Lett. 117, 032501 (2016), arXiv:1603.02943 [hep-ph] .
  * Bitaghsir Fadafan  _et al._ [2020] K. Bitaghsir Fadafan, J. Cruz Rojas, and N. Evans, Phys. Rev. D 101, 126005 (2020), arXiv:1911.12705 [hep-ph] .
  * Järvinen [2022] M. Järvinen, Eur. Phys. J. C 82, 282 (2022), arXiv:2110.08281 [hep-ph] .
  * Hoyos  _et al._ [2021] C. Hoyos, N. Jokela, and A. Vuorinen, (2021), arXiv:2112.08422 [hep-th] .
  * Kovensky  _et al._ [2022] N. Kovensky, A. Poole, and A. Schmitt, Phys. Rev. D 105, 034022 (2022), arXiv:2111.03374 [hep-ph] .
  * Bartolini  _et al._ [2022] L. Bartolini, S. Bjarke Gudnason, J. Leutgeb, and A. Rebhan, arXiv e-prints , arXiv:2202.12845 (2022), arXiv:2202.12845 [hep-th] .
  * Bigazzi  _et al._ [2005] F. Bigazzi, R. Casero, A. L. Cotrone, E. Kiritsis, and A. Paredes, JHEP 10, 012, arXiv:hep-th/0505140 .
  * Casero  _et al._ [2007] R. Casero, E. Kiritsis, and A. Paredes, Nucl. Phys. B 787, 98 (2007), arXiv:hep-th/0702155 .
  * Gursoy and Kiritsis [2008] U. Gursoy and E. Kiritsis, JHEP 02, 032, arXiv:0707.1324 [hep-th] .
  * Gursoy  _et al._ [2008] U. Gursoy, E. Kiritsis, and F. Nitti, JHEP 02, 019, arXiv:0707.1349 [hep-th] .
  * Veneziano [1976] G. Veneziano, Nucl. Phys. B 117, 519 (1976).
  * Veneziano [1979] G. Veneziano, Nucl. Phys. B 159, 213 (1979).
  * Alho  _et al._ [2013] T. Alho, M. Järvinen, K. Kajantie, E. Kiritsis, and K. Tuominen, JHEP 01, 093, arXiv:1210.4516 [hep-ph] .
  * Alho  _et al._ [2014] T. Alho, M. Järvinen, K. Kajantie, E. Kiritsis, C. Rosen, and K. Tuominen, JHEP 04, 124, [Erratum: JHEP 02, 033 (2015)], arXiv:1312.5199 [hep-ph] .
  * Panero [2009] M. Panero, Phys. Rev. Lett. 103, 232001 (2009), arXiv:0907.3719 [hep-lat] .
  * Alho  _et al._ [2015] T. Alho, M. Jarvinen, K. Kajantie, E. Kiritsis, and K. Tuominen, Phys. Rev. D 91, 055017 (2015), arXiv:1501.06379 [hep-ph] .
  * Jokela  _et al._ [2019] N. Jokela, M. Järvinen, and J. Remes, JHEP 03, 041, arXiv:1809.07770 [hep-ph] .
  * Hoyos  _et al._ [2020] C. Hoyos, N. Jokela, M. Jarvinen, J. G. Subils, J. Tarrio, and A. Vuorinen, Phys. Rev. Lett. 125, 241601 (2020), arXiv:2005.14205 [hep-th] .
  * Hoyos  _et al._ [2022] C. Hoyos, N. Jokela, M. Järvinen, J. G. Subils, J. Tarrio, and A. Vuorinen, Phys. Rev. D 105, 066014 (2022), arXiv:2109.12122 [hep-th] .
  * Note [1] By using the combination of the HS(DD2) and APR models together with V-QCD, we obtain EoSs that agree well with all observational constraints. Using only the HS(DD2) model would lead to a slight tension with the measurement of the tidal deformability from GW170817.
  * Akmal  _et al._ [1998] A. Akmal, V. R. Pandharipande, and D. G. Ravenhall, Phys. Rev. C 58, 1804 (1998), arXiv:nucl-th/9804027 .
  * Annala  _et al._ [2020] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Phys. 16, 907 (2020), arXiv:1903.09121 [astro-ph.HE] .
  * Annala  _et al._ [2021] E. Annala, T. Gorda, E. Katerini, A. Kurkela, J. Nättilä, V. Paschalidis, and A. Vuorinen, arXiv e-prints (2021), arXiv:2105.05132 [astro-ph.HE] .
  * Altiparmak  _et al._ [2022] S. Altiparmak, C. Ecker, and L. Rezzolla, arXiv e-prints , arXiv:2203.14974 (2022), arXiv:2203.14974 [astro-ph.HE] .
  * Papenfort  _et al._ [2021] L. J. Papenfort, S. D. Tootle, P. Grandclément, E. R. Most, and L. Rezzolla, Phys. Rev. D 104, 024057 (2021).
  * Grandclement [2010] P. Grandclement, J. Comput. Phys. 229, 3334 (2010), arXiv:0909.1228 [gr-qc] .
  * Pfeiffer and York [2003] H. P. Pfeiffer and J. W. York, Phys. Rev. D 67, 044022 (2003), gr-qc/0207095 .
  * Pfeiffer and York [2005] H. P. Pfeiffer and J. W. York, Phys. Rev. Lett. 95, 091101 (2005), arXiv:gr-qc/0504142 [gr-qc] .
  * Haas and et al. [2020] R. Haas and et al., The Einstein Toolkit (2020), to find out more, visit http://einsteintoolkit.org.
  * Schnetter  _et al._ [2004] E. Schnetter, S. H. Hawley, and I. Hawke, Class. Quantum Grav. 21, 1465 (2004), gr-qc/0310042 .
  * Most  _et al._ [2019b] E. R. Most, L. J. Papenfort, and L. Rezzolla, Mon. Not. R. Astron. Soc. 490, 3588 (2019b), arXiv:1907.10328 [astro-ph.HE] .
  * Bernuzzi and Hilditch [2010] S. Bernuzzi and D. Hilditch, Phys. Rev. D 81, 084003 (2010), arXiv:0912.2920 [gr-qc] .
  * Alic  _et al._ [2012] D. Alic, C. Bona-Casas, C. Bona, L. Rezzolla, and C. Palenzuela, Phys. Rev. D 85, 064040 (2012), arXiv:1106.2254 [gr-qc] .
  * Etienne  _et al._ [2015] Z. B. Etienne, V. Paschalidis, R. Haas, P. Mösta, and S. L. Shapiro, Class. Quantum Grav. 32, 175009 (2015), arXiv:1501.07276 [astro-ph.HE] .
  * Del Zanna  _et al._ [2007] L. Del Zanna, O. Zanotti, N. Bucciantini, and P. Londrillo, Astron. Astrophys. 473, 11 (2007), arXiv:0704.3206 .
  * Gill  _et al._ [2019] R. Gill, A. Nathanail, and L. Rezzolla, Astrophys. J. 876, 139 (2019), arXiv:1901.04138 [astro-ph.HE] .
  * Murguia-Berthier  _et al._ [2020] A. Murguia-Berthier, E. Ramirez-Ruiz, F. De Colle, A. Janiuk, S. Rosswog, and W. H. Lee, arXiv e-prints , arXiv:2007.12245 (2020), arXiv:2007.12245 [astro-ph.HE] .
  * Weih  _et al._ [2020] L. R. Weih, M. Hanauske, and L. Rezzolla, Phys. Rev. Lett. 124, 171103 (2020), arXiv:1912.09340 [gr-qc] .
  * Barsotti  _et al._ [2018] L. Barsotti, S. Gras, M. Evans, and P. Fritschel, LIGO Document No. T1800044-v5 (2018).
  * Punturo  _et al._ [2010] M. Punturo _et al._ ,  _Gravitational waves. Proceedings, 8th Edoardo Amaldi Conference, Amaldi 8, New York, USA, June 22-26, 2009_ , Class. Quant. Grav. 27, 084007 (2010).
  * Alcubierre [2006] M. Alcubierre, Introduction to 3+1 Numerical Relativity (2006).
  * Bishop and Rezzolla [2016] N. T. Bishop and L. Rezzolla, Living Reviews in Relativity 19, 2 (2016), arXiv:1606.02532 [gr-qc] .
  * Ghirlanda  _et al._ [2019] G. Ghirlanda _et al._ , Science 363, 968 (2019), arXiv:1808.00469 [astro-ph.HE] .
  * Takami  _et al._ [2014] K. Takami, L. Rezzolla, and L. Baiotti, Phys. Rev. Lett. 113, 091104 (2014), arXiv:1403.5672 [gr-qc] .
  * Sekiguchi  _et al._ [2011] Y. Sekiguchi, K. Kiuchi, K. Kyutoku, and M. Shibata, Phys. Rev. Lett. 107, 211101 (2011), arXiv:1110.4442 [astro-ph.HE] .
  * Radice  _et al._ [2017] D. Radice, S. Bernuzzi, W. Del Pozzo, L. F. Roberts, and C. D. Ott, Astrophys. J. Lett. 842, L10 (2017), arXiv:1612.06429 [astro-ph.HE] .
  * Tootle  _et al._ [2021] S. D. Tootle, L. J. Papenfort, E. R. Most, and L. Rezzolla, Astrophys. J. Lett. 922, L19 (2021), arXiv:2109.00940 [gr-qc] .
  * Kashyap  _et al._ [2021] R. Kashyap, A. Das, D. Radice, S. Padamata, A. Prakash, D. Logoteta, A. Perego, D. A. Godzieba, S. Bernuzzi, I. Bombaci, F. J. Fattoyev, B. T. Reed, and A. da Silva Schneider, arXiv e-prints , arXiv:2111.05183 (2021), arXiv:2111.05183 [astro-ph.HE] .
  * Kölsch  _et al._ [2021] M. Kölsch, T. Dietrich, M. Ujevic, and B. Bruegmann, arXiv e-prints , arXiv:2112.11851 (2021), arXiv:2112.11851 [gr-qc] .
  * Fujibayashi  _et al._ [2020] S. Fujibayashi, S. Wanajo, K. Kiuchi, K. Kyutoku, Y. Sekiguchi, and M. Shibata, Astrophys. J. 901, 122 (2020), arXiv:2007.00474 [astro-ph.HE] .
  * Bauswein  _et al._ [2021] A. Bauswein, S. Blacker, G. Lioutas, T. Soultanis, V. Vijayan, and N. Stergioulas, Phys. Rev. D 103, 123004 (2021), arXiv:2010.04461 [astro-ph.HE] .
  * Bozzola [2021] G. Bozzola, The Journal of Open Source Software 6, 3099 (2021), arXiv:2104.06376 [gr-qc] .