Maximum Mass Of Differentially Rotating Strange Quark Stars

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.


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 electromag-netic 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 Ruiz et al. 2018;Shibata & Kiuchi 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 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 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 max,stat (Cook et al. 1992(Cook et al. , 1994a, 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 Γ = 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.

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 (or equivalently the total mass density ρ = /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 where a and ρ 0 are constants such that c √ a is the speed of sound and ρ 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 ρ 0 = 4.2785 × 10 14 g/cm 3 .
Once given a barotropic EOS p = p(ρ), 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 ρ and an angular velocity Ω. 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 Ω 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Ã, 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; •Ã = 0 corresponds to rigid rotation, Ω = Ω c in the whole star; • forÃ = 1 the angular velocity at the equatorial surface is around half Ω c ; • the degree of differential rotation is an increasing and monotonic function ofÃ.
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 ratio = r p /r e , with 0 < r ratio ≤ 1 as we consider only stars without a hole at their center, and also to the so-called "rescaled shedding parameter"β, see the precise definition in Ansorg et al. (2009), which is such that 0 <β < 1, withβ → 1 when a star enters into the toroidal regime (r ratio = 0),β → 0 at the massshedding limit andβ = 1/2 for a non-rotating spherical star (r ratio = 1).
The last technical point worth mentioning is that the FlatStar code uses as a primary thermodynamical variable not the mass density ρ, but the corresponding dimensionless relativistic enthalpy H(ρ) defined by the differential relation dH = dp c 2 ρ(p) + p , 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, ρ max , or of the maximum enthalpy H max . Notice that with the EOS (1), one has the relation so that H = 0 corresponds to ρ = ρ 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 ρ on the same axis.

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 ρ max (equivalently of the maximum enthalpy H max ). On Figure 1, we display the results for that test in the case J = 2.16 [GM 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 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 (Ã) and to check how its value changes with an increasing number of grid points. The results for SQSs withÃ = 0.1 are displayed on Figure 2. The relative accuracy of the maximum . Geometrical convergence rate for the maximum baryonic mass MB,max (black) and the corresponding maximum enthalpy Hmax (red) for differentially rotating Strange Quark Stars with a fixed degree of differential rota-tionÃ = 0.1. The plot displays the relative accuracy of the nth spectral approximation with respect to the accuracy of the 38th one.
baryonic mass M B (and of the corresponding H 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.

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 Greens functions (Stergioulas & Friedman 1995 Figure 3. Gravitational mass M as a function of the ratio between polar and equatorial radii rratio = rp/re for sequences of differentially rotating Strange Quark Stars calculated with FlatStar (solid lines) and LORENE (dots or diamonds), for a fixed maximum density (maximum enthalpy) ρmax = 0.82 · 10 15 g/cm 3 (Hmax = 0.2) and two values of the degree of differential rotationÃ = 0.1 and 0.2. As can be seen, while the agreement between the two codes is very good for slow rotation, LORENE fails to converge for rapidly rotating configurations characterized by the lowest values of rratio.
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 ρ c = 2.059·10 15 g cm −3 (central enthalpy H c = 0.4514). To be more specific, the discrepancies between the values of the calculated maximum gravitational mass M stat max 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 ρ 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.

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) ρ max = 0.82 · 10 15 and 1.70 · 10 15 g/cm 3 (H max = 0.2 and 0.4) as well as of the degree of differential rotationÃ = 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 (Ω c = 0 with r ratio = r p /r e = 1) and terminates at the Keplerian limit (with Ω c = 0 and 0 < r 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 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 ratio that could be reached for strange stars with LORENE was ∼ 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Ã = 0.2 and ρ max = 0.82 · 10 15 g/cm 3 (H 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 max = 0.2) and 3 (for H max = 0.4), in which two values of r ratio (0.9 and 0.7) are taken into account.  Gourgoulhon et al. (1999); RNS - Stergioulas et al. (1999)]. For all configurations, the Table displays the gravitational mass M , the baryonic mass MB, the angular velocity Ω, the central mass density ρc, the central enthalpy Hc, the circumferential radius Rcirc, the coordinate equatorial radius re, the ratio between the polar and equatorial radii rratio = rp/re, the angular momentum J and the ratio between the kinetic and potential energies T /|W |. The upper part of the

Building sequences
Our first step in building sequences of configurations was to calculate the structure of spherically symmetric stars with given maximum densities ρ max . Then, we obtained a complete sequence by keeping a value of ρ max , fixing the degree of differential rotationÃ (with, in the case of rigid rotation,Ã = 0), and using a third parameter to spin up the star, e.g. by increasing the central angular velocity Ω c or by decreasing the ratio of polar and equatorial radii r 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 ro-tationÃ. Indeed, when the latter is sufficiently low, the sequences with fixed ρ max terminate at the massshedding limit with a finite value of r 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Ã is large enough, a sequence with a fixed value of ρ max can be extended down to r 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 ρ max , there is a critical valueÃ crit (ρ max ) ofÃ such that a sequence with fixed ρ max andÃ starting from a non-rotating star ends at the mass-shedding limit (with a finite non-zero r ratio ) ifÃ <Ã crit (ρ max ) but enters into the toroidal regime (when r ratio = 0) ifÃ >Ã crit (ρ 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 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) andStudzińska et al. (2016): the structure of the solution space of differentially rotating relativistic stars is even richer since, forÃ ∼Ã crit (ρ max ), it also contains some sequences of configurations that, for fixed ρ max andÃ, 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 ratio . This feature and the global structure of the solution space for a given value of ρ max are nevertheless better captured when one represents the sequences of configurations with fixed ρ max andÃ in a (β, r ratio ) plane, as illustrated by Fig. 9 for ρ max = 0.82 · 10 15 g/cm 3 (H max = 0.2). On this fig-  ure, one easily sees that the plane is divided, by the curve made of configurations with the critical degree of differential rotationÃ crit (ρ max ), in four domains corresponding to the four types of sequences. As was already stated in Ansorg et al. (2009 and Studzińska et al. (2016), type B sequences start at the mass-shedding limit (β = 0) and end at (β = 1, r 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Ã B ≤Ã ≤Ã crit (always when type A exists as well but type C doesn't), and type D forÃ D ≥Ã ≥Ã 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 ρ max (top axis), or equivalently of H max (bottom axis): • as a black curve, the critical valueÃ crit ; • as a blue curve, the valueÃ D such that forÃ D ≥ A ≥Ã crit types D and C coexist, but forÃ >Ã D only type C is left; • as a red curve, the valueÃ B such that forÃ B ≤ A ≤Ã crit types B and A coexist, but forÃ <Ã 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 ρ max is sufficiently small (ρ max ∼ 0.5 · 10 15 g/cm 3 , equivalent to H max ∼ 0.05, and below), the critical value attached to the MIT Bag model A crit is very close to the one for an incompressible fluid, but when ρ 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 γ = d ln P/d ln n, where n is the baryon number density, diverges when ρ → ρ 0 (low density limit) and has a + 1 as a limit for ρ → ∞ (see Haensel et al. 2007). Finally, the specific values ofÃ crit also lead to the conclusion that for most ρ max ,Ã crit (ρ max ) is around 0.3 so that sequences will generally be of type A (or B) for A < 0.3 and of type C forÃ > 0.3 (D being quite rare).
When looking for the maximum mass of a sequence of configurations with fixed ρ max andÃ, as we shall do in the following, having in mind the value ofÃ crit (ρ 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Ã <Ã 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 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 ρ max andÃ. In fact, as will be discussed in the next Subsection, the maximum mass of type C sequences was always found for r 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.

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 ρ max , which can be different from the central one ρ c , the degree of differential rotationÃ and a third quantity that describes how fast the star rotates (e.g. Ω c the central angular velocity, r 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Ã and ρ max , and using them to calculate, for various chosen values ofÃ, the highest mass as a function of ρ max . Finally, the maximum of those functions is the maximum mass as a function ofÃ, taking into account all ρ max and all rates of rotation. The only subtlety of the described procedure is that, as explained earlier, the highly nonlinear nature of the equations makes the solution space quite complicated with several types of configurations possibly coexisting, even for fixedÃ, ρ max and r ratio . This property is especially crucial to have in mind when one calculates configurations of stars close to the critical curveÃ =Ã crit (ρ 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Ã and with a given type.
Applying the preceding process, we represented on Figure 11 the baryonic mass M B as a function of r ratio for three complete evolutionary sequences with the same fixed degree of differential rotationÃ = 0.2 but different values of the maximum density (enthalpy), ρ max = 0.82, 1.70, 3.64·10 15 g/cm 3 (equivalently H max = 0.2, 0.4, 0.6). For suchÃ, 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 2), with a fixed degree of differential rotationÃ = 0.2. For the chosen ρmax, AB ≤ 0.2 ≤Ãcrit, so that the sequences can be of type B, which they are. Consequently they both start at the massshedding limit (β = 0) and end at the entrance into the toroidal regime, rratio → 0 andβ = 1. For all such sequences, the highest mass is found at the mass-shedding limit which also corresponds to the highest value of rratio.
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 B is displayed as a function of r ratio for three complete evolutionary se- , with a fixed degree of differential rotationÃ = 0.5. Since for the chosen ρmax,Ãcrit(ρmax) < 0.5, those sequences, which start from non-rotating configurations, are of type C and consequently end at the entrance in the toroidal regime, rratio → 0. For all of them, the highest mass is found at this specific point.
quences with the same maximum density as on Fig. 11, but this time withÃ = 0.5 which is larger thanÃ crit for those values of ρ max . Consequently, all sequences are of type C and reach the entrance into the toroidal regime at r ratio = 0. As stated in the previous Section, one easily observes that on all those sequences the highest mass is at r 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 ratio < 1.0. On Fig. 12, the baryonic mass M B is displayed as a function of the shedding parameterβ for two complete sequences with the same fixedÃ = 0.2 but different maximum density ρ max = 0.59 and 0.82 · 10 15 g/cm 3 (maximum enthalpy H 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 (β = 0) and at the entrance into the toroidal regime (β = 1 and r 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Ã = 0.3, but with different maximum density ρ max = 0.59 and 0.64 · 10 15 g/cm 3 (maximum enthalpy H max = 0.1 and 0.125). Type ρ max = 0.59 ⋅ 10 15 g/cm 3 ρ max = 0.64 ⋅ 10 15 g/cm 3 Figure 14. Baryonic mass MB as a function ofβ for two values of the maximum density (enthalpy) ρmax = 0.59 and 0.64 · 10 15 g/cm 3 (Hmax = 0.1 and 0.125), with a fixed degree of differential rotationÃ = 0.3. Since for the chosen ρmax, AD ≥ 0.3 ≥Ãcrit, the sequences can be of type D, which they are. Consequently they both start and end at the massshedding limit (β = 0). In practice, we observed that for stars of this type, the highest mass was always found at the mass-shedding limit with the smallest value of rratio.
D sequences start and end at the mass-shedding limit (β = 0), and we observed, as was already the case for other EOSs, that the highest mass is found for the lowest value of r 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 ρ max , one eventually gets curves as those gathered on Fig. 15 which displays the results for various values ofÃ 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 Ω = Ω 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 max , r ratio ) planes (for fixedÃ) 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 ratio = 0 line (and to save some computational time without decreasing the precision of the results we could even end the sequences at r ratio ∼ 0.01).  Figure 15. Highest baryonic mass MB as a function of the maximum enthalpy Hmax (bottom axis) or the maximum density ρmax (top axis) for fixed values of the degree of differential rotationÃ. The black solid line represents nonrotating stars, the red solid line rigid rotation (for which the highest mass is at the mass-shedding limit), while other lines correspond to differentially rotating stars. Solid lines stand for type A sequences for which the highest mass is reached close to the mass-shedding limit, dashed-line are for type C sequences on which the highest mass is obtained at rratio = 0, dashed-dotted lines are for type B sequences for which the highest mass is reached at the mass-shedding limit and finally dotted lines are for type D sequences whose highest mass is found at the mass-shedding limit for their smallest value of rratio. For all values ofÃ, configurations with maximum mass are marked with a black cross. -red and green lines for Γ = 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Ã than for polytropes. Additionally,the maximum mass is an increasing function of the degree of differential rotationÃ when the sequence is of type A (with a highest mass close to the mass-shedding limit), but it becomes a decreasing one onceÃ 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 ratio = 0, and is a supremum Table 4. Properties of configurations with maximum mass of non-rotating (Ã = 0.0, Ω = 0.0), rigidly rotating (Ã = 0.0 -Keplerian limit) and differentially rotating (Ã = 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 MB, the gravitational mass M , the central and equatorial angular velocities Ωc and Ωeq, the corresponding frequencies fc and feq, the central and maximal densities ρc and ρmax, the corresponding enthalpies Hc and Hmax, the circumferential radius Rcirc, the coordinate equatorial radius re, the ratio between the polar and equatorial radii rratio = rp/re, the angular momentum J, and the ratio between the kinetic and potential energies T /|W |. A = 0.0Ã = 0.0Ã = 0.1Ã = 0.2Ã = 0.3Ã = 0.5Ã = 0.7 Ω = 0.  Table are given the baryonic mass MB, the gravitational mass M , the central and equatorial angular velocities Ωc and Ωeq, the corresponding frequencies fc and feq, the central and maximal densities ρc and ρmax, the corresponding enthalpies Hc and Hmax, the circumferential radius Rcirc, the coordinate equatorial radius re, the ratio between the polar and equatorial radii rratio = rp/re, the angular momentum J, and the ratio between the kinetic and potential energies T /|W |. 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Ã and not for the highest ones. More quantitatively, for all values ofÃ considered in this study, the largest increase was found forÃ = 0.15, giving a configuration with a baryonic mass M B = 10.81 M (an increase of around 308% with respect to the non-rotating case), a gravitational mass M = 7.59 M , a circumferential radius R circ = 30.00 km and a maximum density ρ max = 0.46 · 10 15 g 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Ã, 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 , 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 ) 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 rota-tion 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Ã, the relative increase of the maximum baryonic mass with respect to the maximum static mass (Ω = 0). Also, we indicate between parenthesis the corresponding value of r ratio . As can be seen, we were actually able to compare neutron and quark stars for four values of the degree of differential rotationÃ (one of which is the case of rigid rotatioñ A = 0) as the work of Morrison et al. (2004) did not includeÃ = 0.1 nor 0.2, while we did not studyÃ = 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Ã. 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Ã. 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 ( max ), the actual maximum mass being probably larger. This assumption is supported by the fact that the critical value ofÃ is around 0.5 or 0.6 for polytropes with adiabatic indices similar to the effective indices calculated by Morrison et al. (2004) (Γ eff ∼ 2.5 or 3), as shown in Studzińska et al. (2016). This implies that the sequences with highÃ should be of type C or D and the maximum mass is a priori found for configurations with r ratio ∼ 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.

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 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 rratio 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). a A -Reid soft core (Pandharipande 1971).
f APR -A18 + δυ + UIX* (Akmal et al. 1998). Figure 16. Maximum baryonic mass of differentially rotating strange quark stars and neutron stars described by a polytropic EOS normalized by the maximum mass of a nonrotating configuration with the same EOS as a function of the degree of differential rotationÃ. Black lines correspond to our results for four types of differentially rotating SQS, other lines corresponds to differentially rotating polytropes: blue -Γ = 2.0 taken from Gondek-Rosińska et al. (2017), red -Γ = 1.8 and green -Γ = 2.5 respectively taken from Studzińska et al. (2016).
indeed quite high and that it enables a thorough and precise exploration of the solution space for broad ranges Figure 17. Maximum baryonic mass of differentially rotating strange stars normalized by the maximum mass in the non-rotating case as a function of the degree of differential rotationÃ. The lines without points correspond to our results (Fig.15), while lines with symbols are for differentially rotating neutron stars described by different realistic equations of states taken from Morrison et al. (2004).
of the degree of differential rotationÃ and of the maximum density ρ max . We built numerous sequences of configurations with fixedÃ and ρ max , and observed that depending on the value ofÃ, 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Ã crit (ρ max ) such that a sequence starting from a non-rotating configuration finishes at the Keplerian limit ifÃ <Ã crit , but enters into the toroidal regime ifÃ >Ã crit . Also, again as was the case for constant density stars and polytropes, we found two threshold values ofÃ, notedÃ B andÃ 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 ≤Ã <Ã B -type A only; •Ã B <Ã <Ã crit -type A or B; •Ã crit <Ã <Ã D -type C or D; •Ã D <Ã -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Ã 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 ratio < 0.25, the higher mass configurations were found at the mass-shedding limit. Similar behaviors were previously observed for polytropic equations of state 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 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,Ã ∼ 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. 1 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 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).