Neutron star cooling and GW170817 constraint within quark-meson coupling models

  • In the present work, we used five different versions of the quark-meson coupling (QMC) model to compute astrophysical quantities related to the GW170817 event and neutron star cooling process. Two of the models are based on the original bag potential structure and three versions consider a harmonic oscillator potential to confine the quarks. The bag-like models also incorporate the pasta phase used to describe the inner crust of neutron stars. Within the simple method studied in the present work, we show that the pasta phase does not play a significant role. Moreover, the QMC model that satisfies the GW170817 constraints with the lowest slope of the symmetry energy exhibits a cooling profile compatible with observational data.
  • 加载中
  • [1] B. P. Abbott et al., The Astrophysical Journal 848, L12 (2017) doi: 10.3847/2041-8213/aa91c9
    [2] B. P. Abbott et al., The lsc-virgo white paper on instrument science (2016-2017 edition), Technical ReportNo. T1400226, LIGO and Virgo Scientific Collaborations, LIGO Internal Document No. T1400226 (2016).
    [3] J. Aasi et al., Classical and Quantum Gravity 32, 074001 (2015) doi: 10.1088/0264-9381/32/7/074001
    [4] F. Acernese et al., Classical and Quantum Gravity 32, 024001 (2014)
    [5] Y. Aso et al. (The KAGRA Collaboration), Phys. Rev. D 88, 043007 (2013) doi: 10.1103/PhysRevD.88.043007
    [6] R. Essick, S. Vitale, and M. Evans, Phys. Rev. D 96, 084004 (2017) doi: 10.1103/PhysRevD.96.084004
    [7] K. Chamberlain and N. Yunes, Phys. Rev. D 96, 084039 (2017) doi: 10.1103/PhysRevD.96.084039
    [8] M. Punturo, M. Abernathy, F. Acernese et al., Classical and Quantum Gravity 27, 084007 (2010) doi: 10.1088/0264-9381/27/8/084007
    [9] K. D. Kokkotas and G. Schafer, Mon. Not. R. Astron Soc. 275, 301 (1995) doi: 10.1093/mnras/275.2.301
    [10] S. Tsuruta and A. G. W. Cameron, Nature 207, 364 (1965) doi: 10.1038/207364a0
    [11] O. V. Maxwell, The Astrophysical Journal 231, 201 (1979) doi: 10.1086/157181
    [12] D. Page and S. Reddy, Annu. Rev. Nucl. Part. Sci. 56, 327 (2006) doi: 10.1146/annurev.nucl.56.080805.140600
    [13] F. Weber, R. Negreiros, P. Rosenfield, and M. Stejner, Prog. Part. Nucl. Phys. 59, 94 (2007) doi: 10.1016/j.ppnp.2006.12.008
    [14] R. Negreiros, V. A. Dexheimer, and S. Schramm, Phys. Rev. C 82, 035803 (2010) doi: 10.1103/PhysRevC.82.035803
    [15] R. Negreiros, L. Tolos, M. Centelles, A. Ramos, and V. Dexheimer, The Astrophysical Journal 863, 104 (2018) doi: 10.3847/1538-4357/aad049
    [16] D. N. Aguilera, J. A. Pons, and J. A. Miralles, Astron. Astrophys. 486, 255 (2008) doi: 10.1051/0004-6361:20078786
    [17] J. A. Pons, J. A. Miralles, and U. Geppert, Astron. Astrophys. 496, 207 (2009) doi: 10.1051/0004-6361:200811229
    [18] B. Niebergal, R. Ouyed, R. Negreiros, and F. Weber, Phys. Rev. D 81, 1 (2010)
    [19] B. Franzon, R. Negreiros, and S. Schramm, Phys. Rev. D 96, 123005 (2017) doi: 10.1103/PhysRevD.96.123005
    [20] J. Horvath, O. Benvenuto, and H. Vucetich, Phys. Rev. D 44, 3797 (1991) doi: 10.1103/PhysRevD.44.3797
    [21] D. Blaschke, T. Klahn, and D. N. Voskresensky, The Astrophysical Journal 533, 406 (2000) doi: 10.1086/308664
    [22] I. Shovkovy and P. Ellis, Phys. Rev. C 66, 015802 (2002)
    [23] H. Grigorian, D. Blaschke, and D. Voskresensky, Phys. Rev. C 71, 045801 (2005)
    [24] M. Alford, M. Braby, M. Paris, and S. Reddy, The Astrophysical Journal 629, 969 (2005) doi: 10.1086/430902
    [25] K. P. Levenfish and D. G. Yakovlev, Astronomy Letters 20, 43 (1994)
    [26] C. Schaab and D. Voskresensky, Astron. Astrophys. 321, 591 (1997)
    [27] M. Alford, P. Jotwani, C. Kouvaris, J. Kundu, and K. Rajagopal, Phys. Rev. D 71, 114011 (2005) doi: 10.1103/PhysRevD.71.114011
    [28] D. Page, J. M. Lattimer, M. Prakash, and A. W. Steiner, The Astrophysical Journal 707, 1131 (2009) doi: 10.1088/0004-637X/707/2/1131
    [29] R. Negreiros, S. Schramm, and F. Weber, Phys. Rev. D 85, 104019 (2012) doi: 10.1103/PhysRevD.85.104019
    [30] R. Negreiros, S. Schramm, and F. Weber, Phys. Lett. B 718, 1176 (2013) doi: 10.1016/j.physletb.2012.12.046
    [31] R. Negreiros, S. Schramm, and F. Weber, Astron. Astrophys. 603, A44 (2017) doi: 10.1051/0004-6361/201730435
    [32] M. E. Gusakov, A. D. Kaminker, D. G. Yakovlev, and O. Y. Gnedin, Mon. Not. R. Astron Soc. 363, 555 (2005) doi: 10.1111/j.1365-2966.2005.09459.x
    [33] F. Weber, Progress in Particle and Nuclear Physics 54, 193 (2005) doi: 10.1016/j.ppnp.2004.07.001
    [34] T. Frederico, B. V. Carlson, R. A. Rego, and M. S. Hussein, J. Phys. G: Nucl. Part. Phys. 15, 297 (1989) doi: 10.1088/0954-3899/15/3/007
    [35] E. F. Batista, B. V. Carlson, and T. Frederico, Nucl. Phys. A 697, 469 (2002) doi: 10.1016/S0375-9474(01)01250-7
    [36] P. A. M. Guichon, Phys. Lett. B 200, 235 (1988) doi: 10.1016/0370-2693(88)90762-9
    [37] K. Saito and A. W. Thomas, Phys. Lett. B 327, 9 (1994) doi: 10.1016/0370-2693(94)91520-2
    [38] P. K. Panda, A. Mishra, J. M. Eisenberg, and W. Greiner, Phys. Rev. C 56, 3134 (1997)
    [39] A. M. Santos, C. c. Providência, and P. K. Panda, Phys. Rev. C 79, 045805 (2009) doi: 10.1103/PhysRevC.79.045805
    [40] P. Guichon, J. Stone, and A. Thomas, Progress in Particle and Nuclear Physics 100, 262 (2018) doi: 10.1016/j.ppnp.2018.01.008
    [41] J. R. Stone, P. Guichon, H. Matevosyan, and A. Thomas, Nuclear Physics A 792, 341 (2007) doi: 10.1016/j.nuclphysa.2007.05.011
    [42] M. Dutra, O. Lourenço, S. S. Avancini, B. V. Carlson, A. Delfino, D. P. Menezes, C. Providência, S. Typel, and J. R. Stone, Phys. Rev. C 90, 055203 (2014) doi: 10.1103/PhysRevC.90.055203
    [43] M. Oertel, M. Hempel, T. Klähn, and S. Typel, Rev. Mod. Phys. 89, 015007 (2017) doi: 10.1103/RevModPhys.89.015007
    [44] C. Providência, S. S. Avancini, R. Cavagnoli, S. Chiacchiera, C. Ducoin, F. Grill, J. Margueron, D. P. Menezes, A. Rabhi, and I. Vidaña, The European Physical Journal A 50, 44 (2014) doi: 10.1140/epja/i2014-14044-7
    [45] P. K. Panda, A. M. S. Santos, D. P. Menezes, and C. Providência, Phys. Rev. C 85, 055802 (2012)
    [46] R. Cavagnoli, D. P. Menezes, and C. m. c. Providência, Phys. Rev. C 84, 065810 (2011) doi: 10.1103/PhysRevC.84.065810
    [47] G. Grams, A. M. Santos, P. K. Panda, C. Providência, and D. P. Menezes, Phys. Rev. C 95, 055807 (2017) doi: 10.1103/PhysRevC.95.055807
    [48] T. Maruyama, T. Tatsumi, D. N. Voskresensky, T. Tanigawa, and S. Chiba, Phys. Rev. C 72, 015802 (2005) doi: 10.1103/PhysRevC.72.015802
    [49] S. S. Avancini, D. P. Menezes, M. D. Alloy, J. R. Marinelli, M. M. W. Moraes, and C. Providência, Phys. Rev. C 78, 015802 (2008)
    [50] H. Pais and J. R. Stone, Physical Review Letters, 109 (2012) doi: 10.1103/physrevlett.109.151101
    [51] A. S. Schneider, D. K. Berry, M. E. Caplan, C. J. Horowitz, and Z. Lin, Phys. Rev. C 93, 065806 (2016) doi: 10.1103/PhysRevC.93.065806
    [52] A. S. Schneider, M. E. Caplan, D. K. Berry, and C. J. Horowitz, Phys. Rev. C 98, 055801 (2018) doi: 10.1103/PhysRevC.98.055801
    [53] C. C. Barros, D. P. Menezes, and F. Gulminelli, Physical Review C, 101 (2020) doi: 10.1103/physrevc.101.035211
    [54] B. Schuetrumpf, G. Martínez-Pinedo, M. Afibuzzaman, and H. M. Aktulga, “A survey of nuclear pasta in the intermediate density regime i: Shapes and energies,” (2019), arXiv: 1906.08155[nucl-th].
    [55] N. Barik, R. N. Mishra, D. K. Mohanty, P. K. Panda, and T. Frederico, Phys. Rev. C 88, 015206 (2013) doi: 10.1103/PhysRevC.88.015206
    [56] R. N. Mishra, H. S. Sahoo, P. K. Panda, N. Barik, and T. Frederico, Phys. Rev. C 92, 045203 (2015)
    [57] R. N. Mishra, H. S. Sahoo, P. K. Panda, N. Barik, and T. Frederico, Phys. Rev. C 94, 035805 (2016) doi: 10.1103/PhysRevC.94.035805
    [58] J. R. Stone, N. J. Stone, and S. A. Moszkowski, Phys. Rev. C 89, 044316 (2014) doi: 10.1103/PhysRevC.89.044316
    [59] F. J. Fattoyev, J. Carvajal, W. G. Newton, and B.-A. Li, Phys. Rev. C 87, 015806 (2013) doi: 10.1103/PhysRevC.87.015806
    [60] T. Malik, N. Alam, M. Fortin, C. Providência, B. K. Agrawal, T. K. Jha, B. Kumar, and S. K. Patra, Phys. Rev. C 98, 035804 (2018) doi: 10.1103/PhysRevC.98.035804
    [61] N. Hornick, L. Tolos, A. Zacchi, J.-E. Christian, and J. Schaffner-Bielich, Phys. Rev. C 98, 065804 (2018) doi: 10.1103/PhysRevC.98.065804
    [62] B. Kumar, S. K. Biswal, and S. K. Patra, Phys. Rev. C 95, 015801 (2017) doi: 10.1103/PhysRevC.95.015801
    [63] T. Damour and A. Nagar, Phys. Rev. D 80, 084035 (2009) doi: 10.1103/PhysRevD.80.084035
    [64] T. Binnington and E. Poisson, Phys. Rev. D 80, 084018 (2009) doi: 10.1103/PhysRevD.80.084018
    [65] E. E. Flanagan and T. Hinderer, Phys. Rev. D 77, 021502 (2008) doi: 10.1103/PhysRevD.77.021502
    [66] J. Piekarewicz and F. J. Fattoyev, Phys. Rev. C 99, 045802 (2019) doi: 10.1103/PhysRevC.99.045802
    [67] G. Baym, C. Pethick, and P. Sutherland, Astrophys. J. 170, 299 (1971) doi: 10.1086/151216
    [68] A. F. Fantina, S. D. Ridder, N. Chamel, and F. Gulminelli, Astronomy & Astrophysics 633, A149 (2020)
    [69] D. G. Yakovlev, A. D. Kaminker, O. Y. Gnedin, and P. Haensel, Phys. Rep. 354, 1 (2000)
    [70] D. Page, J. M. Lattimer, M. Prakash, and A. W. Steiner, The Astrophysical Journal Suppl. Ser. 155, 623 (2004) doi: 10.1086/424844
    [71] S. Beloin, S. Han, A. W. Steiner, and D. Page, Phys. Rev. C 97, 015804 (2018) doi: 10.1103/PhysRevC.97.015804
    [72] S. Safi-Harb and H. S. Kumar, The Astrophysical Journal 684, 532 (2008) doi: 10.1086/590359
    [73] V. E. Zavlin, J. Trumper, and G. G. Pavlov, The Astrophysical Journal 525, 959 (1999) doi: 10.1086/307919
    [74] V. E. Zavlin, The Astrophysical Journal 665, L143 (2007) doi: 10.1086/521300
    [75] V. E. Zavlin, in Neutron Stars and Pulsars (Springer Berlin Heidelberg, Berlin, Heidelberg, 2009) pp. 181–211.
    [76] G. G. Pavlov, V. E. Zavlin, D. Sanwal, V. Burwitz, and G. P. Garmire, The Astrophysical Journal 552, L129 (2001) doi: 10.1086/320342
    [77] G. G. Pavlov, V. E. Zavlin, D. Sanwal, and J. Trümper, The Astrophysical Journal 569, L95 (2002) doi: 10.1086/340640
    [78] S. Mereghetti, G. F. Bignami, and P. A. Caraveo, The Astrophysical Journal 464, 842 (1996) doi: 10.1086/177370
    [79] E. V. Gotthelf, J. P. Halpern, and R. Dodson, The Astrophysical Journal 567, L125 (2002) doi: 10.1086/340109
    [80] K. E. McGowan, J. A. Kennea, S. Zane, F. A. Cordova, M. Cropper, C. Ho, T. Sasseen, and W. T. Vestrand, The Astrophysical Journal 591, 380 (2003) doi: 10.1086/375332
    [81] K. E. McGowan, S. Zane, M. Cropper, J. A. Kennea, F. A. Cordova, C. Ho, T. Sasseen, and W. T. Vestrand, The Astrophysical Journal 600, 343 (2004) doi: 10.1086/379787
    [82] K. E. McGowan, S. Zane, M. Cropper, W. T. Vestrand, and C. Ho, The Astrophysical Journal 639, 377 (2006) doi: 10.1086/497327
    [83] D. Klochkov, V. Suleimanov, G. Pühlhofer, D. G. Yakovlev, A. Santangelo, and K. Werner, Astron. Astrophys. 573, A53 (2015) doi: 10.1051/0004-6361/201424683
    [84] A. Possenti, S. Mereghetti, and M. Colpi, Astron. Astrophys. 313, 565 (1996)
    [85] J. P. Halpern and F. Y. Wang, The Astrophysical Journal 477, 905 (1997) doi: 10.1086/303743
    [86] J. A. Pons, F. M. Walter, J. M. Lattimer, M. Prakash, R. Neuhauser, and P. An, The Astrophysical Journal 564, 981 (2002) doi: 10.1086/324296
    [87] V. Burwitz, F. Haberl, R. Neuhäuser, P. Predehl, J. Trümper, and V. E. Zavlin, Astron. Astrophys. 399, 1109 (2003) doi: 10.1051/0004-6361:20021747
    [88] D. L. Kaplan, M. H. van Kerkwijk, H. L. Marshall, B. A. Jacoby, S. R. Kulkarni, and D. A. Frail, The Astrophysical Journal 590, 1008 (2003) doi: 10.1086/375052
    [89] W. C. Ho, K. G. Elshamouty, C. O. Heinke, and A. Y. Potekhin, Phys. Rev. C 91, 1 (2015)
    [90] D. Page, U. Geppert, and F. Weber, Nucl. Phys. A 777, 497 (2006) doi: 10.1016/j.nuclphysa.2005.09.019
    [91] F. Weber, Pulsars as astrophysical laboratories for nuclear and particle physics /F. Weber. Bristol, U.K. : Institute of Physics, c1999. QB 464 W42 1999. DA (Institute of Physics, Bristol, U.K., 1999).
    [92] C. Schaab, F. Weber, M. K. Weigel, and N. K. Glendenning, Nucl. Phys. A 605, 531 (1996) doi: 10.1016/0375-9474(96)00164-9
    [93] A. Sedrakian, Phys. Rev. D 93, 065044 (2016) doi: 10.1103/PhysRevD.93.065044
    [94] C. Kouvaris, Phys. Rev. D 77, 023006 (2008) doi: 10.1103/PhysRevD.77.023006
    [95] M. Prakash, M. Prakash, J. M. Lattimer, and C. J. Pethick, The Astrophysical Journal 390, 77 (1992) doi: 10.1086/186376
    [96] B. Friman and O. Maxwell, The Astrophysical Journal 232, 541 (1979) doi: 10.1086/157313
    [97] A. B. Migdal, E. E. Saperstein, M. A. Troitsky, and D. N. Voskresensky, Physics Reports 192, 179 (1990) doi: 10.1016/0370-1573(90)90132-L
    [98] C. Schaab, D. Voskresensky, A. D. Sedrakian, F. Weber, and M. K. Weigel, Astronomy and Astrophysics 321, 591 (1997), arXiv:astro-ph/9605188[astro-ph
    [99] D. Blaschke, H. Grigorian, and D. N. Voskresensky, Astronomy and Astrophysics 424, 979 (2004), arXiv:astroph/0403170[astro-ph doi: 10.1051/0004-6361:20040404
    [100] D. Blaschke, H. Grigorian, and D. N. Voskresensky, Physical Review C 88, 065805 (2013), arXiv:1308.4093 [nucl-th doi: 10.1103/PhysRevC.88.065805
    [101] E. Flowers and N. Itoh, The Astrophysical Journal 206, 218 (1976) doi: 10.1086/154375
    [102] E. Flowers and N. Itoh, Astrophys. J. 230, 847 (1979) doi: 10.1086/157145
    [103] E. Flowers and N. Itoh, Astrophys. J. 250, 750 (1981) doi: 10.1086/159423
    [104] E. Gudmundsson, The Astrophysical Journal 259, L19 (1982) doi: 10.1086/183840
    [105] E. H. Gudmundsson, C. J. Pethick, and R. I. Epstein, The Astrophysical Journal 272, 286 (1983) doi: 10.1086/161292
    [106] B. P. Abbott et al. (The LIGO Scientific Collaboration and the Virgo Collaboration), Phys. Rev. Lett. 121, 161101 (2018) doi: 10.1103/PhysRevLett.121.161101
    [107] E. R. Most, L. R. Weih, L. Rezzolla, and J. SchaffnerBielich, Phys. Rev. Lett. 120, 261103 (2018) doi: 10.1103/PhysRevLett.120.261103
    [108] C. A. Raithel, F. Özel, and D. Psaltis, The Astrophysical Journal 857, L23 (2018) doi: 10.3847/2041-8213/aabcbf
    [109] E. Annala, T. Gorda, A. Kurkela, and A. Vuorinen, Phys. Rev. Lett. 120, 172703 (2018) doi: 10.1103/PhysRevLett.120.172703
    [110] O. Lourenço, M. Dutra, C. H. Lenzi, C. V. Flores, and D. P. Menezes, Phys. Rev. C 99, 045202 (2019) doi: 10.1103/PhysRevC.99.045202
    [111] B. P. Abbott et al. (LIGO Scientific Collaboration and Virgo Collaboration), Phys. Rev. Lett. 119, 161101 (2017) doi: 10.1103/PhysRevLett.119.161101
    [112] S. M. de Carvalho, R. Negreiros, M. Orsaria, G. A. Contrera, F. Weber, and W. Spinella, Phys. Rev. C 92, 035810 (2015) doi: 10.1103/PhysRevC.92.035810
    [113] V. Dexheimer, R. Negreiros, and S. Schramm, Phys. Rev. C 92, 012801 (2015) doi: 10.1103/PhysRevC.92.012801
    [114] A. R. Raduta, A. Sedrakian, and F. Weber, Mon. Not. R. Astron Soc. 475, 4347 (2017)
    [115] B. P. Abbott et al., Phys. Rev. X 9, 011001 (2019)
    [116] V. Dexheimer, R. Negreiros, and S. Schramm, Phys. Rev. C 92, 012801 (2015) doi: 10.1103/PhysRevC.92.012801
    [117] P. S. Shternin, Phys. Rev. D 98, 063015 (2018) doi: 10.1103/PhysRevD.98.063015
    [118] A. Bauswein, N.-U. F. Bastian, D. B. Blaschke, K. Chatziioannou, J. A. Clark, T. Fischer, and M. Oertel, Phys. Rev. Lett. 122, 061102 (2019) doi: 10.1103/PhysRevLett.122.061102
    [119] M. Lucca and L. Sagunski, “The lifetime of binary neutron star merger remnants, ” (arXiv: 1909.08631), arXiv: 1909.08631[astro-ph.HE].
    [120] S. Y. Lau, P. T. Leung, and L.-M. Lin, Phys. Rev. D 95, 101302 (2017) doi: 10.1103/PhysRevD.95.101302
    [121] E. Annala, T. Gorda, A. Kurkela, J. Nättilä, and A. Vuorinen, Nature Physics, (2020) doi: 10.1038/s41567-020-0914-9
    [122] NICER and visited on October 2, https://heasarc.gsfc.nasa.gov/docs/nicer/ (2019).
    [123] J.-E. Christian, A. Zacchi, and J. Schaffner-Bielich, Phys. Rev. D 99, 023009 (2019) doi: 10.1103/PhysRevD.99.023009
  • 加载中

Figures(14) / Tables(2)

Get Citation
Odilon Lourenço, César H. Lenzi, Mariana Dutra, Tobias Frederico, M. Bhuyan, Rodrigo Negreiros, César V. Flores, Guilherme Grams and Débora P. Menezes. Neutron star cooling and GW170817 constraint within quark-meson coupling models[J]. Chinese Physics C.
Odilon Lourenço, César H. Lenzi, Mariana Dutra, Tobias Frederico, M. Bhuyan, Rodrigo Negreiros, César V. Flores, Guilherme Grams and Débora P. Menezes. Neutron star cooling and GW170817 constraint within quark-meson coupling models[J]. Chinese Physics C. shu
Milestone
Article Metric

Article Views(10)
PDF Downloads(2)
Cited by(0)
Policy on re-use
To reuse of subscription content published by CPC, the users need to request permission from CPC, unless the content was published under an Open Access license which automatically permits that type of reuse.
通讯作者: 陈斌, bchen63@163.com
  • 1. 

    沈阳化工大学材料科学与工程学院 沈阳 110142

  1. 本站搜索
  2. 百度学术搜索
  3. 万方数据库搜索
  4. CNKI搜索

Email This Article

Title:
Email:

Neutron star cooling and GW170817 constraint within quark-meson coupling models

  • Departamento de Física, Instituto Tecnológico de Aeronáutica, DCTA, 12228-900, São José dos Campos, SP, Brazil
  • Department of Physics, Faculty of Science, University of Malaya, Kuala Lumpur 50603, Malaysia Institute of Research and Development, Duy Tan University, Da Nang 550000, Vietnam
  • Instituto de Física, Universidade Federal Fluminense 20420, Niterói, RJ, Brazil
  • Centro de Ciências Exatas, Naturais e Tecnológicas, UEMASUL, Rua Godofredo Viana 1300, Centro CEP: 65901- 480, Imperatriz, Maranhão, Brazil
  • Departamento de Física, Universidade Federal de Santa Catarina, Florianópolis, SC, CP 476, CEP 88.040-900, Brazil Univ. Lyon, Univ. Claude Bernard Lyon 1, CNRS/IN2P3, IP2I Lyon, UMR 5822, F-69622, Villeurbanne, France
  • Departamento de Física, Universidade Federal de Santa Catarina, Florianópolis, SC, CP 476, CEP 88.040-900, Brazil

Abstract: In the present work, we used five different versions of the quark-meson coupling (QMC) model to compute astrophysical quantities related to the GW170817 event and neutron star cooling process. Two of the models are based on the original bag potential structure and three versions consider a harmonic oscillator potential to confine the quarks. The bag-like models also incorporate the pasta phase used to describe the inner crust of neutron stars. Within the simple method studied in the present work, we show that the pasta phase does not play a significant role. Moreover, the QMC model that satisfies the GW170817 constraints with the lowest slope of the symmetry energy exhibits a cooling profile compatible with observational data.

    HTML

    I.   INTRODUCTION
    • The observation of the binary neutron star system GW170817 by the LIGO-VIRGO scientific collaboration and also in the X-ray, ultraviolet, optical, infrared, and radio bands gave rise to the new era of multimessenger astronomy. All these joint observations supported the GW170817 neutron star merger event in NGC4993. The observation of the electromagnetic counterpart of GW170817 by the Fermi Gamma-ray observatory corroborated that binary neutron star mergers are associated with short gamma-ray bursts (sGRB) and kilonova emissions probably powered by the radioactive decay of r-process nuclei synthesized in the ejecta (Abbott et al. (2017) [1]). Therefore multi-messenger observations are an excellent tool to extract the information of compact stars under extreme conditions. In the next years, further upgrades of dubbed A+ and LIGO Voyager (Abbott et al. (2016) [2] and the collaboration of KAGRA and LIGO India (Aasi et al. (2014) [3]; Acernese et al. (2014) [4]; Aso et al. (2013) [5]) will certainly report new results. Besides, planned third-generation observatories such as Einstein Telescope (ET) and Cosmic Explorer (CE) may bring us the detection of binary neutron stars at cosmological distances (Essick et al. (2017) [6]; Chamberlain & Yunes (2017) [ 7]; Punturo et al. (2010) [8].

      The observation of binary systems in different channels is of outstanding importance to establish reliable constraints on neutron star physics. As a matter of fact neutron stars have been, for many decades, objects of intense astrophysical research. Electromagnetic observations were used to establish some limits on the mass and radius of neutron stars. In light of the recent detection, the gravitational wave channel has become a new way to infer neutron star (NS) properties. The dimensionless tidal deformability/polarizability (TP) and/or Love number are related to the induced deformation that a neutron star undergoes by the influence of the tidal field of its neutron star companion in the binary system (Kokkotas et al. (1995) [9]). That influence is expected to be detected in the low frequency range of the inspiral stage, where the effect is a small correction in the waveform phase. As each different neutron star composition has a characteristic response to the tidal field, the TP can be used to discriminate between different equations of state (EOS). In addition, thermal evolution studies are a complementary way of confronting the properties of EOS describing neutron stars with existing data. The investigation of the cooling of compact stars has been proved as a promising method for exploring the properties of neutron stars, as it is strongly connected to both micro and macroscopic realms (Tsuruta & Cameron (1965) [ 10]; Maxwell (1979) [11]; Page & Reddy (2006) [ 12]; Weber et al. (2007) [13]; Negreiros et al. (2010) [14]; (2018) [15]). There is a wealth of literature on the cooling of neutron stars and related phenomena such as magnetic field (Aguilera et al. (2008) [16]; Pons et al. (2009) [17]; Niebergal et al. (2010) [18]; Franzon et al. (2017) [19]), deconfined quark matter (Horvath et al. (1991) [20]; Blaschke et al. (2000) [21]; Shovkovy & Ellis (2002) [ 22]; Grigorian et al. (2005) [23]; Alford et al. (2005) [24]), superfluidity (Levenfish & Yakovlev (1994) [ 25]; Schaab & Voskresensky (1997) [ 26]; Alford et al. (2005b) [27]; Page et al. (2009) [28]), rotation (Negreiros et al. (2012) [29]; (2013) [30]; (2017) [31]), among others (Negreiros et al. (2010) [14]; Alford et al. (2010) [27]; Gusakov et al. (2005) [32]; Weber (2005) [33]).

      The observation of binary systems in different channels is of outstanding importance to establish reliable constraints on neutron star physics. As a matter of fact neutron stars have been, for many decades, objects of intense astrophysical research. Electromagnetic observations were used to establish some limits on the mass and radius of neutron stars. In light of the recent detection, the gravitational wave channel has become a new way to infer neutron star (NS) properties. The dimensionless tidal deformability/polarizability (TP) and/or Love number are related to the induced deformation that a neutron star undergoes by the influence of the tidal field of its neutron star companion in the binary system (Kokkotas et al. (1995) [9]). That influence is expected to be detected in the low frequency range of the inspiral stage, where the effect is a small correction in the waveform phase. As each different neutron star composition has a characteristic response to the tidal field, the TP can be used to discriminate between different equations of state (EOS). In addition, thermal evolution studies are a complementary way of confronting the properties of EOS describing neutron stars with existing data. The investigation of the cooling of compact stars has been proved as a promising method for exploring the properties of neutron stars, as it is strongly connected to both micro and macroscopic realms (Tsuruta & Cameron (1965) [ 10]; Maxwell (1979) [11]; Page & Reddy (2006) [ 12]; Weber et al. (2007) [13]; Negreiros et al. (2010) [14]; (2018) [15]). There is a wealth of literature on the cooling of neutron stars and related phenomena such as magnetic field (Aguilera et al. (2008) [16]; Pons et al. (2009) [17]; Niebergal et al. (2010) [18]; Franzon et al. (2017) [19]), deconfined quark matter (Horvath et al. (1991) [20]; Blaschke et al. (2000) [21]; Shovkovy & Ellis (2002) [ 22]; Grigorian et al. (2005) [23]; Alford et al. (2005) [24]), superfluidity (Levenfish & Yakovlev (1994) [ 25]; Schaab & Voskresensky (1997) [ 26]; Alford et al. (2005b) [27]; Page et al. (2009) [28]), rotation (Negreiros et al. (2012) [29]; (2013) [30]; (2017) [31]), among others (Negreiros et al. (2010) [14]; Alford et al. (2010) [27]; Gusakov et al. (2005) [32]; Weber (2005) [33]).

      In the present work, we study five different versions of the quark-meson coupling (QMC) model equations of state chosen to describe neutron star (NS) matter. The first two are respectively the original QMC model with an underlying bag structure and its counterpart with the inclusion of an interaction between the meson $ \omega $ and $ \rho $ fields. In both cases, the pasta phase is also considered in the description of the NS inner crust. The other three versions are modified QMC (MQMC), where the parameters are adjusted so that the constituent quarks are confined to a flavor-independent harmonic oscillator potential (Frederico et al. (1989) [34]; Batista et al. (2002) [35]). These models are confronted with recent astrophysical constraints, including the ones obtained from the detection of GW170817. We also perform cooling simulations for all models and compare the results with observed data. This study will allow us (in conjunction with the aforementioned studies) a better evaluation of the quality of the underlying microscopic models adopted.

      Here we analyze the macroscopic properties obtained from a relativistic model by taking the quark degrees of freedom into account, which is the modified version of the commonly used in quantum hadrodynamics (also referred to as Walecka-type) models. Furthermore, the inclusion of other internal phase transitions and their consequences on the thermal evolution of neutron stars will be investigated shortly.

      The formalism used for the EOS is shown in Sec. II, the necessary equations to estimate the quantities related to the gravitational wave (GW) constraints are given in Sec. III and for the cooling process in Sec. IV. We present and discuss our results in Sec. V and in the last section, we make our brief remarks.

    II.   QMC AND QMC $ \omega\rho $ MODELS
    • In this section, we present the formalism of the QMC model, its counterpart with the inclusion of the $ \omega\rho $ interaction and in the sequel, the MQMC model.

      In the QMC model, the nucleon in a nuclear medium is treated as a collection of non-overlapping static spherical MIT bags. The interactions arise through the coupling of quarks to various meson fields. Here we have considered the quark-mesons interactions for the scalar ( $ \sigma $ ), vector ( $ \omega $ ), and iso-vector vector ( $ \rho $ ) fields and those are treated as classical fields in the mean-field approximation (MFA) (Guichon (1988) [36]; Saito & Thomas (1994) [ 37]; Panda et. al. (1997)[38]). The quark field $ \psi_{q_{N}} $ satisfies the equation of motion within the bag, and is given by,

      $\begin{aligned}[b] \Bigg\{i \gamma . \partial -& (m_q-g_\sigma^q \sigma)-g_\omega^q\, \omega\,\gamma^0 \\ +&\frac{1}{2} g^q_\rho \tau_z \rho_{03}\gamma^0 \Bigg\} \,\psi_{q_{N}}(x) = 0. \quad q = u,d \end{aligned} $

      (1)

      Here $ m_q $ stands for the current quark mass, and $ \tau_z $ is the 3 $ ^{rd} $ component of the iso-spin. The constants $ g_\sigma^q $ , $ g_\omega^q $ , and $ g_\rho^q $ are the quark-meson coupling for $ \sigma- $ , $ \omega- $ , and $ \rho- $ field, respectively. The normalized ground state for a quark in the bag is given by,

      $ \begin{aligned}[b] \psi_{q_{N}}({\bf{r}}, t) =& {\cal N}_{q_{N}} \exp \left(-i\epsilon_{q_{N}} t/R_N \right) \\ &\times \left( \begin{array}{c} j_{0_{N}}\left(x_{q_{N}} r/R_N\right)\\ i\beta_{q_{N}} \vec{\sigma} \cdot \hat r j_{1_{N}}\left(x_{q_{N}} r/R_N\right) \end{array}\right) \frac{\chi_q}{\sqrt{4\pi}}, \end{aligned} $

      (2)

      where, $ j_0 $ and $ j_1 $ are spherical Bessel functions and $ \chi_q $ is the quark spinor. The quantity $ \epsilon_{q_{N}} $ is given by:

      $ \epsilon_{q_{N}} = \Omega_{q_{N}}+R_N\left(g_\omega^q\, \omega+ \frac{1}{2} g^q_\rho \tau_z \rho_{03} \right), $

      (3)

      where,

      $ \beta_{q_{N}} = \sqrt{\frac{\Omega_{q_{N}}-R_N\, m_q^*}{\Omega_{q_{N}}\, +R_N\, m_q^* }}\ , $

      (4)

      with the normalization factor reading,

      $ {\cal N}_{q_{N}}^{-2} = 2R_N^3 j_0^2(x_q)\left[\Omega_q(\Omega_q-1) + R_N m_q^*/2 \right] \Big/ x_q^2, $

      (5)

      and $ \Omega_{q_{N}}\equiv \sqrt{x_{q_{N}}^2+(R_N\, m_q^*)^2} $ , $ m_q^* = m_q-g_\sigma^q\, \sigma $ , $ R_N $ is the bag radius of nucleon N. The eigenvalue for the nucleon N ( $ x_{q_{N}} $ ) is determined by

      $ j_{0_{N}}(x_{q_{N}}) = \beta_{q_{N}}\, j_{1_{N}}(x_{q_{N}}), $

      (6)

      which is the boundary condition at the bag surface. The value of $ x_{q_{N}} $ depends on the flavor through the mass m $ _q $ , and satisfies the boundary condition given in Eq. (6) for a converged solution.

      The energy of a static bag that describes nucleon N, which in turn, consists of three quarks in the ground state is expressed as

      $ E^{\rm bag}_N = \sum\limits_q n_q \, \frac{\Omega_{q_{N}}}{R_N}-\frac{Z_N}{R_N} +\frac{4}{3}\, \pi \, R_N^3\, B_N\ , $

      (7)

      where $ Z_N $ , and $ B_N $ are respectively the zero-point motion of the nucleon and the bag constant. The quark density is $ n_q $ . We choose hereafter the symbol n to denote a generic density, with two exceptions. The first one is when it appears as a subscript of a certain quantity. For these cases, the quantity refers to the neutron one, as for instance, in $ \mu_n $ (neutron chemical potential). The second one is concerning the appearing of n in the definitions of the direct Urca, modified Urca and Bremssthralung processes (see Sec. IV). For these cases, n refers to the neutron.

      The set of parameters used in the present work is determined by enforcing the stability of the nucleon (the "bag"), as for instance in (Guichon (1988) [36], Santos et al. (2009) [39]) for identical proton and neutron masses. The effective mass of a nucleon bag at rest is then given by $ M_N^* = E_N^{\rm bag} $ and the equilibrium condition for the bag is obtained by minimizing the effective mass, $ M_N^* $ with respect to the bag radius, i.e.,

      $ \frac{d\, M_N^*}{d\, R_N^*} = 0, $

      (8)

      with N denoting protons or neutrons. The unknown parameters ZN = 4.0050668 and B $ ^{1/4}_N $ = 210.85 MeV are determined by fixing the bag radius RN = 0.6 fm and the bare nucleon mass MN = 939 MeV. The desired values for the binding energy and saturation density shown in Table I are obtained for g $ _\sigma^q $ = 5.9810, g $ _{\omega} $ = 3g $ ^q_\omega $ = 8.9817, and g $ _\rho $ = g $ ^q_\rho $ , the last one also given in Table I. The meson masses are m $ _{\sigma} $ = 550 MeV, m $ _{\omega} $ = 783 MeV and m $ _{\rho} $ = 770 MeV. At this point, some comments are in order: the QMC model was developed in such a way that it describes both finite nuclei and nuclear matter properties. A nice review can be found in Guichon et al. (2018) [40]. In the present work, we use a mean field approximation (RMF), which is the simplest variation of the model. Other parametrizations with larger values of RN are found in the literature, as in the review just mentioned, Guichon (1988) [36] and Stone et al. (2007) [41]. In Guichon (1988) [36] a simpler version of the model is used ( $ \rho $ mesons are not included) and $ R_N $ left to vary from 0.6 to 1.0 fm, the lower values (0.6 and 0.7 fm) being the ones that could better reproduce $ \sigma $ and $ \omega $ coupling values fitted from the Bonn potential. Moreover, with these low $ R_N $ values, the fraction of volume occupied by the bags is small, indicating that the overlap between the bags is negligible. On the other hand, in Stone et al. (2007) [41], the authors show that the nucleon magnetic moments cannot be reproduced with values lower than 0.8 fm. Nevertheless, it is a consensus that the equation of state, seen by plotting the binding energy in Guichon (1988) [36] and the pressure as a function of energy density in Stone et al. (2007) [41] is insensitive to the chosen $ R_N $ value.

      Quantity Model
      QMC QMC $ \omega\rho $
      $ \Lambda_v $ 0.00 0.03
      $ g_\rho $ 8.6510 9.0078
      $ M_N^*/M_N $ 0.77 0.77
      J (MeV) 34.50 30.92
      $ L_0 $ (MeV) 90.00 69.17
      $ K_0 $ (MeV) 295 295
      $ M_{max} $ (M $ _{\odot} $ ) 2.14 2.07
      $ R_{M_{max}} $ (km) 11.51 10.96
      $ R_{1.4} $ (km) 13.55 12.83

      Table 1.  Nuclear matter and stellar properties obtained with the QMC and QMC $ \omega\rho $ models. For all characterizations, one has $ B/A = -16.4 $ MeV, $ n_0 $ = 0.15 fm $ ^{-3} $ .

      The quantity $ \Lambda_v $ in Table I, and also in Eqs. (9)-(10) is the nonlinear cross coupling between the isoscalar-vector and isovector-vector mesons ( $ \rho-\omega $ ). The detailed expressions for this cross-coupling can be found in Dutra et al. (2014) [42] and references therein. The more important saturation properties of nuclear matter (i.e., compressibility $ K_0 $ , symmetry energy J, and slope of the symmetry energy $ L_0 $ ), and stellar matter observables (i.e., maximum star mass M $ _{max} $ , maximum radius R $ _{M_{max}} $ , and star radius corresponding to 1.4M $ _{\odot} $ ) for QMC and QMC $ \omega \rho $ are listed in Table I. It is worth mentioning that the parameter $ \Lambda_v $ that gives the strength of the $ \omega-\rho $ interaction plays an important role on the symmetry energy value and its related quantities (see e.g., Dutra et al. (2014) [42]; Oertel et al. (2017) [43]; Providência et al. (2014) [44]; Panda et al. (2012) [45]; Cavagnoli et al. (2011) [46]). In other words, the larger the value of $ \Lambda_v $ , the lower the value of the symmetry energy and vice-verse (Panda et al. (2012) [45]; Cavagnoli et al. (2011) [46]). In the present work $ \Lambda_v $ = 0.03, as listed in Table I and we have followed the calculations mentioned in Grams et al. (2017) [47], where an $ \omega-\rho $ interaction strength that results in a symmetry energy equal to 22 MeV at 0.1 fm $ ^{-3} $ was used, with a consequent change in the $ g_\rho $ coupling constant, whose values are shown in Table I. The new values of the symmetry energy and its slope at saturation are also given in Table I. The total energy density of the nuclear matter within the relativistic mean field (RMF) approach is given by,

      $ \begin{aligned}[b] \varepsilon =& \frac{1}{2}m^{2}_{\sigma}\sigma^2+ \frac{1}{2}m^{2}_{\omega}\omega^{2}_{0}+ \frac{1}{2}m^{2}_{\rho}\rho^{2}_{03} +3\Lambda_v g_{\omega}^{2}g_{\rho}^{2}\omega^{2}_{0}\rho^{2}_{03} \\ & +\sum\limits_{N} \frac{1}{\pi^{2}}\int^{k_{N}}_{0}k^{2}dk[k^{2}+M^{*2}_{N}]^{1/2}. \end{aligned} $

      (9)

      The pressure is also expressed as,

      $ \begin{aligned}[b] p =& -\frac{1}{2}m^{2}_{\sigma}\sigma^2+ \frac{1}{2}m^{2}_{\omega}\omega^{2}_{0}+ \frac{1}{2}m^{2}_{\rho}\rho^{2}_{03} +\Lambda_v g_{\omega}^{2}g_{\rho}^{2}\omega^{2}_{0}\rho^{2}_{03} \\ & +\sum\limits_{N} \frac{1}{\pi^{2}}\int^{k_{N}}_{0}k^{4}dk/[k^{2}+M^{*2}_{N}]^{1/2}. \end{aligned} $

      (10)

      The vector mean field $ \omega_0 $ and $ \rho_{03} $ are determined through

      $ \omega_0 = \frac{g_\omega (n_p +n_n)}{m^{*^{2}}_{\omega}}, \; \rho_{03} = \frac{g_\rho (n_p -n_n)}{2 m^{*^{2}}_{\rho}}, $

      (11)

      where

      $ n_B = n_p+n_n = \sum\limits_N \frac{2 k_{N}^3}{3 \pi ^2} $

      (12)

      is the baryonic density ( $ n_p $ and $ n_n $ are the proton and neutron densities, respectively), and the effective masses of the meson fields are

      $ m^{{*}^{2}}_{\omega} = m^{2}_{\omega}+2\Lambda_v g_{\omega}^{2}g_{\rho}^{2}\rho^2_{03} $

      (13)

      and

      $ m^{{*}^{2}}_{\rho} = m^{2}_{\rho}+2\Lambda_v g_{\omega}^{2}g_{\rho}^{2}\omega^2_{0}. $

      (14)

      Finally, the mean field $ \sigma $ is fixed by imposing that

      $ \frac{\partial \varepsilon}{\partial \sigma} = 0. $

      (15)

      In stellar matter, the $ \beta $ -equilibrium, i.e., $ \mu_p = \mu_n-\mu_e, \quad \mu_e = \mu_\mu $ , and the charge neutrality $ n_p = n_e + n_\mu $ conditions are enforced. These conditions imply that a free gas of leptons (electrons and muons) need to be included in the energy density and pressure of the system.

      Pasta phases: The pasta phases are constructed within the QMC and QMC $ \omega\rho $ models by using the coexisting phases method (Maruyama et al. (2005) [48]; Avancini et al. (2008) [49]). Although extensively used in quantum hadrodynamical models, it was obtained for the first time for the QMC model in Grams et al. (2017) [47] and the main steps and equations are discussed next. For a given total density $ n_B $ the pasta structures are built with different geometrical forms, usually called sphere (bubble), cylinder (tube), and slab, in three, two, and one dimension, respectively. Following Gibbs conditions, this imposes that both phases have the same pressure and chemical potentials for proton and neutron. For stellar matter, the following equations must be solved simultaneously,

      $ P^I = P^{II}, $

      (16)

      $ \mu_p^I = \mu_p^{II}, $

      (17)

      $ \mu_n^I = \mu_n^{II}, $

      (18)

      $ f(n^{I}_p-n^{I}_e)+(1-f)(n^{II}_p-n^{II}_e) = 0. $

      (19)

      where I ( $ II $ ) represents the high (low) density region, $ n_p $ is the global proton density and f is the volume fraction of the phase I, that reads

      $ f = \frac{n_B-n_B^{II}}{n_B^I-n_B^{II}}. $

      (20)

      The total hadronic matter energy reads:

      $ \varepsilon_{matter} = f\varepsilon^I+(1-f)\varepsilon^{II}+\varepsilon_e. $

      (21)

      The total energy can be obtained by adding the surface and Coulomb terms to the matter energy in Eq. (21),

      $ \varepsilon = \varepsilon_{matter} +\varepsilon_{surf}+\varepsilon_{Coul}. $

      (22)

      Minimizing $ \varepsilon_{surf}+\varepsilon_{Coul} $ with respect to the size of the droplet/bubble, cylinder/tube or slabs, we obtain (Maruyama et al. (2005) [48]) $ \varepsilon_{surf} = 2 \varepsilon_{Coul} $ where

      $ \varepsilon_{Coul} = \frac{2\alpha}{4^{2/3}}(e^{2}\pi \Phi)^{1/3}\left[ \sigma^{surf} D(n^{I}_p -n^{II}_p)\right] ^{2/3}, $

      (23)

      with $ \alpha = f $ for droplets, rods and slabs, and $ \alpha = 1-f $ for tubes and bubbles. The quantity $ \Phi $ is given by

      $ \Phi = \left\{ {\begin{array}{*{20}{l}} {\left( {\dfrac{{2 - D{\alpha ^{1 - 2/D}}}}{{D - 2}} + \alpha } \right)\dfrac{1}{{D + 2}},& D = 1,3}\\ {\dfrac{{\alpha - 1 - \ln \alpha }}{{D + 2}}},& D = 2 \end{array}} \right. $

      (24)

      where, $ \sigma^{surf} $ is the surface tension, which measures the energy per area necessary to create a planar interface between the two regions. The surface tension is a crucial quantity in the pasta calculation and in this work we follow the recipe used in Grams et al. (2017) [47], where it is calculated by using an adapted geometric approach that relies on the fitting of the energy per unit area necessary to create a planar interface between two phases. Within this approach, the pasta phase obtained for matter in $ \beta $ - equilibrium is constituted only by droplets and at zero temperature ends below 0.07 fm $ ^{-3} $ , as seen in Fig. 6 of Grams et al. (2017) [47]). Hence, the pasta phase is only present at the low-density regions of the neutron stars, and in this region, muons are not present, although they are present in the EOS that describes the homogeneous region.

      Figure 6.  Surface temperature as a function of age for stars within the QMC $\omega \rho$ model. Each curve represents the cooling of stars with the indicated mass.

      There are many properties of the pasta phase that remain to be better explored in astrophysical scenarios. For example, its transport properties, as the electron conductivity and its consequences on the shear and bulk viscosities related to the instability of rotating compact stars. However, to the best of our knowledge, the present analysis is the first one that considers a pasta-phase structure in the calculation of both tidal deformability and NS cooling. Hence, the phase coexistence method just described may be taken as the first step towards more sophisticated treatments. The present procedure gives only a naive picture of the real system, which has many coexisting geometries as the ones obtained in Skyrme-Hartree-Fock plus BCS treatment (Pais & Stone (2012) [ 50]), in molecular dynamics simulations (Schneider et al. (2016) [51] (2018) [52]) and very recently also reproduced by allowing statistical fluctuations on the pasta composition (Barros et al. (2020) [53]). Other even more exotic structures are the ones nominated triple periodic minimal surface pasta configurations (Schuetrumpf et al. (2019) [54]). These extremely sophisticated pasta configurations can only be obtained at a very expensive computational cost and, hence are out of the scope of the present work. The more viable approach for future already expensive calculations that cannot use yet other heavy codes as input is the one that takes into account coexisting pasta shapes based on the equilibrium distribution of the different geometries mentioned in Barros et al. (2020) [ 53].

    • A.   MQMC model

    • In this subsection, we present the improved version of the modified quark-meson coupling (MQMC) model given in Barik et al. (2013) [55], Mishra et al. (2015) [56] and (2016) [57]. The Lagrangian density is given by,

      $ \begin{aligned}[b] {\cal{L}}_{\rm{ MQMC}} =& \overline{\psi}_q[i\gamma^\mu\partial_\mu - m_q - U(r)]\psi_q + g_\sigma^q\sigma\overline{\psi}_q\psi_q \\ & - g_\omega^q\overline{\psi}_q\gamma^\mu\omega_\mu\psi_q - \frac{g_\rho^q}{2}\overline{\psi}_q\gamma^\mu\vec{\rho}_\mu\vec{\tau}\psi_q + \frac{1}{2} \partial^\mu \sigma \partial_\mu \sigma \\ & - \frac{1}{2} m^2_\sigma\sigma^2 -\frac{1}{4}F^{\mu\nu}F_{\mu\nu} + \frac{1}{2}m^2_\omega\omega_\mu\omega^\mu \\ &- \frac{1}{4}\vec{B}^{\mu\nu}\vec{B}_{\mu\nu} + \frac{1}{2}m^2_\rho\vec{\rho}_\mu\vec{\rho}^\mu. \end{aligned} $

      (25)

      All the quantities/parameters are defined as usual in the literature and some of them were already defined earlier in Sec. II. The field tensors $ F_{\mu\nu} $ and $ \vec{B}_{\mu\nu} $ are given by $ F_{\mu\nu} = \partial_\nu\omega_\mu-\partial_\mu\omega_\nu $ and $ \vec{B}_{\mu\nu} = \partial_\nu\vec{\rho}_\mu-\partial_\mu\vec{\rho}_\nu $ . It is worth mentioning that the MQMC model contains the same kind of nucleon interactions as in the original QMC model from Guichon (1988) [36]. However, the interaction between quarks inside the nucleon is taken into account via a confining harmonic oscillator potential given by $ U(r) = \frac{1}{2}(1+\gamma^0)V(r) $ , with $ V(r) = ar^2+V_0 $ , instead of the bag-like treatment. The potential intensity and depth are related, respectively, to the constants a and $ V_0 $ .

      Using Euler-Lagrange equations to the Lagrangian density in Eq. (25), one can write the Dirac equation for the quarks as,

      $ \begin{aligned}[b] \Bigg\{\vec{\alpha}\cdot\vec{k} &+ \gamma^0\left[m_q-V_\sigma+\frac{V(r)}{2}\right] \\ &+ \frac{V(r)}{2} + V_\omega + \frac{\tau_zV_\rho}{2}\Bigg\}\psi_q = \varepsilon_q\psi_q \end{aligned}$

      (26)

      with

      $ \vec \alpha \cdot \vec k = \left( {\begin{array}{*{20}{c}} 0&{\vec \sigma \cdot \vec k}\\ {\vec \sigma \cdot \vec k}&0 \end{array}} \right),\quad {\psi _q} = \left( {\begin{array}{*{20}{c}} \varphi \\ \chi \end{array}} \right), $

      (27)

      and

      $ V_\sigma = g_\sigma^q\sigma,\quad V_\omega = g_\omega^q\omega_0,\quad V_\rho = g_\rho^q\rho_{03}, $

      (28)

      where $ \sigma $ , $ \omega_0 $ and $ \rho_{03} $ are the classical meson fields at mean-field limit. By replacing the small component ( $ \chi $ ) of the Dirac field in terms of the larger one ( $ \varphi $ ), and rewriting the equations as:

      $ \left[\frac{k^2}{2(\varepsilon_q^* + m_q^*)}+\frac{a}{2}r^2\right]\varphi(\vec{r}) = \frac{\varepsilon_q^* - m_q^* - V_0}{2}\varphi(\vec{r}), $

      (29)

      with

      $ \varepsilon_q^* = \varepsilon_q - V_\omega - \frac{1}{2}\tau_zV_\rho $

      (30)

      and

      $ m_q^* = m_q - V_\sigma. $

      (31)

      The Schrödinger equation of the tridimensional harmonic oscillator in Eq. (29) with the lowest order energy in the right-hand-side satisfy,

      $ \varepsilon_q^* - m^*_q - V_0 = 3\sqrt{\frac{a}{\varepsilon_q^*+m^*_q}}, $

      (32)

      by taking into account that $ \dfrac{3}{2}\omega = \dfrac{\varepsilon_q^* - m^*_q - V_0}{2} $ with $ \omega = \sqrt{\frac{a}{\varepsilon_q^* + m^*_q}} $ , in units of $ \hbar = c = 1 $ .

      We next follow the procedure used in Batista et al. (2002) [35] to extract the center of mass corrections from the nucleon observables. The center of mass energy is then given by,

      $ \varepsilon_{\rm cm} = \frac{3}{2}\frac{\alpha}{(\varepsilon_q^* + m_q^*)} \frac{(3+23\beta/6)}{(1+3\beta/2)^2}, $

      (33)

      with

      $ \alpha = \sqrt{a}(\varepsilon_q^* + m_q^*)^{1/2} $

      (34)

      and

      $ \beta = \frac{\alpha}{(\varepsilon_q^* + m_q^*)^2} = \sqrt{a}(\varepsilon_q^* + m_q^*)^{-3/2}. $

      (35)

      The effective nucleon mass in the medium as the center of mass correction energy of the three independent quarks is taken into account is expressed as in Batista et al. (2002) [35],

      $ M_N^* = 3\varepsilon_q^* - \varepsilon_{\rm cm}, $

      (36)

      and the mean squared nucleon radius, also corrected for center of mass effects, is written as in Batista et al. (2002) [35],

      $ \left<r_N^2\right> = \frac{1+5\beta/2}{\alpha(1+3\beta/2)}. $

      (37)

      The harmonic oscillator parameters a and $ V_0 $ are determined by imposing the vacuum values for $ M_N^* $ and $ \left<r_N^2\right> $ . Here we adopt $ M_N^*(n_B = 0) = 939 $ MeV and $ \left<r_N^2\right>(n_B = 0) = 0.8^2 $ fm $ ^2 $ .

      The equations of state (EoS) and field equations of the MQMC model are given as in the QMC one, by taking $ \Lambda_v = 0 $ . More specifically, the energy density and pressure of the MQMC model are given by Eqs. (9) and (10), respectively. The mean fields $ \omega_0 $ , $ \rho_{03} $ and $ \sigma $ are obtained as indicated in Eqs. (11) and (15), all of which with the restriction that the $ \Lambda_v $ parameter is set equal to zero. The free parameters $ G_\sigma^{q2}\equiv (g_\sigma^q/m_\sigma)^2 $ and $ G_\omega^2\equiv (g_\omega/m_\omega)^2 $ are found by imposing the nuclear matter saturation at $ n_B = n_0 = 0.15 $ fm $ ^{-3} $ with a binding energy of $ B/A = -16.4 $ MeV. Finally, $ G_\rho^2\equiv (g_\rho/m_\rho)^2 $ is determined by fixing a particular value for $ J\equiv{\cal{S}}(n_B = n_0) $ , with the symmetry energy given by Mishra et al. (2015) [56] and (2016) [57],

      $ {\cal{S}} = \frac{k_F^2}{6(k_F^2 + M_N^{*2})^{1/2}} + \frac{1}{8}G_\rho^2n_B, $

      (38)

      with $ k_F $ being the Fermi momentum. The quark mass $ m_q $ is used here to fix the incompressibility at the saturation density $ K_0 = K(n_B = n_0) $ , with $ K = 9\partial p/\partial n_B $ . Here we restrict the MQMC model to present the same J and $ K_0 $ values as those from the QMC models in Sec. II, see Table I. Such parametrizations are named as MQMC1 and MQMC2. We also generate a third one, namely, MQMC3 in which one has $ J = 25 $ MeV. Notice that all these parametrizations present J inside the range of $ 25\,\rm{MeV}\leqslant J \leqslant 35\,\rm{MeV} $ (Dutra et al. (2014) [42]). With regard to the incompressibility, all parametrizations present the same value of $ K_0 = 295 $ MeV, obtained by using the quark mass $ m_q = 210.61 $ MeV. Notice that $ K_0 $ is within the limits of $ 250\,\rm{MeV}\leqslant K_0 \leqslant 315\,\rm{MeV} $ (Stone et al. (2014) [58]).

      For the stellar matter calculations, we proceed as described in Sec. II concerning the $ \beta $ -equilibrium conditions on the chemical potentials and densities. In particular, the nucleon chemical potentials in the MQMC model are given by,

      $ \mu_{p,n} = (k_F^2 + M_N^{*2})^{1/2} + G_\omega^2(n_p+n_n)\,\pm \frac{1}{4}G_\rho^2(n_p-n_n) $

      (39)

      with the upper (lower) sign for protons (neutrons). The nuclear matter and stellar properties, along with the free parameters obtained from the parametrizations of the MQMC model, are given in Table. II.

    III.   GW170817 CONSTRAINTS
    • The gravitational Love number depends directly on the detailed structure of the neutron star (NS). Therefore the calculation of Love numbers can offer us paramount information on the NS composition. In fact this physics has triggered intense research recently (see e.g., Fattoyev et al. (2013) [59]; Malik et al. (2018) [60]; Hornick et al. (2018) [61]; Kumar et al. (2017) [62]). When one of the neutron stars in a binary system gets close to its companion just before merging, a mass quadrupole develops as a response to the tidal field induced by the companion. This is known as tidal polarizability (Damour & Nagar (2009) [ 63]; Binnington & Poisson (2009) [ 64]) and can be used to constrain the macroscopic properties of neutron star (Flanagan & Hinderer (2008) [ 65]), which in turn, are obtained from appropriate equations of state.

      In a binary system the induced quadrupole moment $ Q_{ij} $ in one neutron star due to the external tidal field $ {\cal E}_{ij} $ created by a companion compact object can be written as in Flanagan & Hinderer (2008) [ 65],

      $ Q_{ij} = -\lambda {\cal E}_{ij}, $

      (40)

      where, $ \lambda $ is the tidal deformability parameter, which can be expressed in terms of dimensionless $ l = 2 $ quadrupole tidal Love number $ k_2 $ as

      $ \lambda = \frac{2}{3} {k_2}R^{5}. $

      (41)

      To obtain $ k_{2} $ , we have to simultaneously solve the TOV equations and find the value of y from the following differential equation

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

      (42)

      with its coefficients given by

      $ F(r) = \frac{r - 4\pi r^3(\varepsilon - p)}{r-2m} $

      (43)

      and

      $ \begin{aligned}[b] Q(r) =& \dfrac{4\pi r \left[5\varepsilon + 9p + \dfrac{(\varepsilon + p)}{\partial p/\partial \varepsilon}-\dfrac{6}{4\pi r^2}\right]}{r-2m} \\ &- 4\left( \dfrac{m+4\pi r^3 p}{r^2-2mr} \right)^2, \end{aligned} $

      (44)

      where $ \varepsilon $ and p are the energy density and pressure profiles inside the star. Then we can compute the Love number $ k_2 $ , which is given by

      $ \begin{aligned}[b] k_2 =& \frac{8C^5}{5}(1-2C)^2[2+2C(y_R-1)-y_R] \Big\{2C [6-3y_R \\& + 3C(5y_R-8)]+4C^3[13-11y_R+C(3y_R-2) \\ & + 2C^2(1+y_R)] +3(1-2C)^2 \\ & \times [2-y_R+2C(y_R-1)]{\rm ln}(1-2C)\Big\}^{-1}, \end{aligned} $

      (45)

      where $ y_{R} = y(R) $ , and $ C = M/R $ is the compactness of the star of radius R and mass $ M = m(R) $ . The tidal deformability $ \Lambda $ (i.e., the dimensionless version of $ \lambda $ ) is connected with the compactness parameter C by,

      $ \Lambda = \frac{2k_2}{3C^5}. $

      (46)

      One can obtain the values of $ \Lambda_1 $ and $ \Lambda_2 $ using Eq. (46) by using corresponding quantities for each companion neutron star in the binary system.

    • A.   The inner and outer crust effects on the tidal polarizability

    • We would expect that the crust thickness and constitution affect the second Love number and consequently, the tidal polarizability. The neutron star crust is usually divided into two different parts: the outer crust and the inner crust. Piekarewicz and Fattoyev (2008) [66], investigated the impact of the crust by considering simple expressions. For the outer crust, the region where all neutrons are bound to finite nuclei, a crystal lattice calculation that depends on the masses of different nuclei was performed. For the inner crust, where the pasta phase is expected to exist, the authors used a polytropic EoS that interpolates between the homogeneous core and the outer crust. Their work concludes that, for fixed compactness, the second Love number is sensitive to the inner crust, but as the tidal polarizability scales as the fifth power of the compactness parameter, the overall impact is minor. In the present work, we use the BPS EoS from Baym et al. (1971) [67] for the outer crust and the pasta phase for the inner crust and test their effects on the deformability of the NS. We can see the differences between the total equations of state (outer crust + inner crust + liquid core) presented here and the one used in Piekarewicz & Fattoyev (2008) [ 66] in Fig. 1. One can see that the outer crusts are coincident and the liquid cores of all EoS are very similar. We would like to stress that as the QMC model is not used to calculate the outer crust, we employ an EoS generally found in the literature for this region. A better prescription for the outer crust is given in Fantina et al. (2020) [68] and we discuss the differences with our work when the results are presented in Section V. Hence, most of the differences seen in Fig. 1 reside in the inner crust and around the crust-core transition region. Notice that the pasta phase is present only in the QMC and QMC $ \omega\rho $ models. An inset with details of the inner crust, where the pasta phase lies, and the transition region from the pasta to the homogeneous phase is also included. More specifically, the pasta region is inside the ranges of $ n_B = (0.40-6.2)\times 10^{-2} $ fm $ ^{-3} $ (QMC model), and $ n_B = (0.45-6.6)\times 10^{-2} $ fm $ ^{-3} $ (QMC $ \omega\rho $ model). In the next section, we obtain the results for the second Love number and the tidal deformability and discuss them.

      Figure 1.  Stellar matter EoSs described by QMC and MQMC models in comparison with the FSUGarnet model [66]. The region in between the two stars on the orange curve represents the inner crust. Before (after) that region, one has the outer crust (liquid core).

    IV.   NEUTRON STAR COOLING
    • We now turn our attention to the thermal evolution of neutron stars whose composition is described by the models discussed in this manuscript. The cooling of neutron stars is governed by the emission of neutrinos from their core, and of photons from the surface. All thermal properties of the neutron star from the neutrino emission to the heat transport, depend on their microscopic composition, as different compositions lead to different cooling processes, or alter the rate at which certain processes take place. Global properties of the star such as size, mass, rotation also play an important role on the thermal evolution (Franzon et al. (2017) [19]; Negreiros et al. (2012) [29]; (2013) [30]; Yakovlev et al. (2000) [69]; Page et al. (2004) [70]), this makes cooling studies a superb way of bridging the gap between the micro and macroscopic realms of neutron stars.

      In recent years there has been great advances in the observation of thermal properties of compact objects (see e.g., Beloin et al. (2018) [71]; Safi-Harb & Kumar (2008) [ 72]; Zavlin (1999) [73], (2007) [74], (2009) [75]; Pavlov et al. (2001) [76], (2002) [77]; Mereghetti et al. (1996) [78]; Gotthelf et al. (2002) [79]; MvGowan et al. (2003) [80], (2004) [81], (2006) [82]; Klochkov et al. (2015) [83]; Possenti et al. (1996) [84]; Halpern et al. (1997) [85]; Pons et al. (2002) [86]; Burwitz et al. (2003) [87]; Kaplan et al. (2003) [88]; Zavlin et al. (2009); Ho et al. (2015) [89]). The equations that describe thermal energy balance and transport inside a spherically symmetric general relativistic star are given in Page et al. (2006) [90], Weber (1999) [91] and Schaab et al. (1996) [92], with ( $ G = c = 1 $ )

      $ \frac{ \partial (l e^{2\phi})}{\partial m} = -\frac{1}{ \varepsilon \sqrt{1 - 2m/r}} \left( \epsilon_\nu e^{2\phi} + c_v \frac{\partial (T e^\phi) }{\partial t} \right), $

      (47)

      $ \frac{\partial (T e^\phi)}{\partial m} = - \frac{(l e^{\phi})}{16 \pi^2 r^4 \kappa \varepsilon \sqrt{1 - 2m/r}}, $

      (48)

      In Eqs. (47)-(48) l is the luminosity, $ \phi $ is the metric function, m is the mass as a function of the radial distance r and T is the temperature. In addition to the macroscopic properties obtained by the solution of the Tolman-Opphenheimer-Volkoff (TOV) equations, one also needs to calculate thermal properties that regulate the cooling process. These are the neutrino emissivity ( $ \epsilon_\nu $ ), specific heat ( $ c_v $ ), and thermal conductivity ( $ \kappa $ ). In this work, we consider all established cooling mechanism for the thermal evolution of neutron stars and do not consider any exotic process such as axion cooling (Sedrakian (2016) [93]) and dark matter decay, for instance (Kouvaris (2008) [94]). We now describe the relevant thermal quantities and also refer the reader to Yakovlev et al. (2000) [69] for a thorough review on the subject.

      The neutrino emissivities are given by the following processes:

      ● Direct Urca (DU) Process:

      $ n \rightarrow p + e^- + \bar{\nu} $ ,

      $ p + e^- \rightarrow n + e^- + \nu $ .

      It consists of neutron beta decay and electron capture by the proton. This process is extremely efficient in extracting thermal energy from the neutron star and its emissivity is $ \epsilon_{DU} \sim 10^{27} T_9^6 $ ergs cm $ ^{-3} $ s $ ^{-1} $ (where $ T_9 = T/10^9 $ K) (Prakash et al. (1992) [95]). It can only take place when the three triangle inequalities $ k_{fi} +k_{fj} \geq k_{fk} $ (where i, j, and k are the baryons and leptons, respectively) are satisfied, as to guarantee momentum conservation (Prakash et al. (1992) [95]).

      ● Modified Urca (MU) Process:

      $ n + n' \rightarrow p + n' + e^- + \bar{\nu} $ ,

      $ p+n' + e^- \rightarrow n + n' + \nu $ . The MU is similar to the DU, except it contains a bystander particle (denoted by prime) that facilitates momentum conservation. The phase space reduction also limits the efficiency of such process (when compared to the DU). The emissivity of the MU is $ \epsilon_{MU} \sim 10^{21} T_9^8 $ ergs cm $ ^{-3} $ s $ ^{-1} $ . We note here that we adopt the MU emissivities as calculated in (Friman & Maxwell (1979) [ 96]). It should be noted however that it has been discussed that medium effects due to pions at higher densities may alter such emissivities (Migdal et al. (1990) [97]; Schaab et al. (1997) [98]; Blaschke et al. (2004) [99]; (2013) [100]).

      ● Bremssthralung (BR) process:

      $ n + n' \rightarrow n + n' +\nu + \bar{\nu} $ ,

      $ n + p' \rightarrow n + p' +\nu + \bar{\nu} $ ,

      $ p + p' \rightarrow p + p' +\nu + \bar{\nu} $ .

      These scattering processes have also been calculated in Friman & Maxwell (1979) [ 96] and have emissivities $ \epsilon_{nn} \sim 10^{19} T_9^8 $ ergs cm $ ^{-3} $ s $ ^{-1} $ , $ \epsilon_{np} \sim 10^{20} T_9^8 $ ergs cm $ ^{-3} $ s $ ^{-1} $ and $ \epsilon_{pp} \sim 10^{20} T_9^8 $ ergs cm $ ^{-3} $ s $ ^{-1} $ .

      ● Pair Breaking and Formation (PBF) process:

      $ \{ BB \} \rightarrow B + B +\nu + \bar{\nu} $ ,

      $ B + B \rightarrow \{ BB \} +\nu + \bar{\nu} $ .

      (B stands for the baryon in question)

      This process is associated with the breaking and formation of pairs as the temperature nears the critical temperature for pairing formation. This process is transient by definition, only taking place near the pairing onset. The emissivity depends on the type of pairing, with $ \epsilon_{PBF} \sim 10^{21} T_9^7 $ , $ \epsilon_{PBF} \sim 10^{22} T_9^7 $ , and $ \epsilon_{PBF} \sim 10^{21} T_9^7 $ for proton singlet ( $ ^1S_0 $ ), neutron singlet ( $ ^1S_0 $ ), and neutron triplet ( $ ^3P_2 $ ) pairing, respectively (Yakovlev et al. (2000) [69]).

      It is important to note that the processes discussed above are subjected to suppression once the temperature is low enough for pairing to set in. This is due to the presence of an energy gap ( $ \Delta $ ) in the particle energy spectrum. Generally the emissivities are suppressed by a factor $ \sim e^{-\Delta(T)/T} $ (Levenfish et al. (1994) [25]).

      As for the specific heat, one calculates the contribution of all constituents by employing traditional microscopic calculations, which leads to the specific heat of species i as

      $ c_{i,v} = \frac{m*_i n_i}{k_{f,i}^2}\pi^2 T, $

      (49)

      where $ m*_i $ is the effective mass of the i-species. It is important to note that the specific heat is also subject to suppression due to pairing, much like neutrino emission processes, albeit with a more complicated mathematical form (Levenfish et al. (1994) [25]).

      The last ingredient is the thermal conductivity for neutron star matter, that has been thoroughly discussed in Flowers & Itoh (1976) [ 101]; (1979) [102]; (1981) [103], whose fit for the high density regime is $ \kappa \sim 10^{23} \rho_{14}/T_8 $ ergs cm $ ^{-1} $ s $ ^{-1} $ K $ ^{-1} $ (where $ \rho_{14} = \rho/10^{14} $ g/cm $ ^3 $ ).

      Once we have all the microscopic and macroscopic ingredients we solve Eqs. (47)-(48) numerically. In order to do that we need two boundary conditions:

      ● 1. $ l (r = 0) = 0 $ - which is necessary as there is no heat flow at the star core.

      ● 2. $ T_e = T_e(T_b(R)) $ which is the dependence of the effective surface temperature at the star envelope as a function of the mantle ( $ T_b $ ) temperature at $ r = R $ . Such relationship has been thoroughly studied by Gudmundsson et al. (1982) [104], (1983) [105]), accounting for the properties of the star envelope, as well as the possibility of accreted matter due to fallback at the neutron star formation process.

    V.   RESULTS AND DISCUSSION
    • In the present section, we study the effects of the different versions of the QMC and MQMC models on quantities related to gravitational wave observables and on the neutron star cooling.

    • A.   GW170817 constraint results

    • The results of the tidal deformability $ \Lambda $ and the Love number $ k_2 $ obtained from the QMC/QMC $ \omega\rho $ models, with and without pasta phases, and the three versions of the MQMC are shown in Fig. 2, with the range for $ \Lambda_{1.4} $ [106] also displayed. In the top and bottom panel of Fig. 2, we display the Love number $ k_2 $ as a function of the compactness of the neutron star, and the variation of $ \Lambda $ with the mass of the star, respectively. As already pointed out in Piekarewicz & Fattoyev (2019) [ 66], the curves do not collapse into one single curve because of their different dependence on y. When the same models are used to compute the tidal polarizability, their differences are enhanced and only two models, namely, QMC $ \omega\rho $ and MQMC3 give results that are consistent with LIGO results for the canonical star. The pasta phase plays a very modest role and its influence is practically unnoticed, what can be seen if one compares the curves obtained from QMC with and without pasta and QMC $ \omega\rho $ with and without pasta. As far as the outer crust is concerned, all prescriptions, namely, BPS, crystal lattice calculation used in Piekarewicz & Fattoyev (2019) [ 66] and the state-of-the-art calculation done in Fantina et al. (2020) [68] (tested, but not shown in this work) result on practically the same tidal deformabilities.

      Figure 2.  (top) $k_2$ as a function of the neutron star compactness, and (bottom) $\Lambda$ as a function of M. Full circle: recent result of $\Lambda_{1.4}=190_{-120}^{+390}$ obtained by LIGO and Virgo Collaboration [106].

      In the recent literature, a lot of effort has been put to constrain the radius of the canonical neutron star, see for instance, Malik et al. (2008) [60], Most et al. (2018) [107], Raithel et al. (2018) [108] and Annala et al. (2018) [109] for studies in which the outcomes of the GW170817 event are used to infer the radius. In Malik et al. (2008) [60], the authors have discussed this constraint in the context of the Skyrme and RMF models and their calculations suggest the range of $ 11.82\; \rm{km} \leqslant R_{1.4}\leqslant 13.72\; \rm{km} $ . In Most et al. (2018) [107], the constraint on $ R_{1.4} $ was established from a large number of EoS with pure hadronic matter without any kind of phase transition. The authors found the value of $ R_{1.4} $ inside the range of $ 12.00\; \rm{km} \leqslant R_{1.4} \leqslant 13.45\; \rm{km} $ , with the most likely value of $ R_{1.4} = 12.39\; \rm{km} $ . The upper limit of this range is compatible with the findings of Raithel et al. (2018) [108], namely, the radius of a neutron star cannot be larger than around $ 13 $ km. A similar result was obtained in Annala et al. (2018) [109], in which the authors conclude that the maximal radius of a canonical neutron star is $ 13.6 $ km. Concerning the parametrizations that satisfy the constraint for $ \Lambda_{1.4} $ shown in the bottom panel of Fig. 2, namely, QMC $ \omega\rho $ and MQMC3, it is found that they predict $ R_{1.4} = 12.83 $ km and $ 12.94 $ km, respectively (see Tables I and II). These values are compatible with the aforementioned results.

      Quantity Model
      MQMC1 MQMC2 MQMC3
      J (MeV) 34.50 30.92 25.00
      $L_0$ (MeV) 93.20 82.46 64.70
      a fm $^{-3}$ 0.95 0.95 0.95
      $V_0$ (MeV) -92.27 -92.27 -92.27
      $G_\sigma^{q2}$ 5.13 5.13 5.13
      $G_\omega^2$ 8.33 8.33 8.33
      $G_\rho^2$ 14.67 12.18 8.08
      $M_{max}$ (M $_{\odot}$ ) 1.97 1.97 1.97
      $R_{M_{max}}$ (km) 11.43 11.34 11.18
      $R_{1.4}$ (km) 13.55 13.32 12.94

      Table 2.  Nuclear matter and stellar properties obtained from the parametrizations of the MQMC model. The free parameters are also given. $G_\sigma^{q2}$ , $G_\omega^2$ and $G_\rho^2$ are given in $10^{-5}$ MeV $^{-2}$ . For all parametrizations, one has $B/A=-16.4$ MeV, $n_0=0.15$ fm $^{-3}$ , $m_q = 210.61$ MeV, $K_0 = 295$ MeV, and $M^*_N/M_N = 0.84$ .

      In Fig. 3 we show the tidal deformabilities of each neutron star in the binary system. $ \Lambda_1 $ is associated to the neutron star with mass $ m_1 $ , which corresponds to the integration of every EoS in the range $ 1.37 \leqslant m/M_{\odot} \leqslant1.60 $ obtained from GW170817. The mass $ m_2 $ of the companion star is determined by solving $ {\cal{M}}_c = 1.188 M_{\odot} = \frac{(m_1 m_2)^{3/5}}{(m_1+m_2)^{1/5}} $ (Abbott et al. (2017) [111]). We notice that only the QMC $ \omega\rho $ model (with and without pasta phases), and the MQMC3 parametrization show values in between the confidence lines. These specific models satisfy the predictions from LIGO and Virgo collaboration due to the $ R_{1.4} $ dependence of $ \Lambda_{1.4} $ . Different studies point out that $ \Lambda_{1.4}\propto R_{1.4}^\alpha $ but with $ \alpha\ne 5 $ as Eq. (46) could suggest. In Lourenço et al. (2019) [110], for instance, where a previous analysis of 35 relativistic mean-field models was performed, the authors found $ \alpha = 6.58 $ . Therefore, parametrizations presenting lower values of $ R_{1.4} $ , in general, lead to a decreasing of $ \Lambda_{1.4} $ . In the case of the QMC/MQMC models, such a decreasing of $ \Lambda_{1.4} $ favors the parametrizations to satisfy the constraint shown in Fig. 2 (this is the case of QMC $ \omega\rho $ and MQMC3 models). Regarding Fig. 3, we marked the points (circles and squares) in which one of the companion stars of the binary system presents the value of $ m_1 = 1.4M_\odot $ . Therefore, in these points, one has $ \Lambda_1 = \Lambda_{1.4} $ . We can also see from this figure that lower values of $ \Lambda_{1.4} $ , caused by lower values of $ R_{1.4} $ , also induce the curves to be consistent with the LIGO and Virgo band constraint.

      Figure 3.  Tidal deformabilities obtained from the QMC, QMC $\omega\rho$ and MQMC models for both components of the binary system related to the GW170817 event. The confidence lines (50% and 90%) are the recent results of LIGO and Virgo collaboration taken from Abbott et al. (2018) [106]. The shaded region represents the results obtained with consistent relativistic mean field models in Lourenço et al. (2019) [110]. Circles (QMC model) and squares (MQMC model): points in which $m_1=1.4M_\odot$ and, consequently, $\Lambda_1=\Lambda_{1.4}$ .

    • B.   Cooling results

    • We now present the results of our thermal evolution studies, obtained by the numerical solution of Eqs. (47) and (48). We begin by showing the cooling of a set of stars covering a wide range of masses for the QMC model. Those are shown in Fig. 4. As shown in Fig. 4 all stars exhibit fast cooling, characterized by a sharp drop in their surface temperature at the age of $ \sim 100 $ years. Such fast cooling is a manifestation of a prominent presence of the direct Urca process (DU) in the stellar core, which is indeed the case for these stars. The DU process is extremely efficient in exhausting the star thermal energy, leading to the manifested fast cooling (Prakash et al. (1992) [95]). This, in turn indicates that such thermal evolution is mostly incompatible with observed data (unless the DU process is suppressed), as we discuss in the following.

      Figure 4.  Thermal evolution of neutron stars described by the QMC model. Y-axis represents the red-shifted surface temperature and x-axis the age. Each curve represents the thermal evolution of a neutron star with a different mass. All objects exhibit fast cooling due to prominent DU process in their core.

      We now show in Fig. 5 the thermal evolution of QMC neutron stars with the pasta phase describing the inner crust. We can see that the pasta phase has little effect on the overall cooling behavior of the star. This is expected as the pasta phase occupies only a small region. The more careful analysis leads us to the conclusion that the pasta phase had the minor effect of aiding in the core-crust thermalization, leading to a core-crust thermal relaxation on average $ \sim 3.5 $ years faster. This seems reasonable, as the pasta phase smooths the core-crust transition, thus, one can expect that will also smooth out heat propagation in between these regions. We note, however, that these are rough estimates and a more detailed calculation is warranted, one in which the thermal and conductivity properties of the pasta phase are explored in more detail.

      Figure 5.  Same as in Fig. 4 but for stars with pasta phase in the inner crust.

      We now repeat the calculations above, but for the QMC $ \omega \rho $ model, without pasta (Fig. 6) and with pasta (Fig. 7).

      Figure 7.  Same as Fig. 6 but for stars with a pasta phase.

      We can see that, as happens in the previous case (QMC), the pasta phase has little effect on the overall cooling - serving only to delay the thermalization by a few years. Differently than the previous model, however, in this case, all stars exhibit slow cooling, meaning that there is no sharp drop in temperature when the core and crust become thermally coupled. The reason for that is that in this particular model, the proton fraction at the core of the star, even at larger densities of heavier stars, is low enough to prevent the DU process from taking place, thus leading to a substantially slower thermal evolution.

      Figure 9.  Same as Fig. 8 but for the MQMC2 model.

      Now we investigate the cooling of neutron stars described by the MQMC model to see if the different prescriptions used in the EoS calculations play any role in the cooling. We show in Figs. 8 - 10 the cooling of neutron stars of different masses calculated for the different MQMC models studied.

      Figure 8.  Thermal evolution of different neutron stars under the MQMC1 model.

      Figure 10.  Same as Fig. 8 but for the MQMC3 model.

      We can see from the results shown in Figs. 8 - 10 that the overall cooling behavior associated with each model is, qualitatively, the same. For the MQMC models 1 and 2, we see that all stars exhibit fast cooling, with MQMC 3 displaying slow cooling (except for high mass stars). It is worth mentioning that we see a clear correlation between the slope of the nuclear symmetry energy ( $ L_0 $ ) and a fast/slow cooling behavior, with stars with higher values of $ L_0 $ having faster cooling (QMC, MQMC1, and MQMC2) whereas lower values of $ L_0 $ leading to slower cooling scenario (QMC $ \omega \rho $ and MQMC3).

      So far in our thermal investigation, we have completely ignored the effects of pairing, as we were interested in probing possible differences in the thermal behavior of the different models studied, and the inclusion of pairing could potentially murk the results. As it turns out all models behave similarly, with major differences connected with the different values of the symmetry energy slope rather than more specific minutia about the models. We now need to make sure that our model is in agreement with observed data, or at the very least is comparable to other results in the recent literature (Negreiros et al. (2018) [15]; Beloin et al. (2018) [71]; Carvalho et al. (2018) [112]; Dexheimer et al. (2015) [113]; Raduta et al. (2017) [114]) - for that we need to include pairing. Pairing in neutron stars may potentially occur in three different manners, neutron-neutron singlets ( $ ^1S_0 $ ), neutron-neutron triplet ( $ ^3P_2 $ ) and proton-proton singlet ( $ ^1S_0 $ ). Neutron-neutron singlets are known to form mainly in the neutron-free region of the crust, whereas its triplet counterpart should form mainly in the lower densities of the core. As for the proton-proton pairs, there is still great uncertainty as to how deep into the star it may take place. For a review of pairing in neutron stars, we refer the reader to Yakovlev et al. (2000) [69]. Our implementation of pairing follows the same principles as in Beloin et al. (2018) [71], i.e., given the current uncertainties of pairing at high-density neutron star matter, we use thermal observation data of neutron stars to constrain the neutron and proton pairing at the neutron star core. This allows us to calculate the pairing effects in the neutron star cooling processes [suppression of the DU process, the appearance of Pair-Breaking-Formation (PBF) process near $ T_c $ , and the modification of the specific heat of paired particles (Yakovlev et al. (2000)[69])].

      For completeness sake, we have shown the cooling of all models studied so far. For our final cooling (with the inclusion of pairing), and comparison with observed data we next focus solely on the models that have met the observational constraints discussed so far, namely the observed mass and tidal deformability, which are QMC- $ \omega \rho $ and MQMC3. As shown in Figs. 7 and 10, these models show very different thermal evolution, with QMC $ \omega \rho $ displaying a total absence of DU, and thus a very slow cooling, whereas MQMC3 has more balanced behavior, with low mass stars displaying slow cooling while heavier objects cooling down faster.

      Fig. 11 shows the thermal evolution of the 1.4 and 2.0 $ M_\odot $ stars under model MQMC3, with nucleon pairing taken into account. Intermediate mass stars have their thermal evolution in between these curves. Also shown are the most prominent observed data regarding cooling of neutron stars (Beloin et al. (2018) [71]; Safi-Harb & Kumar (2008) [ 72]; Zavlin (1999) [73], (2007) [74], (2009) [75]; Pavlov et al. (2001) [76], (2002) [77]; Mereghetti et al. (1996) [78]; Gotthelf et al. (2002) [79]; MvGowan et al. (2003) [80], (2004) [81], (2006) [82]; Klochkov et al. (2015) [83]; Possenti et al. (1996) [84]; Halpern et al. (1997) [85]; Pons et al. (2002) [86]; Burwitz et al. (2003) [87]; Kaplan et al. (2003) [88]; Zavlin et al. (2009); Ho et al. (2015) [89]). As expected nucleon pairing leads to slower cooling rates (associated with the suppression of neutrino emissivities, as discussed above). Thus, when compared with stars without pairing (see Fig. 10)) one gets a warmer object at the same age. From Fig. 11 we can note that to explain most observed data, within the confines of the cooling description adopted in this work, namely making use of the neutrino emissivities as described by Friman and Maxwell [96], the mass of the stars must be within the interval 1.4– 1.8 $ M_\odot $ (the latter indicated by the cyan curve in Fig. 12). One also notices, however, a few stars below such curve, in which case their mass would be slightly above 1.8 $ M_\odot $ , which would indicate a group of observed stars within a small mass window, which may be unlikely. This poses a challenge for the description of the data making use of the Friman and Maxwell description for the Urca process.The objects that fall above the expected range can be explained by considering fallback in the early stages of neutron star evolution (Gudmundsson et al. (1983) [105]). Such fall back alters the composition of the atmosphere from pure iron, to other elements, which affect the thermal evolution, keeping the star warmer for longer.

      Figure 11.  Thermal evolution of a set of stars with different masses stars of model MQMC3, under the effect of nucleon pairing, as well as observed thermal data of cooling neutron stars.

      Figure 12.  Thermal evolution of a set of stars with different masses under MQMC3 model with pairing and by taking into account different atmosphere compositions (represented by the amount of accreted matter $\Delta M/M$ ). The blue hashed region represents the different possible thermal evolution that encompasses different masses and different atmospheric compositions.

      Fig. 12 shows the thermal evolution of stars under the MQMC3 model, with pairing, as well as considering different possibilities of fallback accretion ( $ \Delta M/M = 10^{-9} $ and $ 10^{-11} $ ). The hashed band represents the different possible thermal evolution that encompasses different masses and different possible atmospheres.

      We can now conclude, based on the results of Fig. 11 and 12 that the MQCM3 model not only complies with the constraints set by neutron star mass measurements, and those extracted from tidal deformability, as it is also in agreement with the thermal data of neutron stars.

      We now turn our attention to the QMC $ \omega \rho $ (pasta) model, whose cooling of 1.4 and 2.0 $ M_\odot $ stars, under the effect of pairing, is shown in Fig. 13. As was the case for stars without pairing (Fig. 6) stars in this model exhibit slow cooling. The presence of pairing does broaden the cooling "spectrum" (meaning the different cooling tracks associated with stars of different masses) but it still leaves out a large portion of observed data.

      Figure 13.  Thermal evolution of the 1.4 and 2.0 $M_\odot$ stars of model QMC $\omega \rho$ (pasta), under the effect of nucleon pairing, as well as observed thermal data of cooling neutron stars.

      For completeness, we also investigate the effect of different atmospheric compositions of the stars of the QMC $ \omega \rho $ (pasta) model under the effect of pairing, which is shown in Fig. 14. As was the case for the pure iron atmosphere, the thermal evolution of stars in this model is too slow to explain satisfactory all data, failing to describe the observed temperature of colder stars.

      Figure 14.  Thermal evolution of stars under QMC $\omega \rho$ (pasta) model with pairing and by taking into account different atmosphere compositions (represented by the amount of accreted matter $\Delta M/M$ ). The blue hashed region represents the different possible thermal evolution that encompasses different masses and different atmospheric compositions.

    VI.   FINAL REMARKS
    • In the present work, we have used five different versions of the quark-meson coupling model to compute the Love number and tidal polarizability and confront them with GW170817 constraints and also to investigate the cooling process they describe. Two of the models are based on the original bag potential structure (QMC model) and three versions consider a harmonic oscillator potential to confine the quarks (MQMC model). The bag-like models also incorporate the pasta phase used to describe the inner crust of neutron stars. Although based on different concepts, all versions of the quark-meson coupling model used in the present work were chosen because they satisfy nuclear bulk properties within acceptable ranges. We have compared our EoS with the FSUGarnet (Piekarewicz & Fattoyev (2019) [ 66]) and checked that the outer crusts are coincident and the liquid core of all EoS are very similar, most of the differences residing on the inner crust and around the crust-core transition region, where the pasta phase was included. Our results point to the fact that the pasta phase plays only a minor role in the calculation of the dimensionless tidal polarizability and just two of the models (QMC $ \omega\rho $ and MQMC3) give results that lie within the expected range of the canonical star $ \Lambda_{1.4} $ and, at the same time, appear inside the region delimited by the confidence lines. At this point, we emphasize that there is a correlation between the model symmetry energy slope and the fact that it satisfies (or not) the constraints investigated, such that the lower values seem to be more consistent with the observational data. Furthermore, we can see that the two models that satisfy the LIGO and Virgo predictions, present similar patterns for the radii and polarisability of the canonical star. It was already known (Lourenço et al. (2019) [110]) that, in general, parametrizations presenting lower values of $ R_{1.4} $ lead to a decreasing of $ \Lambda_{1.4} $ , and among the models examined, QMC $ \omega\rho $ and MQMC3 models, which have lower values of $ \Lambda_{1.4} $ are the ones consistent with the LIGO and Virgo band constraint. It is worth mentioning that the chirp mass value was updated to $ {\cal{M}}_c = 1.186 M_{\odot} $ (Abbott et al. (2019) [115]) after this work was completed, but new calculations for $ m_2 $ would not affect our main results and conclusions.

      We have then used the very same models to study the cooling process and have verified that, in this case, there is a clear correlation between the slope of the symmetry energy and the velocity of the cooling process, i.e., models with higher (lower) values of the slope produce fast (slow) cooling. We note that this result seems to agree with an analysis done in Dexheimer et al. (2015) [116], in which the authors update the non-linear $ \sigma $ model. In this work a similar trend was found, i.e. models with larger values for the slope of the symmetry energy exhibited faster cooling than their counterparts with lower values. The pasta phase, once again, plays a minor role in the cooling process by speeding up the process in about 3.5 years if superconductivity is not considered.

      We have further explored the thermal evolution, under a more realistic scenario in which nucleon pairing as well as different atmospheric models were considered. By taking such phenomena into account, we calculated the thermal evolution of the two models (MQMC3 and QMC $ \omega \rho $ ) that were in agreement with the constraints set by neutron star mass measurements, as well as the tidal deformability estimates from GW170817. Our results were confronted with observed thermal data on cooling neutron stars, which showed that the MQMC3 model can successfully explain the whole set of observed data (by encompassing the different possible masses and atmospheric compositions). Model QMC $ \omega \rho $ , on another hand, can only match the observations of higher temperature neutron stars, leaving the colder objects out of its prediction range. It is important to note, however, that there is still room for further investigation, by considering, for instance, in medium corrections to the modified Urca process ((Migdal et al. (1990) [97]; Schaab et al. (1997) [98]; Blaschke et al. (2004) [99]; (2013) [100])) as well as the effects of pairing on the electron thermal conductivity [117]. We intend to pursue such investigations in future research.

      Our studies show that the only model out of the five chosen ones in the present paper that can, at the same time, describe astrophysical quantities as massive NS and tidal polarizabilities and a perfect description of the possible cooling processes is the MQMC3 model. However, more observational data is desirable. It is also important to note that such a conclusion was drawn taking into account neutrino emissivities as calculated by Friman and Maxwell [96]. A different approach for the Urca emissivity, taking into account in medium corrections, as mentioned above, could potentially lead to different conclusions. We would like to point out that the adoption of a purely nucleonic EoS, such as the one we adopted in this work is fairly usual (see the recent works from Bauswein et al. (2019)[118] and Lucca & Sagunski (2019) [ 119], for instance), which is not to say that the consideration of other phases of matter such as quark matter is unusual. There are many other avenues concerning the strongly interacting matter that may still be considered in the context of our work, for instance, stellar matter consisting of hadronic and quark EoS bearing or not a mixed phase, and at higher densities a deconfined quark matter in the CFL and CCS phases (Lau et al. (2017) [120]). It has recently been pointed out that the existence of quark cores inside massive stars should be considered the standard pattern and not an exotic alternative (Annala et al. (2020) [121]). This conclusion was based on a model-independent analysis of the sound velocity in both hadronic and quark matter. Hence, once the possibilities of these different phases are considered, the internal constitution would certainly affect the quantities studied in the present work and every case is worth investigating and confronting with observational constraints. The NICER telescope [122], launched in 2017 is expected to test nuclear physics theory by exploring the exotic states of matter within neutron stars through rotation-resolved X-ray spectroscopy and, hence, provide precise values for masses and radii. If associated with future GW results, the identification of the inner constituents of neutron stars in the binary system (neutron-neutron, neutron-hybrid or hybrid-hybrid) (Christian et al. (2019) [123]) becomes a real possibility. The inclusion of other internal phase transitions and their consequences on the thermal evolution of neutron stars is then postponed to future studies. In that sense, the present work is the first step towards calculations including a pasta-like phase possibly present at higher stellar densities, in between hadronic and quark matters existing in hybrid stars.

    ACKNOWLEDGMENTS
    • D.P.M. thanks Dr. Prafulla K. Panda for useful comments and discussions on the parametrizations of the QMC model and Thomas Carreau and Francesca Gulminelli for providing their outer crust EOS used for tests and comparisons with the results presented in this work.

Reference (123)

目录

/

DownLoad:  Full-Size Img  PowerPoint
Return
Return