Maximum Mass Of Differentially Rotating Strange Quark Stars

TL;DR

This paper presents relativistic simulations showing that differential rotation can significantly increase the maximum mass of strange quark stars, with potential implications for astrophysical observations of massive compact objects.

Contribution

First fully relativistic numerical calculations of differentially rotating strange quark stars using the MIT Bag model, exploring their maximum mass limits across various rotation degrees.

Findings

Maximum mass increases with differential rotation and solution type.

Maximum mass can be over four times the non-rotating value.

Differential rotation allows strange stars to reach masses much larger than nuclear matter stars.

Abstract

We present the first fully relativistic numerical calculations of differentially rotating Strange Quark Stars models for broad ranges of the maximum density and of the degree of differential rotation. Our simulations are performed with the very accurate and stable multi-domain spectral code FlatStar and use the MIT Bag model for describing strange quark matter. Our calculations based on a thorough exploration of the solution space show that the maximum mass of strange stars depends on both the degree of differential rotation and a type of solution, similarly to neutron stars described by a polytropic equations of state. The highest increase of the maximum mass (compared to the value for a non-rotating star) is obtained for models with a low degree of differential rotation. This highest mass is over four times larger than that of the equivalent non rotating configuration. Comparing our results with calculations done for realistic models of neutron stars, we conclude that with the help of differential rotation, strange stars can sustain masses much larger than stars made from nuclear matter for low degree of differential rotation which reinforces the hope of demonstrating, or of ruling out, the existence of strange matter through the observation by gravitational waves, by gamma rays or by neutrinos of the massive material object born from the merger of a compact binary system or during some supernova events.

Figures

Tables

Table 1. Table 1: Comparison of three codes (FlatStar, RotSeq and RNS) for non-rotating and rigidly rotating Strange Quark Stars [RotSeq - Gourgoulhon et al. ( 1999 ) ; RNS - Stergioulas et al. ( 1999 ) ]. For all configurations, the Table displays the gravitational mass M 𝑀 M , the baryonic mass M B subscript 𝑀 B M_{\mathrm{B}} , the angular velocity Ω Ω \Omega , the central mass density ρ c subscript 𝜌 c \rho_{\mathrm{c}} , the central enthalpy H c subscript 𝐻 c H_{\mathrm{c}} , the circumferential radius R circ subscript 𝑅 circ R_{\mathrm{circ}} , the coordinate equatorial radius r e subscript 𝑟 𝑒 r_{e} , the ratio between the polar and equatorial radii r ratio = r p / r e subscript 𝑟 ratio subscript 𝑟 𝑝 subscript 𝑟 𝑒 r_{\mathrm{ratio}}=r_{p}/r_{e} , the angular momentum J 𝐽 J and the ratio between the kinetic and potential energies T / | W | 𝑇 𝑊 T/|W| . The upper part of the Table contains static configurations ( Ω = 0 Ω 0 \Omega=0 ) with maximum mass, while the lower part presents rigidly rotating configurations with a central mass density ρ r ​ m ​ c subscript 𝜌 𝑟 𝑚 𝑐 \rho_{rmc} such that the star is the one with maximum mass for the RotSeq code. As explained in the text, those stars are not exactly those of maximum mass for RNS and FlatStar.

	$M\left[M_{\odot }\right]$	$M_B\left[M_{\odot }\right]$	$\Omega \left[\mathrm{rad}{s}^{-1}\right]$	${\rho }_c\left[{10}^{15}g{\mathrm{cm}}^{-3}\right]$	$H_c$	$R_{\mathrm{circ}}\left[\mathrm{km}\right]$	$r_e\left[\mathrm{km}\right]$	$r_{\mathrm{ratio}}$	$J\left[G{M}_{\odot }^2/c\right]$	$T/\left|W\right|$
$M_{\max }^{\mathrm{stat}}$
FlatStar	1.9637	2.6489	0	2.059	0.4514	10.71	7.532	1	0	0
RotSeq	1.9638	2.6255	0	2.059	0.4514	10.71	7.532	1	0	0
RNS	1.9676	2.6330	0	2.059	0.4514	10.71	7.527	1	0	0
$M_{\max }^{\mathrm{rig}}$
FlatStar	2.8188	3.778	9560.19	1.261	0.32	16.463	11.308	0.4643	7.013	0.209
RotSeq	2.8310	3.751	9547.0	1.261	0.32	16.54	11.37	0.4618	7.084	0.21
RNS	2.8339	3.759	9562.61	1.261	0.32	16.425	11.263	0.466	7.093	0.21

Table 2. Table 2: Comparison of differentially rotating Strange Quark Stars obtained with FlatStar and with LORENE, for two values of r ratio subscript 𝑟 ratio r_{\mathrm{ratio}} ( 0.9 0.9 0.9 and 0.7 0.7 0.7 ) as well as of A ~ ~ 𝐴 \tilde{A} ( 0.1 0.1 0.1 and 0.2 0.2 0.2 ), and for fixed ρ max = 0.82 ⋅ 10 15 ​ g / cm 3 subscript 𝜌 max ⋅ 0.82 superscript 10 15 g superscript cm 3 \rho_{\mathrm{max}}=0.82\cdot 10^{15}\mathrm{g/cm^{3}} ( H max = 0.2 subscript 𝐻 max 0.2 H_{\mathrm{max}}=0.2 ). In the Table are displayed the gravitational mass M 𝑀 M , the baryonic mass M B subscript 𝑀 𝐵 M_{B} , the circumferential radius R c ​ i ​ r ​ c subscript 𝑅 𝑐 𝑖 𝑟 𝑐 R_{circ} , the angular momentum J 𝐽 J , the central angular velocity Ω c subscript Ω 𝑐 \Omega_{c} , the corresponding central rotational frequency f c = Ω c / ( 2 ​ π ) subscript 𝑓 𝑐 subscript Ω 𝑐 2 𝜋 f_{c}=\Omega_{c}/(2\pi) and the ratio between the kinetic and potential energies T / | W | 𝑇 𝑊 T/|W| .

		$\tilde{A}=0.1$		$\tilde{A}=0.2$
		FlatStar	LORENE	FlatStar	LORENE
	$M\left[M_{\odot }\right]$	1.61758	1.61767	1.61839	1.61853
	$M_B\left[M_{\odot }\right]$	2.09699	2.07821	2.09828	2.07954
	$R_{\mathrm{circ}}\left[\mathrm{km}\right]$	11.4927	11.4938	11.4857	11.4871
$r_{\mathrm{ratio}}=0.9$	$J\left[G{M}_{\odot }^2/c\right]$	0.93263	0.93272	0.93015	0.93051
	${\Omega }_c\left[\mathrm{rad}{s}^{-1}\right]$	4178.87	4178.44	4327.73	4328.22
	$f_c\left[Hz\right]$	665.087	665.019	688.779	688.858
	$T/\left|W\right|$	0.03338	0.03338	0.03317	0.03318
	$M\left[M_{\odot }\right]$	1.90198	1.90213	1.90668	1.90681
	$M_B\left[M_{\odot }\right]$	2.47319	2.45111	2.48060	2.45835
	$R_{\mathrm{circ}}\left[\mathrm{km}\right]$	12.9953	12.9969	12.9561	12.9585
$r_{\mathrm{ratio}}=0.7$	$J\left[G{M}_{\odot }^2/c\right]$	2.38766	2.38804	2.38863	2.38921
	${\Omega }_c\left[\mathrm{rad}{s}^{-1}\right]$	7111.69	7110.62	7421.57	7421.32
	$f_c\left[Hz\right]$	1131.86	1131.69	1181.18	1181.14
	$T/\left|W\right|$	0.11146	0.11145	0.11082	0.11086

Table 3. Table 3: Same as in Tab.  2 but for ρ max = 1.70 ⋅ 10 15 ​ g / cm 3 subscript 𝜌 max ⋅ 1.70 superscript 10 15 g superscript cm 3 \rho_{\mathrm{max}}=1.70\cdot 10^{15}\mathrm{g/cm^{3}} ( H max = 0.4 subscript 𝐻 max 0.4 H_{\mathrm{max}}=0.4 ).

		$\tilde{A}=0.1$		$\tilde{A}=0.2$
		FlatStar	LORENE	FlatStar	LORENE
	$M\left[M_{\odot }\right]$	2.05075	2.05085	2.05070	2.05081
	$M_B\left[M_{\odot }\right]$	2.76025	2.73548	2.76040	2.73566
	$R_{\mathrm{circ}}\left[\mathrm{km}\right]$	11.3862	11.3873	11.3710	11.3722
$r_{\mathrm{ratio}}=0.9$	$J\left[G{M}_{\odot }^2/c\right]$	1.49360	1.49381	1.48696	1.48721
	${\Omega }_c\left[\mathrm{rad}{s}^{-1}\right]$	5250.26	5250.07	5510.20	5510.08
	$f_c\left[Hz\right]$	835.604	835.575	876.976	876.956
	$T/\left|W\right|$	0.03501	0.03501	0.03476	0.03477
	$M\left[M_{\odot }\right]$	2.34800	2.34792	2.35005	2.34975
	$M_B\left[M_{\odot }\right]$	3.15643	3.12785	3.16054	3.13140
	$R_{\mathrm{circ}}\left[\mathrm{km}\right]$	12.7302	12.7315	12.6439	12.6448
$r_{\mathrm{ratio}}=0.7$	$J\left[G{M}_{\odot }^2/c\right]$	3.58145	3.58078	3.57413	3.57317
	${\Omega }_c\left[\mathrm{rad}{s}^{-1}\right]$	8894.29	8892.15	9486.42	9486.16
	$f_c\left[Hz\right]$	1415.57	1415.23	1509.81	1509.77
	$T/\left|W\right|$	0.11547	0.11544	0.11516	0.11519

Table 4. Table 4: Properties of configurations with maximum mass of non-rotating ( A ~ = 0.0 ~ 𝐴 0.0 \tilde{A}=0.0 , Ω = 0.0 Ω 0.0 \Omega=0.0 ), rigidly rotating ( A ~ = 0.0 ~ 𝐴 0.0 \tilde{A}=0.0 - Keplerian limit) and differentially rotating ( A ~ = 0.1 , 0.2 , 0.3 , 0.5 , 0.7 ~ 𝐴 0.1 0.2 0.3 0.5 0.7 \tilde{A}=0.1,0.2,0.3,0.5,0.7 ) Strange Quark Stars of types A and C. In the Table are given the baryonic mass M B subscript 𝑀 B M_{\mathrm{B}} , the gravitational mass M 𝑀 M , the central and equatorial angular velocities Ω c subscript Ω 𝑐 \Omega_{c} and Ω eq subscript Ω eq \Omega_{\mathrm{eq}} , the corresponding frequencies f c subscript 𝑓 𝑐 f_{c} and f eq subscript 𝑓 eq f_{\mathrm{eq}} , the central and maximal densities ρ c subscript 𝜌 𝑐 \rho_{c} and ρ max subscript 𝜌 max \rho_{\mathrm{max}} , the corresponding enthalpies H c subscript 𝐻 𝑐 H_{c} and H max subscript 𝐻 max H_{\mathrm{max}} , the circumferential radius R circ subscript 𝑅 circ R_{\mathrm{circ}} , the coordinate equatorial radius r e subscript 𝑟 𝑒 r_{e} , the ratio between the polar and equatorial radii r ratio = r p / r e subscript 𝑟 ratio subscript 𝑟 𝑝 subscript 𝑟 𝑒 r_{\mathrm{ratio}}=r_{p}/r_{e} , the angular momentum J 𝐽 J , and the ratio between the kinetic and potential energies T / | W | 𝑇 𝑊 T/|W| .

	$\tilde{A}=0.0$		$\tilde{A}=0.0$	$\tilde{A}=0.1$	$\tilde{A}=0.2$		$\tilde{A}=0.3$	$\tilde{A}=0.5$	$\tilde{A}=0.7$
	$\Omega =0.0$		Keplerian limit	Type A			Type C
$M\left[M_{\odot }\right]$	1.96372		2.81847	2.91624	3.26356		5.25630	4.24179	3.67839
$M_B\left[M_{\odot }\right]$	2.64888		3.78386	3.90858	4.38617		7.25072	5.77867	4.97946
${\Omega }_c\left[\mathrm{rad}{s}^{-1}\right]$	0		9655.09	10066.3	11037.5		13703.2	19512.4	26019.3
${\Omega }_{\mathrm{eq}}\left[\mathrm{rad}{s}^{-1}\right]$	0		9655.09	9561.17	9039.06		7669.96	7712.25	7555.35
$f_c\left[Hz\right]$	0		1536.66	1602.09	1756.67		2180.93	3105.50	4141.09
$f_{\mathrm{eq}}\left[Hz\right]$	0		1536.66	1521.71	1438.61		1220.71	1227.44	1202.47
${\rho }_c\left[{10}^{15}g{\mathrm{cm}}^{-3}\right]$	2.059		1.323	1.362	1.261		0.428	0.428	0.428
${\rho }_{\max }\left[{10}^{15}g{\mathrm{cm}}^{-3}\right]$	2.059		1.323	1.362	1.261		0.690	0.779	0.866
$H_c$	0.451396		0.33308	0.34105	0.32008		0.00043	0.00033	0.00023
$H_{\max }$	0.451396		0.33308	0.34105	0.32008		0.14913	0.18473	0.21510
$R_{\mathrm{circ}}\left[\mathrm{km}\right]$	10.7104		16.3338	16.5381	17.7295		21.9281	18.5692	16.6423
$r_e\left[\mathrm{km}\right]$	7.53167		11.17434	11.1727	11.5620		10.6879	9.97624	9.46258
$r_{\mathrm{ratio}}$	1		0.46792	0.45191	0.38979		0.01	0.01	0.01
$J\left[G{M}_{\odot }^2/c\right]$	0		6.96846	7.59663	10.1705		28.1603	17.4939	12.4682
$T/\left|W\right|$	0		0.20721	0.21652	0.25269		0.33465	0.30538	0.27652

Table 5. Table 5: Properties of configurations with maximum mass of differentially rotating ( A ~ ~ 𝐴 \tilde{A} = 0.015, 0.2, 0.23, 0.25, 0.27, 0.29, 0.3, 0.33, 0.35) Strange Quark Stars of types B and D. In the Table are given the baryonic mass M B subscript 𝑀 B M_{\mathrm{B}} , the gravitational mass M 𝑀 M , the central and equatorial angular velocities Ω c subscript Ω 𝑐 \Omega_{c} and Ω eq subscript Ω eq \Omega_{\mathrm{eq}} , the corresponding frequencies f c subscript 𝑓 𝑐 f_{c} and f eq subscript 𝑓 eq f_{\mathrm{eq}} , the central and maximal densities ρ c subscript 𝜌 𝑐 \rho_{c} and ρ max subscript 𝜌 max \rho_{\mathrm{max}} , the corresponding enthalpies H c subscript 𝐻 c H_{\mathrm{c}} and H max subscript 𝐻 max H_{\mathrm{max}} , the circumferential radius R circ subscript 𝑅 circ R_{\mathrm{circ}} , the coordinate equatorial radius r e subscript 𝑟 𝑒 r_{e} , the ratio between the polar and equatorial radii r ratio = r p / r e subscript 𝑟 ratio subscript 𝑟 𝑝 subscript 𝑟 𝑒 r_{\mathrm{ratio}}=r_{p}/r_{e} , the angular momentum J 𝐽 J , and the ratio between the kinetic and potential energies T / | W | 𝑇 𝑊 T/|W| .

	$\tilde{A}=0.015$	$\tilde{A}=0.2$	$\tilde{A}=0.23$	$\tilde{A}=0.25$	$\tilde{A}=0.27$	$\tilde{A}=0.29$		$\tilde{A}=0.3$	$\tilde{A}=0.33$	$\tilde{A}=0.35$
	Type B							Type D
$M\left[M_{\odot }\right]$	7.59282	6.92478	6.48297	6.16558	5.83180	5.47377		5.27581	4.59474	4.01156
$M_B\left[M_{\odot }\right]$	10.8074	9.62320	8.88792	8.37731	7.84149	7.26888		6.96280	5.93151	5.08301
${\Omega }_c\left[\mathrm{rad}{s}^{-1}\right]$	9447.75	9271.92	9181.15	9054.78	8882.55	8634.43		8500.81	7937.89	7427.45
${\Omega }_{\mathrm{eq}}\left[\mathrm{rad}{s}^{-1}\right]$	6823.31	6443.00	6269.71	6144.96	6014.84	5865.32		5791.96	5516.91	5286.19
$f_c\left[Hz\right]$	1503.66	1475.67	1461.23	1441.11	1413.70	1374.21		1352.95	1263.36	1182.12
$f_{\mathrm{eq}}\left[Hz\right]$	1085.964	1025.44	997.855	978.001	957.292	933.495		921.818	878.043	841.323
${\rho }_c\left[{10}^{15}g{\mathrm{cm}}^{-3}\right]$	0.462	0.510	0.529	0.535	0.537	0.534		0.531	0.512	0.495
${\rho }_{\max }\left[{10}^{15}g{\mathrm{cm}}^{-3}\right]$	0.588	0.552	0.543	0.539	0.537	0.534		0.531	0.512	0.495
$H_c$	0.025510	0.05689	0.06837	0.07224	0.07340	0.07120		0.06940	0.05850	0.04760
$H_{\max }$	0.10100	0.08200	0.07650	0.07409	0.07340	0.07120		0.06940	0.05850	0.04760
$R_{\mathrm{circ}}\left[\mathrm{km}\right]$	29.9610	30.0117	29.8564	29.7765	29.6938	29.6490		29.5937	29.4891	29.3296
$r_e\left[\mathrm{km}\right]$	12.9534	15.1022	16.2005	16.9673	17.7549	18.6191		19.0441	20.5223	21.6422
$r_{\mathrm{ratio}}$	0.06203	0.09897	0.11268	0.11949	0.12494	0.12848		0.12997	0.12994	0.12594
$J\left[G{M}_{\odot }^2/c\right]$	60.6172	51.7050	46.0668	42.2461	38.3987	34.5055		32.4296	25.8446	20.7697
$T/\left|W\right|$	0.36215	0.36655	0.36815	0.36915	0.37021	0.37174		0.37252	0.37675	0.38172

Table 6. Table 6: Maximum baryonic mass increase (in %, compared to the static configuration) for differentially rotating stars described by different EOSs: Strange Quark Stars (SQS) and six realistic nuclear matter EOS (A, D, L, UT, FPS, APR) taken from Morrison et al. ( 2004 ) . The corresponding value of r ratio subscript 𝑟 ratio r_{\mathrm{ratio}} for the configuration with maximum mass is indicated between parenthesis. See the precise references of the articles that describe the nuclear EOS in Morrison et al. ( 2004 ) .

	$M_{B,\max }$	$\delta M_B$ ($r_p/r_e$)	$\delta M_B$ ($r_p/r_e$)	$\delta M_B$ ($r_p/r_e$)	$\delta M_B$ ($r_p/r_e$)	$\delta M_B$ ($r_p/r_e$)	$\delta M_B$ ($r_p/r_e$)	$\delta M_B$ ($r_p/r_e$)
	$\left[M_{\odot }\right]$	[%]	[%]	[%]	[%]	[%]	[%]	[%]
EoS	$\Omega =0.0$	$\tilde{A}=0.0$	$\tilde{A}=0.1$	$\tilde{A}=0.2$	$\tilde{A}=0.3$	$\tilde{A}=0.5$	$\tilde{A}=0.7$	$\tilde{A}=1.0$
SQS	2.65	43 (0.468)	48 (0.452)	263 (0.10)	174 (0.01)	118 (0.01)	88 (0.01)	-
A aaA - Reid soft core (Pandharipande, 1971).	1.92	17 (0.565)	-	-	30 (0.46)	40 (0.415)	51 (0.22)	22 (0.31)
D bbD - Model V (Bethe & Johnson, 1974).	1.89	18 (0.565)	-	-	28 (0.505)	56 (0.385)	60 (0.275)	57 (2.965)
L ccL - Mean field (Pandharipande & Smith, 1975).	3.23	20 (0.55)	-	-	40 (0.425)	48 (0.32)	56 (0.01)	26 (0.025)
UT ddUT - UV14 + TNI (Wiringa et al., 1988).	2.17	18 (0.565)	-	-	32 (0.46)	34 (0.49)	50 (0.22)	37 (0.01)
FPS eeFPS - UV14 + TNI (Lorenz et al., 1993).	2.10	17 (0.565)	-	-	29 (0.475)	46 (0.36)	45 (0.275)	42 (0.01)
APR ffAPR - A18 + $\delta \upsilon$ + UIX* (Akmal et al., 1998).	2.67	16 (0.58)	-	-	31 (0.445)	25 (0.37)	19 (0.28)	17 (0.01)

Equations

$$p=a{c}^2\left(\rho -{\rho }_0\right),$$

$$dH=\frac{dp}{c^2\rho (p)+p},$$

$$\rho (H)=\frac{{\rho }_0}{1+a}\left(e^{H(a+1)/a}+a\right),$$

Full text

equation of state — gravitation — methods: numerical — stars: neutron, quark, strange — stars: rotation

††journal: ApJ††software: FlatStar for rigid rotation (Ansorg et al., 2003a, b; Schöbel & Ansorg, 2003; Ansorg, 2005; Ansorg & Petroff, 2005), FlatStar for differential rotation (Ansorg et al., 2009; Gondek-Rosińska et al., 2017), RotSeq (Gourgoulhon et al., 2016, http://www.lorene.obspm.fr), RNS (Stergioulas & Friedman, 1995)

1 Introduction

The first detection of gravitational waves (GW) from coalescing neutron stars in a binary system, GW170817, by the Advanced Virgo and LIGO collaboration (Abbott et al., 2017) proved that we have a new tool to investigate the properties of matter at extreme densities and, as a consequence, to put constraints on the equation of state (EOS) of compact stars, since both the electromagnetic and the GW signals are strongly affected by the specific nature of their source. The merger of two neutron stars could lead to the formation of a hot, massive differentially rotating remnant-a neutron star or even a strange quark star (Drago et al., 2016; Most et al., 2019) or to a black hole. The outcome mainly depends on the EOS, the total mass and the mass ratio of compact stars in a binary system (see Baiotti & Rezzolla, 2017, for a review) and has important implications on the gravitational wave signal of such a merger. Many attempts have been made to interpret the observations of GW170817 (Rezzolla et al., 2018; Ruiz et al., 2018; Shibata & Kiuchi, 2017; Shibata et al., 2017; Most et al., 2019) but there is still no consensus what is the fate of the post-merger object.

In our paper we investigate properties of differentially rotating strange quark stars, in particular we focus on their maximum masses.

The possibility of the existence of stable deconfined quark matter was indeed predicted in the early seventies by Bodmer (1971). Later, the hypothesis that strange (quark) matter, with a non-vanishing strangeness per unit baryon number, could constitute the absolute ground state of matter at zero pressure and temperature was made famous by Witten (1984). He pointed out that three-quark-flavored baryonic matter could be energetically favorable compared to the two-flavored one, meaning that at zero pressure strange quark matter could have a lower energy per baryon than iron ${}^{56}\mathrm{Fe}$. If that assumption was correct, normal nuclei could decrease their energy by transforming into strange matter (e.g. Bombaci & Datta, 2000). Among the potential environments that could result in the actual formation of strange quark matter, the early Universe and the interior of compact stars are the most likely, which naturally led to the theoretical study of strange nuclei called strangelets (Farhi & Jaffe, 1984) but also of strange quark stars (SQS) entirely built of strange quark matter with a small addition of electrons (Alcock et al., 1986; Haensel et al., 1986).

Due to its composition quite different at the microscopic level from that of a usual neutron star (NS), a SQS is macroscopically described by a rather distinct EOS, which implies easily distinguishable properties. For instance, while NSs have a strong density contrast between their center and their surface, the density profile of SQSs is almost constant with a non-vanishing surface value. This feature means that even if static SQSs and NSs can sustain a quite similar maximum mass (whose specific value depends on the precise EOS, see Chamel et al., 2013, for a review), it is no longer the case when the stars are rotating. Indeed, rotation can stabilize stars with masses larger than the maximum mass of non-rotating star $M_{\mathrm{max,stat}}$ (Baumgarte et al., 2000) and even leads to the occurrence of configurations that cannot exist in the static case. However, while it was shown that for rigidly rotating NSs the maximum allowed mass can be 14% - 22% larger than $M_{\mathrm{max,stat}}$ (Cook et al., 1992, 1994a, 1994b), it is up to 44% larger for SQSs (Gourgoulhon et al., 1999; Gondek-Rosińska et al., 2000).

Such a difference can be crucial in deciding whether SQSs exist, and consequently in learning more about nuclear matter at high density, as the value of the maximum mass of a relativistic star helps to determine if an observed compact object is a black hole or a material celestial body. Moreover, this quantity is also a key-factor for establishing the life span of the short-lived remnant born from the merger of a compact binary system or of a hot neutron star born from a supernova explosion, situations in which the estimation of the maximum mass is further complicated by the fact that such objects are differentially rotating (Baumgarte et al., 2000; Shibata & Uryū, 2002; Shibata et al., 2005; Giacomazzo et al., 2011).

The maximum mass of differentially rotating neutron stars has already been the topic of many studies (e.g. Baumgarte et al., 2000; Lyford et al., 2003; Morrison et al., 2004; Bozzola et al., 2018; Weih et al., 2018; Uryū et al., 2017), but was more recently approached with the relativistic, highly accurate double-domain pseudo-spectral FlatStar code (based on the so-called “AKM-method”, Ansorg et al., 2003a) which led to several noticeable developments. Indeed, the AKM-method, which was formerly shown to enable the calculation of very extremal configurations of rigidly rotating relativistic stars (Ansorg et al., 2003b; Schöbel & Ansorg, 2003) or rings (Ansorg, 2005; Ansorg & Petroff, 2005), allowed a deeper understanding of the general structure of the solution space for constant density stars and $N=1$ (equivalently $\Gamma=2$) polytropes (Ansorg et al., 2009), demonstrating for instance the existence of new types of configurations that were not considered in previous studies. One of the key features of the code that made such a thorough exploration possible is the resort to the Newton-Raphson method with the help of which any physical quantity can be used either as a free parameter to vary or as a fixed one.

This property turned out to be crucial for the works Gondek-Rosińska et al. (2017) and Studzińska et al. (2016) where for the first time, the maximum mass and various other astrophysical parameters were calculated for all types of differentially rotating neutron stars modelized by polytropes of various adiabatic indices. It was found that the maximum mass of differentially rotating neutron stars depends on the stiffness of the equation of state, on the degree of differential rotation but also on the type of configuration. More precisely, the highest increase with respect to the maximum mass for non rotating stars with the same equation of state was reached for a modest degree of differential rotation and a moderate stiffness.

In this article, we present the first fully relativistic calculation of the structure of stationary and axisymmetric models of differentially rotating strange quark stars for broad ranges of the maximum density and of the degree of differential rotation, focusing on the determination of their maximum mass. The paper is organized as follows. First, we briefly describe the numerical model (chosen EOS, rotation law, etc.), sending the reader to previous articles (Ansorg et al., 2009; Gondek-Rosińska et al., 2017) for more details on the algorithms. Then, in Section 3, we summarize some numerical tests (consistency checks and comparisons with other codes) that were done to validate our results. The latter are displayed in Section 4, emphasis being put on the comparison with the values obtained for realistic EOSs of NSs by Morrison et al. (2004). Finally, Section 5 provides a summary of our work and a discussion of its possible extensions.

2 Numerical model and assumptions

The FlatStar code, based on spectral methods and on the Newton-Raphson algorithm (Ansorg et al., 2009; Gondek-Rosińska et al., 2017; Studzińska et al., 2016), relies on the assumptions that spacetime and the distribution of matter are axisymmetric and stationary. Additionally, the star is supposed to be made from a perfect fluid which, as far as thermodynamics is concerned, is fully characterized by an EOS giving the pressure in the fluid frame $p$ as a function of the total energy density in that same frame $\epsilon$ (or equivalently the total mass density $\rho=\epsilon/c^{2}$, with $c$ the speed of light).

In the case of strange matter, it was shown (Gondek-Rosińska et al., 2000; Zdunik, 2002) that, even if the physics of quarks at high density is highly non-linear and complicated, the EOS could be very precisely approximated by a linear relation

[TABLE]

where $a$ and $\rho_{0}$ are constants such that $c\sqrt{a}$ is the speed of sound and $\rho_{0}$ is the density at zero pressure (and as a consequence also at the surface of the SQS in the absence of a crust, as assumed in the following). In this work, as a first step in the study of differentially rotating SQSs, we only consider the simplest form of the MIT Bag model, in which $a=1/3$ and $\rho_{0}=4.2785\times 10^{14}\mathrm{g/cm^{3}}$.

Once given a barotropic EOS $p=p(\rho)$, a configuration of a rigidly rotating relativistic star is uniquely determined by solving the axisymmetric and stationary Einstein equations taking into account appropriate boundary conditions but also providing a central (or maximal) mass (or energy) density $\rho$ and an angular velocity $\Omega$. In the case of differential rotation, the situation is more intricate and one can freely choose a rotation law relating a kind of specific relativistic angular momentum to the local angular velocity. As in Ansorg et al. (2009); Studzińska et al. (2016); Gondek-Rosińska et al. (2017), we adopt here the simple but realistic and astronomically motivated law introduced by Komatsu et al. (1989) and considered by many authors (among which Cook et al., 1992; Baumgarte et al., 2000; Lyford et al., 2003). With that law, a rotation profile is specified by $\Omega_{c}$, the limit of the angular velocity close to the rotation axis, and a parameter which describes the degree of differential rotation. In the present article as in previous ones, we follow Baumgarte et al. (2000) and use for the latter the dimensionless parameter $\tilde{A}$, whose exact definition can be found in Ansorg et al. (2009). Notice that

• •

the Komatsu et al. (1989) law implies an angular velocity that monotonously decreases from the axis of rotation to the equator;

• •

$\tilde{A}=0$ corresponds to rigid rotation, $\Omega=\Omega_{c}$ in the whole star;

• •

for $\tilde{A}=1$ the angular velocity at the equatorial surface is around half $\Omega_{c}$;

• •

the degree of differential rotation is an increasing and monotonic function of $\tilde{A}$.

Another key difference between rigidly and differentially stationary rotating stars is that the solution space for the latters is much richer. Hence, to uniquely specify a configuration, it is not always sufficient to set the value of the (e.g.) central density, and some geometrical quantities turn out to be quite convenient. In the following, we shall refer to the ratio between the polar and equatorial radii $r_{\mathrm{ratio}}=r_{p}/r_{e}$, with $0<r_{\mathrm{ratio}}\leq 1$ as we consider only stars without a hole at their center, and also to the so-called ”rescaled shedding parameter” $\tilde{\beta}$, see the precise definition in Ansorg et al. (2009), which is such that $0<\tilde{\beta}<1$, with $\tilde{\beta}\rightarrow 1$ when a star enters into the toroidal regime ($r_{\mathrm{ratio}}=0$), $\tilde{\beta}\rightarrow 0$ at the mass-shedding limit and $\tilde{\beta}=1/2$ for a non-rotating spherical star ($r_{\mathrm{ratio}}=1$).

The last technical point worth mentioning is that the FlatStar code uses as a primary thermodynamical variable not the mass density $\rho$, but the corresponding dimensionless relativistic enthalpy $H(\rho)$ defined by the differential relation

[TABLE]

with the property $H(p=0)=0$. As a consequence, in the following our results are equivalently described for given values of the maximum mass density in the star, $\rho_{\mathrm{max}}$, or of the maximum enthalpy $H_{\mathrm{max}}$. Notice that with the EOS (1), one has the relation

[TABLE]

so that $H=0$ corresponds to $\rho=\rho_{0}$, another consequence of which being that in the following sections, a linear scale for $H$ on one axis of a figure is equivalent to a logarithmic-like scale for $\rho$ on the same axis.

3 Numerical tests

3.1 FlatStar testing

The first internal consistency test of the code we did simply consisted in verifying, for rigidly rotating stars, the agreement of calculated configurations with the first law of thermodynamics for relativistic stars (see e.g. Friedman & Stergioulas, 2013), which implies that in a sequence with fixed angular momentum $J$ the maxima of gravitational $M$ and baryonic $M_{B}$ masses occur for the same value of the maximum density $\rho_{\mathrm{max}}$ (equivalently of the maximum enthalpy $H_{\mathrm{max}}$). On Figure 1, we display the results for that test in the case $J=2.16\ [G\mathrm{M_{\odot}^{2}}/c]$. It can be seen that indeed, with a high degree of precision, the maxima of $M_{B}$ and of $M$ are reached for the same maximum enthalpy $H_{\mathrm{max}}=0.43971$.

The second test we performed concerned the dependence of the accuracy of the code on the number of grid points [see also the Appendix in Gondek-Rosińska et al. (2017)] for differentially rotating stars. One way to proceed was to search for the maximum mass of stars with a fixed degree of differential rotation ($\tilde{A}$) and to check how its value changes with an increasing number of grid points. The results for SQSs with $\tilde{A}=0.1$ are displayed on Figure 2. The relative accuracy of the maximum baryonic mass $M_{\mathrm{B}}$ (and of the corresponding $H_{\mathrm{max}}$) obtained for 28 grid points and for 38 points is $10^{-6}$, which is sufficient for our study taking into account all the simplifications made in the model. As a high number of grid points results in a longer computing time (up to several hours for one single configuration) without a noticeable improvement of the precision, most of the results presented in the following rely on calculations made for 28 grid points.

3.2 FlatStar, LORENE, RNS - comparison

3.2.1 Static and rigidly rotating strange stars

Finally, we compared the FlatStar code with other publicly available codes allowing the numerical calculation of the structure of static and rigidly rotating stars and with which were made the first fully relativistic calculations for rigidly rotating strange quark stars described by the MIT Bag model (Gourgoulhon et al., 1999; Stergioulas et al., 1999). The first one, RotSeq, is part of the numerical library LORENE [Langage Objet pour la RElativité NumériquE, http://www.lorene.obspm.fr, Gourgoulhon et al. (2016)], developed by a numerical relativity team from Paris Observatory in Meudon and based on spectral methods (Bonazzola & Marck, 1986). The second code, RNS, was written by Nikolaos Stergioulas and constructs models of rigidly rotating relativistic compact stars following the Komatsu et al. (1989) method, in which the field equations are converted into integral equations using appropriate Green’s functions (Stergioulas & Friedman, 1995).

We started with static (non-rotating) strange star configurations, searching for the one with maximum mass. All three codes were found to be in very good agreement, the maxima being associated to the same value of the central density $\rho_{\mathrm{c}}=2.059\cdot\mathrm{10^{15}g\ cm^{-3}}$ (central enthalpy $H_{\mathrm{c}}=0.4514$). To be more specific, the discrepancies between the values of the calculated maximum gravitational mass $M_{\mathrm{max}}^{\mathrm{stat}}$ were only of 0.005% and 0.2% with LORENE and with RNS, respectively (see Table 1).

In the case of rigid rotation, getting a precise value of the maximum mass was more difficult, mainly due to the lower accuracy of the RNS code. In order to get a meaningful comparison of all codes, we consequently adopted the strategy already used in Stergioulas et al. (1999) (Table 2 in the article): we looked for the configuration with maximum mass using LORENE, and then performed calculations for stars with the corresponding $\rho_{c}$ using the RNS and FlatStar codes. The differences between the outputs of the codes were larger than for static stars but still in good agreement: depending on the physical quantity, the values differ of about 0.1% - 1.0% with LORENE and of 0.02% - 1.1% with RNS, as can be seen in the second part of Table 1.

3.2.2 Differentially rotating stars

The LORENE library also includes a tool to perform calculations of differentially rotating compact stars, RotDiff. In order to compare results of RotDiff and FlatStar we calculated sequences of differentially rotating strange stars for two fixed values of the maximum density (maximum enthalpy) $\rho_{\mathrm{max}}=0.82\cdot 10^{15}$ and $1.70\cdot 10^{15}\,\mathrm{g/cm^{3}}$ ($H_{\mathrm{max}}=0.2$ and $0.4$) as well as of the degree of differential rotation $\tilde{A}=0.1$ and $0.2$. As will be explained in more detail in Section 4, each of those sequences starts from a non-rotating spherical configuration ($\Omega_{c}=0$ with $r_{\mathrm{ratio}}=r_{p}/r_{e}=1$) and terminates at the Keplerian limit (with $\Omega_{c}\neq 0$ and $0<r_{\mathrm{ratio}}<1$), when matter from the equator can no longer be retained by the gravitational attraction if the star spins faster. The results of these simulations are summarized on Figures 3 and 4 which display the evolution of the gravitational mass $M$ as a function of $r_{\mathrm{ratio}}$ along such sequences.

As can be observed, it was not possible with LORENE to build complete sequences of configurations for the chosen values of the parameters. While the agreement with FlatStar is very good for slowly rotating stars, the code RotDiff cannot deal with very flattened stars, as was already observed for proto-neutron stars in Villain et al. (2004). More precisely, the lowest value of $r_{\mathrm{ratio}}$ that could be reached for strange stars with LORENE was $\sim 0.52$. On the other hand, FlatStar could follow sequences up to their physical endpoint, as is illustrated by Fig. 5 which depicts the meridian cross section of a SQS at the Keplerian limit with $\tilde{A}=0.2$ and $\rho_{\mathrm{max}}=0.82\cdot 10^{15}\mathrm{g/cm^{3}}$ ($H_{\mathrm{max}}=0.2$). The detailed results of the comparison between FlatStar and LORENE for differentially rotating strange stars are gathered in Tabs. 2 (for $H_{\mathrm{max}}=0.2$) and 3 (for $H_{\mathrm{max}}=0.4$), in which two values of $r_{\mathrm{ratio}}$ ($0.9$ and $0.7$) are taken into account.

4 Differentially rotating Strange Quark Stars

4.1 Building sequences

Our first step in building sequences of configurations was to calculate the structure of spherically symmetric stars with given maximum densities $\rho_{\mathrm{max}}$. Then, we obtained a complete sequence by keeping a value of $\rho_{\mathrm{max}}$, fixing the degree of differential rotation $\tilde{A}$ (with, in the case of rigid rotation, $\tilde{A}=0$), and using a third parameter to spin up the star, e.g. by increasing the central angular velocity $\Omega_{c}$ or by decreasing the ratio of polar and equatorial radii $r_{\mathrm{ratio}}$. We remind the reader that the algorithmic of our code makes it possible to fix or vary any physical parameter (see Ansorg et al., 2009; Gondek-Rosińska et al., 2017, for more detail).

As was explained in Ansorg et al. (2009) for homogeneous stars or $N=1$ polytropes and by Studzińska et al. (2016) for other adiabatic indices such sequences end in different way depending on the degree of differential rotation $\tilde{A}$. Indeed, when the latter is sufficiently low, the sequences with fixed $\rho_{\mathrm{max}}$ terminate at the mass-shedding limit with a finite value of $r_{\mathrm{ratio}}$, as sequences of rigidly rotating stars do. As stated above, the fact that this limit is reached shows itself by the appearance of cusps along the equator of the star, giving it a shape similar to the one that can be seen on Fig. 5. Such sequences are referred to as of type A. However, when $\tilde{A}$ is large enough, a sequence with a fixed value of $\rho_{\mathrm{max}}$ can be extended down to $r_{\mathrm{ratio}}=0$, situation which indicates the entrance into the toroidal regime and at which we arbitrarily ended the sequence (as in our previous studies). In that case, the shape of the star looks like the one illustrated by Fig. 7 and the sequence is said to be of type C. More precisely, for each value of $\rho_{\mathrm{max}}$, there is a critical value $\tilde{A}_{\mathrm{crit}}(\rho_{\mathrm{max}})$ of $\tilde{A}$ such that a sequence with fixed $\rho_{\mathrm{max}}$ and $\tilde{A}$ starting from a non-rotating star ends at the mass-shedding limit (with a finite non-zero $r_{\mathrm{ratio}}$) if $\tilde{A}<\tilde{A}_{\mathrm{crit}}(\rho_{\mathrm{max}})$ but enters into the toroidal regime (when $r_{\mathrm{ratio}}=0$) if $\tilde{A}>\tilde{A}_{\mathrm{crit}}(\rho_{\mathrm{max}})$. This phenomenon, discovered for homogeneous stars and $N=1$ polytropes in Ansorg et al. (2009) as well as for polytropes with other indices in Studzińska et al. (2016), was confirmed for SQS as we shall see below.

Before discussing in more detail the critical value of $\tilde{A}$ which describes the transition between the A and C types of sequences, one shall remind the reader of another conclusion of Ansorg et al. (2009) and Studzińska et al. (2016): the structure of the solution space of differentially rotating relativistic stars is even richer since, for $\tilde{A}\sim\tilde{A}_{\mathrm{crit}}(\rho_{\mathrm{max}})$, it also contains some sequences of configurations that, for fixed $\rho_{\mathrm{max}}$ and $\tilde{A}$, do not have a non-rotating limit. Those configurations, called of types B and D, which coexist with type A and type C respectively. They are much more complex to obtain with a numerical code, which is why they were ignored before Ansorg et al. (2009). For strange matter as well, we were able to numerically calculate their properties, and their typical shapes are illustrated by meridian cross sections on Figs. 6 (type B) and 8 (type D).

By looking at those representative examples, it can be deduced that stars of those types are characterized by small values of $r_{\mathrm{ratio}}$. This feature and the global structure of the solution space for a given value of $\rho_{\mathrm{max}}$ are nevertheless better captured when one represents the sequences of configurations with fixed $\rho_{\mathrm{max}}$ and $\tilde{A}$ in a $(\tilde{\beta},r_{\mathrm{ratio}})$ plane, as illustrated by Fig. 9 for $\rho_{\mathrm{max}}=0.82\cdot 10^{15}\,\mathrm{g/cm^{3}}$ ($H_{\mathrm{max}}=0.2$). On this figure, one easily sees that the plane is divided, by the curve made of configurations with the critical degree of differential rotation $\tilde{A}_{\mathrm{crit}}(\rho_{\mathrm{max}})$, in four domains corresponding to the four types of sequences. As was already stated in Ansorg et al. (2009), Gondek-Rosińska et al. (2017) and Studzińska et al. (2016), type B sequences start at the mass-shedding limit $(\tilde{\beta}=0)$ and end at ($\tilde{\beta}=1,r_{\mathrm{ratio}}=0$), while type D sequences both start and end at the mass-shedding limit.

Another fact that can be observed on Fig. 9 is that stars of types B and D exist in reduced ranges of degree of differential rotation (especially type D), type B occurring for $\tilde{A}_{B}\leq\tilde{A}\leq\tilde{A}_{\mathrm{crit}}$ (always when type A exists as well but type C doesn’t), and type D for $\tilde{A}_{D}\geq\tilde{A}\geq\tilde{A}_{\mathrm{crit}}$ (always when type C exists as well but type A no longer does). This property is illustrated by Fig. 10. On the latter are plotted, for strange quark stars modeled by the MIT Bag model EOS (1), and as functions of $\rho_{\mathrm{max}}$ (top axis), or equivalently of $H_{\mathrm{max}}$ (bottom axis):

• •

as a black curve, the critical value $\tilde{A}_{\mathrm{crit}}$;

• •

as a blue curve, the value $\tilde{A}_{D}$ such that for $\tilde{A}_{D}\geq\tilde{A}\geq\tilde{A}_{\mathrm{crit}}$ types D and C coexist, but for $\tilde{A}>\tilde{A}_{D}$ only type C is left;

• •

as a red curve, the value $\tilde{A}_{B}$ such that for $\tilde{A}_{B}\leq\tilde{A}\leq\tilde{A}_{\mathrm{crit}}$ types B and A coexist, but for $\tilde{A}<\tilde{A}_{B}$ only type A is left.

Comparing with the results of Ansorg et al. (2009) and of Studzińska et al. (2016), it can be noticed on this figure that when $\rho_{\mathrm{max}}$ is sufficiently small ($\rho_{\mathrm{max}}\sim 0.5\cdot 10^{15}\mathrm{g}/\mathrm{cm^{3}}$, equivalent to $H_{\mathrm{max}}\sim 0.05$, and below), the critical value attached to the MIT Bag model $\tilde{A}_{\mathrm{crit}}$ is very close to the one for an incompressible fluid, but when $\rho_{\mathrm{max}}$ increases, the critical value for the compressible quark fluid becomes larger and reaches values close to those obtained for (polytropic) neutron stars [see Figs. 4, 6 of Ansorg et al. (2009) and of Studzińska et al. (2016)]. This behaviour is in agreement with the fact that for the MIT Bag model EOS (1) the adiabatic index $\gamma=d\ln P/d\ln n$, where $n$ is the baryon number density, diverges when $\rho\to\rho_{0}$ (low density limit) and has $a+1$ as a limit for $\rho\to\infty$ (see Haensel et al., 2007). Finally, the specific values of $\tilde{A}_{\mathrm{crit}}$ also lead to the conclusion that for most $\rho_{\mathrm{max}}$, $\tilde{A}_{\mathrm{crit}}(\rho_{\mathrm{max}})$ is around $0.3$ so that sequences will generally be of type A (or B) for $\tilde{A}<0.3$ and of type C for $\tilde{A}>0.3$ (D being quite rare).

When looking for the maximum mass of a sequence of configurations with fixed $\rho_{\mathrm{max}}$ and $\tilde{A}$, as we shall do in the following, having in mind the value of $\tilde{A}_{\mathrm{crit}}(\rho_{\mathrm{max}})$ and the type of the sequence which is followed is quite important. Indeed, while for a type A sequence (that ends at the mass-shedding limit and exists for $\tilde{A}<\tilde{A}_{\mathrm{crit}}$) the maximum mass on the sequence is a well-defined concept, it is not for those of types B and C which enter into the toroidal regime. As already explained, this difference is due to the fact that terminating type B and C sequences at $r_{\mathrm{ratio}}=0$ is an arbitrary choice and the mass could a priori become even larger for other non-simply-connected configurations with identical values of $\rho_{\mathrm{max}}$ and $\tilde{A}$. In fact, as will be discussed in the next Subsection, the maximum mass of type C sequences was always found for $r_{\mathrm{ratio}}=0$, like for polytropic EOS Gondek-Rosińska et al. (2017) and Studzińska et al. (2016) showing that it is more a supremum than an extremum.

4.2 Maximum mass of strange quark stars

As stated earlier, for a given EOS, a configuration of a differentially rotating star with the Komatsu et al. (1989) rotation law is prescribed by the values of three parameters. One usually sets the maximum density $\rho_{\mathrm{max}}$, which can be different from the central one $\rho_{c}$, the degree of differential rotation $\tilde{A}$ and a third quantity that describes how fast the star rotates (e.g. $\Omega_{\mathrm{c}}$ the central angular velocity, $r_{\mathrm{ratio}}$ the ratio between the radii, or $J$ the angular momentum). As a consequence, the first step in the usual way of proceeding to determine the maximum mass of differentially rotating stars consists in building sequences with fixed $\tilde{A}$ and $\rho_{\mathrm{max}}$, and using them to calculate, for various chosen values of $\tilde{A}$, the highest mass as a function of $\rho_{\mathrm{max}}$. Finally, the maximum of those functions is the maximum mass as a function of $\tilde{A}$, taking into account all $\rho_{\mathrm{max}}$ and all rates of rotation. The only subtlety of the described procedure is that, as explained earlier, the highly non-linear nature of the equations makes the solution space quite complicated with several types of configurations possibly coexisting, even for fixed $\tilde{A}$, $\rho_{\mathrm{max}}$ and $r_{\mathrm{ratio}}$. This property is especially crucial to have in mind when one calculates configurations of stars close to the critical curve $\tilde{A}=\tilde{A}_{\mathrm{crit}}(\rho_{\mathrm{max}})$ visible as a bold line in Fig. 9, situation in which an insufficiently robust numerical code can jump from one type of sequences to another. In the following, taking advantage of the high precision of the FlatStar code, we determine maximum masses for sequences with fixed $\tilde{A}$ and with a given type.

Applying the preceding process, we represented on Figure 11 the baryonic mass $M_{B}$ as a function of $r_{\mathrm{ratio}}$ for three complete evolutionary sequences with the same fixed degree of differential rotation $\tilde{A}=0.2$ but different values of the maximum density (enthalpy), $\rho_{\mathrm{max}}=0.82,1.70,3.64\cdot 10^{15}\mathrm{g/cm^{3}}$ (equivalently $H_{\mathrm{max}}=0.2,0.4,0.6$). For such $\tilde{A}$, Fig. 10 shows that the sequences are of type A and end at the mass-shedding limit. As can be seen on Fig. 11, the highest values of the mass are obtained for configurations (marked with a cross) close to but not at the mass-shedding limit, as was already observed for polytropes in Gondek-Rosińska et al. (2017). Notice that the star represented on Fig. 5 lies at the end (at the mass-shedding limit) of the sequence plotted with a black line on Fig 11.

On Fig. 13, the baryonic mass $M_{\mathrm{B}}$ is displayed as a function of $r_{\mathrm{ratio}}$ for three complete evolutionary sequences with the same maximum density as on Fig. 11, but this time with $\tilde{A}=0.5$ which is larger than $\tilde{A}_{\mathrm{crit}}$ for those values of $\rho_{\mathrm{max}}$. Consequently, all sequences are of type C and reach the entrance into the toroidal regime at $r_{\mathrm{ratio}}=0$. As stated in the previous Section, one easily observes that on all those sequences the highest mass is at $r_{\mathrm{ratio}}=0$. Notice that the configuration from Fig. 7 lies at the end of the sequence plotted with a black line on Fig 13.

As already explained, configurations of types B and D do not have any direct connection with a non-rotating configuration, which means that they exist only for $r_{\mathrm{ratio}}<1.0$. On Fig. 12, the baryonic mass $M_{B}$ is displayed as a function of the shedding parameter $\tilde{\beta}$ for two complete sequences with the same fixed $\tilde{A}=0.2$ but different maximum density $\rho_{\mathrm{max}}=0.59$ and $0.82\cdot 10^{15}\mathrm{g/cm^{3}}$ (maximum enthalpy $H_{\mathrm{max}}=0.1$ and $0.2$). Both are of type B, as can be realized by the fact that their extremities are at the mass-shedding limit ($\tilde{\beta}=0$) and at the entrance into the toroidal regime ($\tilde{\beta}=1$ and $r_{\mathrm{ratio}}=0$). We observed that sequences of this type have their highest mass at the Keplerian limit. For instance, the configuration depicted on Fig. 6 is the one with the highest mass of the sequence plotted as a black line on Fig. 12.

Figure 14 shows curves equivalent to those of Fig. 12, but for two complete sequences of type D. They were calculated with the same fixed $\tilde{A}=0.3$, but with different maximum density $\rho_{\mathrm{max}}=0.59$ and $0.64\cdot 10^{15}\mathrm{g/cm^{3}}$ (maximum enthalpy $H_{\mathrm{max}}=0.1$ and $0.125$). Type D sequences start and end at the mass-shedding limit ($\tilde{\beta}=0$), and we observed, as was already the case for other EOSs, that the highest mass is found for the lowest value of $r_{\mathrm{ratio}}$. The meridian cross-section shown on Fig. 8 corresponds to the maximum mass configuration of the sequence plotted as a red line on Fig. 14.

Doing similar calculations for broad ranges of $\rho_{\mathrm{max}}$, one eventually gets curves as those gathered on Fig. 15 which displays the results for various values of $\tilde{A}$ with the convention that dotted lines are for type D sequences, dashed lines for type C, dashed-dotted lines for type B and solid lines for type A (for comparison, we plotted as well the curves associated to static spherical stars – solid black line – and to rigidly rotating stars – solid red line–, for which the highest mass is obtained at the mass-shedding limit $\Omega=\Omega_{K}$).

In fact, in order to calculate with a good accuracy the maximum masses (and other physical properties of the configurations with maximum mass), we proceeded to a more meticulous exploration of the solution space than what could suggest Figs. 11, 12, 13 and 14. For instance, in the case of type A sequences, we explored $(H_{\mathrm{max}},r_{\mathrm{ratio}})$ planes (for fixed $\tilde{A}$) with a fine mesh, looking for the configuration with maximum mass using a 2D code, while for type C sequences it was sufficient to search for it on the $r_{\mathrm{ratio}}=0$ line (and to save some computational time without decreasing the precision of the results we could even end the sequences at $r_{\mathrm{ratio}}\sim 0.01$). See the Appendix of Gondek-Rosińska et al. (2017) for more details on the method.

The main results of our study of the maximum mass of differentially rotating strange stars are gathered in Tables 4 and 5, as well as on Fig. 16 and 17. On Fig. 16, we compare the relative maximum mass increase of SQS with its values for polytropes with three different $\Gamma$ indices [results taken from Gondek-Rosińska et al. (2017) - blue lines for $\Gamma=2.0$ - and Studzińska et al. (2016) - red and green lines for $\Gamma=1.8,2.5$, respectively]. We see that for SQS, the maximum mass depends on both the degree of differential rotation and on the type of solution, as was the case for neutron stars described by polytropic equations of state. However, for strange stars, types B, C and D occur at much smaller degrees of differential rotation $\tilde{A}$ than for polytropes. Additionally,the maximum mass is an increasing function of the degree of differential rotation $\tilde{A}$ when the sequence is of type A (with a highest mass close to the mass-shedding limit), but it becomes a decreasing one once $\tilde{A}$ is sufficiently large for the sequences to be of type C (in which case the maximum mass is reached at the entrance in the toroidal regime, for $r_{\mathrm{ratio}}=0$, and is a supremum but not an actual maximum), or of types B and D (both with a maximum mass at the mass-shedding limit). Taking into account simultaneously those two behaviors, the intuitive conclusion is challenged since the largest increase of the maximum mass of differentially rotating strange stars is obtained for a low degree of differential rotation $\tilde{A}$ and not for the highest ones. More quantitatively, for all values of $\tilde{A}$ considered in this study, the largest increase was found for $\tilde{A}=0.15$, giving a configuration with a baryonic mass $M_{\mathrm{B}}=10.81\ \mathrm{M_{\odot}}$ (an increase of around 308% with respect to the non-rotating case), a gravitational mass $M=7.59\ \mathrm{M_{\odot}}$, a circumferential radius $R_{\mathrm{circ}}=30.00\ \mathrm{km}$ and a maximum density $\rho_{\mathrm{max}}=0.46\cdot 10^{15}\mathrm{g}\ \mathrm{cm}^{-3}$ (see Table 5 for more precise values and additional properties). Such an increase of the maximum mass allowed by differential rotation, for small values of $\tilde{A}$, is much larger than what was found for neutron stars described by polytropic EOSs in previous studies. If we assume that GW170817 merger have led to the quite quick formation of a black hole (Rezzolla et al., 2018), such a high value is of course questionable. Direct natural conclusions could be that quark matter does not exist, or that differential rotation is dissipated on a very short timescale. However, one should also not forget that our code does nothing more than finding equilibrium configurations: it does not allow us to say anything about their dynamical stability. As a consequence, one could still reasonably think that maybe GW170817 hosted a phase transition from nuclear to quark matter, as proposed in Drago et al. (2016); Most et al. (2019), but that it was very fast followed by the collapse of the stellar corpse, even though its mass was smaller than the highest possible one. Only a deeper analysis based on dynamical numerical simulations with realistic microphysics inputs can solve this question. In addition we should keep in mind that the picture on the fate of the post-merger object is still uncertain. For example taking into account electromagnetic observations of GW170817 the long-lived massive neutron star surrounded by a torus was used (Shibata et al., 2017) to interpret the data or short lived high mass neutron star (Margalit & Metzger, 2017; Sapountzis & Janiuk, 2019).

A first step towards realistic calculations consists in checking whether there is such a difference between quark stars and neutron stars when these latter are also modeled using realistic EOSs. A first comparison between realistic differentially rotating neutron stars and quark stars can be done using our results and those of Morrison et al. (2004) (who studied the influence of various realistic nuclear EOSs on the maximum mass of differentially rotating neutron stars using the same rotation law as us), which is done for the maximum baryonic mass in Table 6 and on Figure 17.

More precisely, in the Table we gathered, for stars with various degree of differential rotation $\tilde{A}$, the relative increase of the maximum baryonic mass with respect to the maximum static mass $(\Omega=0)$. Also, we indicate between parenthesis the corresponding value of $r_{\mathrm{ratio}}$. As can be seen, we were actually able to compare neutron and quark stars for four values of the degree of differential rotation $\tilde{A}$ (one of which is the case of rigid rotation $\tilde{A}=0$) as the work of Morrison et al. (2004) did not include $\tilde{A}=0.1$ nor $0.2$, while we did not study $\tilde{A}=1$.

Table 6 shows that for each fixed degree of differential rotation the highest increase of the maximum mass, comparing to non-rotating configurations, is reached by strange quark stars, as was already known for rigid rotation. However, in Morrison et al. (2004), there is a graph (Fig. 2) which is quite analogue to our Figure 15, and which depicts the highest mass as a function of the maximum density for fixed $\tilde{A}$. Comparing those figures, and taking also into account the results presented for polytropes in Studzińska et al. (2016), it seems reasonable to conclude that while for low degrees of differential rotation the code of Morrison et al. (2004) did allow the determination of the actual maximum mass, things are not so clear for high $\tilde{A}$. Indeed the curves they present give the impression that the authors did not get proper sequences of configurations for high degrees of differential rotation and high maximal energy densities ($\epsilon_{\mathrm{max}}$), the actual maximum mass being probably larger. This assumption is supported by the fact that the critical value of $\tilde{A}$ is around $0.5$ or $0.6$ for polytropes with adiabatic indices similar to the effective indices calculated by Morrison et al. (2004) ($\Gamma_{\textrm{eff}}\sim 2.5$ or $3$), as shown in Studzińska et al. (2016). This implies that the sequences with high $\tilde{A}$ should be of type C or D and the maximum mass is a priori found for configurations with $r_{\mathrm{ratio}}\sim 0$ that their code could maybe not calculate precisely. All this would suggest that even if the relative increase of the maximum mass is in all likelihood larger for SQSs, the quantitative difference with its value for realistic models of neutron stars is still an open question.

5 Conclusions

We presented the first study of the maximum mass of differentially rotating relativistic stars made from quark matter, using the MIT Bag model. Our results were obtained with the highly accurate spectral code FlatStar, which allows to determine the structure of strongly deformed configurations. By making internal tests and by comparing our code with others in the public domain (RNS and LORENE), we showed that its accuracy is indeed quite high and that it enables a thorough and precise exploration of the solution space for broad ranges of the degree of differential rotation $\tilde{A}$ and of the maximum density $\rho_{\mathrm{max}}$.

We built numerous sequences of configurations with fixed $\tilde{A}$ and $\rho_{\mathrm{max}}$, and observed that depending on the value of $\tilde{A}$, as was already discovered for homogeneous stars and polytropes in Ansorg et al. (2009); Gondek-Rosińska et al. (2017) and Studzińska et al. (2016), there can be several types of sequences with different properties. There is indeed a critical value $\tilde{A}_{\mathrm{crit}}(\rho_{\mathrm{max}})$ such that a sequence starting from a non-rotating configuration finishes at the Keplerian limit if $\tilde{A}<\tilde{A}_{\mathrm{crit}}$, but enters into the toroidal regime if $\tilde{A}>\tilde{A}_{\mathrm{crit}}$. Also, again as was the case for constant density stars and polytropes, we found two threshold values of $\tilde{A}$, noted $\tilde{A}_{B}$ and $\tilde{A}_{D}$, such that sequences without a non-rotating limit (said of B or D type) can appear. More precisely, we have the following restrictions of the possible coexistence of the various types:

• •

$0\leq\tilde{A}<\tilde{A}_{B}$ - type A only;

• •

$\tilde{A}_{B}<\tilde{A}<\tilde{A}_{\mathrm{crit}}$ - type A or B;

• •

$\tilde{A}_{\mathrm{crit}}<\tilde{A}<\tilde{A}_{D}$ - type C or D;

• •

$\tilde{A}_{D}<\tilde{A}$ - type C only.

We find that the maximum mass of SQS depends on both the degree of differential rotation and on the type of solution. It is an increasing function of $\tilde{A}$ for type A solutions and a decreasing one for types B, C and D. For type A sequences, which end at the mass-shedding limit, the configuration with maximum mass was found close to that limit (but not at it). For type C sequences, it was obtained for the toroidal-like configurations with a vanishing ratio between the polar and equatorial radii. Since such sequences were arbitrarily terminated at that point, the conclusion is that in that case the maximum mass is just a supremum and even higher values could be obtained for non-simply connected configurations. For type B and D, which exists only for small values of $r_{\mathrm{ratio}}<0.25$, the higher mass configurations were found at the mass-shedding limit. Similar behaviors were previously observed for polytropic equations of state (Gondek-Rosińska et al., 2017; Studzińska et al., 2016).

Gathering all results, it turned out that the largest increase of the maximum mass with respect to its value in the non-rotating case is reached for type B sequences. This conclusion is again the same as for neutron stars described by a polytropic EOS (Gondek-Rosińska et al., 2017; Studzińska et al., 2016). We found that the maximum mass of differentially rotating SQS is obtained for a relatively low degree of differential rotation, $\tilde{A}\sim 0.15$, and could even be four times higher than the maximum mass of non rotating strange quark stars described by the same EOS. This increase is much larger than for neutron stars modeled with realistic EOSs as was shown by comparing our results with those of Morrison et al. (2004). A natural extension of our work would be to use the FlatStar code for realistic nuclear EOSs in order to explore the corresponding solution space in detail and make the comparison with strange matter more precise.111Recently Espino & Paschalidis (2019) confirmed that the universal relation, discovered by Gondek-Rosińska et al. (2017), between maximum mass and the degree of differential rotation is also valid for realistic EOS. Our conclusion that for small degrees of differential rotation the maximum mass of SQS is much larger than the maximum mass for differentially rotating NS still holds taking into account the results of Espino & Paschalidis (2019). Also, another important issue, to be dealt with using dynamical simulations, is the stability of massive differentially rotating stars that could be formed during core collapse or in mergers of neutron stars binaries. Indeed, knowing the theoretical life span of such bodies is crucial to put constraints on observations that could enable to demonstrate the existence of quark matter in the Universe and also to better understand the last stages of the dynamics of compact binaries. For instance, even if some analyses of GW170817 suggest a quick collapse into a black hole (Rezzolla et al., 2018; Ruiz et al., 2018; Shibata & Kiuchi, 2017), it cannot be ruled out that it was maybe predated by a transition from nuclear to quark matter (Most et al., 2019).

This work is dedicated to the memory of our deeply regretted friend and colleague Marcus Ansorg who sadly passed away before this article could be finished. Nothing could have be done without him. We also would like to thank the anonymous referee for their relevant comments. This work was supported by National Science Center grant UMO-2015/17/N/ST9/01605; by POMOST/2012-6/11 Program of Foundation for Polish Science co-financed by the European Union within the European Regional Development Fund, by the grant of the National Science Centre UMO-2014/14/M/ST9/00707 and DEC-2013/08/M/ST9/00664, by the COST Actions MP1304, CA16104, CA16214 and by the HECOLS International Associated Laboratory programme. Calculations were performed on the PIRXGW computer cluster at University of Zielona Gora funded by the Foundation for Polish Science within the FOCUS program.