We have developed and tested the novel N-body simulation code BIFROST
specialising in star cluster simulations with a large fraction of binary and
multiple star systems. Based on our earlier GPU-accelerated hierarchical
fourth-order forward integrator code FROST the new code includes both
regularised (LOGH, MSTAR) and secular integration (binary, triple) methods for
compact subsystems.

In the regularised subsystems, the equations of motion of the simulation
particles contain post-Newtonian terms up to PN3.5 order including three-body
Einstein–Infeld–Hoffmann terms in PN1.0. Spin-dependent PN terms can be
enabled by the user as well. In the secular integration the highest orbit-
averaged post-Newtonian term is PN2.5. For the FSI integration we have the
possibility to include the so-called global PN1.0 term which enables post-
Newtonian dynamics also outside subsystems.

In addition to the new subsystem integrator, BIFROST now includes single
stellar evolution using the SSE tracks and prescriptions of various particle
merger processes. We include gravitational-wave inspiral mergers for compact
binary systems, tidal disruption events between compact objects and stars, and
finally a prescription for stellar mergers. For binary black hole mergers we
include relativistic mass loss and gravitational-wave recoil kicks using
fitting functions from numerical relativity simulations.

We have presented the timing and scaling tests for the BIFROST code simulating
star cluster models between $N=10^{4}$ and $N=5\times 10^{6}$ particles. The
binary fractions of the star clusters range from $f_{\mathrm{bin}}=10\%$ to
$f_{\mathrm{bin}}=100\%$ with the highest number of binary systems being
$N_{\mathrm{bin}}=2.5\times 10^{6}$. For particle numbers above $N\gtrsim
10^{5}$ the binary fraction $f_{\mathrm{bin}}$ has only a small effect on the
scaling and running times of the simulations. We confirm that the strong
scaling of the novel BIFROST code is very similar as the scaling of its
precursor code FROST. The maximum number of GPUs to be used in a simulation
with $N$ particles is still approximately $N_{\mathrm{GPU}}\sim 40\times
N/10^{6}$. The maximum number of simulation particles for BIFROST is
approximately a few million, depending on the density and compactness of the
system as well as the desired simulation time.

We have tested the accuracy of our new code in a series of star cluster
simulations. The various star cluster simulations examine the effect of the
stellar IMF and the binary fraction on the amount of accumulated energy error.
In general, simulations with more massive particles and higher binary
fractions are less accurate due to more frequent close and strong encounters.
By carefully choosing the user-given code accuracy parameters the total
relative energy error was always less than $|\Delta E/E_{\mathrm{tot}}|\sim
6\times 10^{-7}$ in the runs. We note that the energy error can be further
decreased by decreasing the time-step accuracy parameters $\eta$ at the cost
of increased code running times.

Finally, we have performed more simulations of more extreme stellar-dynamical
scenarios. We modelled the inspiral and rapid merger of an IMBH binary
embedded in a star cluster. The merger remnant black hole received a strong
gravitational-wave recoil kick of over $v=200$ km/s and soon escaped the star
cluster. The total relative energy error of the simulation was approximately
$|\Delta E_{\mathrm{tot}}/E|\sim 10^{-6}$. In order to further test of our
code we simulated the core collapse of a star cluster with $N=20000$,
$r_{\mathrm{h}}=0.5$ pc and $f_{\mathrm{bin}}=10\%$. The total relative energy
error of the core collapse simulation was $|\Delta E_{\mathrm{tot}}/E|\sim
3\times 10^{-4}$, a result comparable to the accuracy of other collisional
simulation codes in the literature (e.g. Konstantinidis & Kokkotas 2010;
Pelupessy et al. 2012; Wang et al. 2020b. We also performed a run of the same
initial conditions using NBODY6++GPU in order to directly compare the global
properties of the star cluster simulated with BIFROST with another code. The
running times and energy errors of the particular star cluster model with the
two codes are very similar. Most importantly, star cluster properties
including the Lagrangian radii and the binary star population evolve in an
analogous manner.

Hierarchical fourth-order forward symplectic integrators are a promising
alternative for both Hermite integrators in traditional direct-summation
N-body simulations and second-order leapfrog-type symplectic integrators,
especially when supplemented with accurate subsystem integrators. Efficient
parallelization of our BIFROST code allows for simulating stellar systems with
arbitrary binary fractions up to $f_{\mathrm{bin}}=100\%$, a feature which
only few other simulation codes include (Wang et al., 2020b). In addition to
the arbitrary binary fractions we highlight the capability of BIFROST to
simulate relativistic binary systems. This is possible as both the secular and
regularised subsystem integrators in BIFROST contain post-Newtonian terms in
their equations of motion.

The version of BIFROST presented in this study is an efficient and versatile
tool for simulating stellar systems from few-body systems to massive star
clusters and low-mass nuclear star clusters including the inner parsecs of the
Milky Way centre. The immediate near-future code updates including binary
stellar evolution and hydrodynamics will upgrade BIFROST to an even more
powerful tool for studying the formation, evolution and possible disruption of
star clusters in a wide range of mass scales. In addition, incorporating a
tree or FMM gravity solver will eliminate the current wall-clock time
bottleneck placed by the direct summation approach in the near future. This
will enable simulating the evolution of the most massive star clusters in the
Universe, and even full low-mass isolated dwarf galaxies and their mergers.

## Data availability statement

The data relevant to this article will be shared on reasonable request to the
corresponding author.

## Acknowledgments

The numerical simulations were performed using facilities hosted by the Max
Planck Computing and Data Facility (MPCDF) and the Leibniz Supercomputing
Centre (LRZ), Germany. TN acknowledges support from the Deutsche
Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s
Excellence Strategy - EXC-2094 - 390783311 from the DFG Cluster of Excellence
"ORIGINS". F.P.R. acknowledges the support by the European Research Council
via ERC Consolidator Grant KETJU (no. 818930).

## References

  * Aarseth (1966) Aarseth S. J., 1966, MNRAS, 132, 35
  * Aarseth (1971) Aarseth S. J., 1971, Astrophysics & Space Science, 14, 118
  * Aarseth (2003) Aarseth S. J., 2003, Gravitational N-Body Simulations. Cambridge University Press
  * Aarseth (2012) Aarseth S. J., 2012, MNRAS, 422, 841
  * Aarseth & Zare (1974) Aarseth S. J., Zare K., 1974, Celestial Mechanics, 10, 185
  * Abbott et al. (2016) Abbott B. P., et al., 2016, Phys. Rev. Lett., 116, 061102
  * Abbott et al. (2017) Abbott B. P., et al., 2017, Phys. Rev. Lett., 119, 161101
  * Amaro-Seoane (2018) Amaro-Seoane P., 2018, Living Reviews in Relativity, 21, 4
  * Avramov et al. (2021) Avramov B., Berczik P., Just A., 2021, arXiv e-prints, p. arXiv:2109.08469
  * Banerjee et al. (2020) Banerjee S., Belczynski K., Fryer C. L., Berczik P., Hurley J. R., Spurzem R., Wang L., 2020, A&A, 639, A41
  * Barnes & Hut (1986) Barnes J., Hut P., 1986, Nature, 324, 446
  * Belczynski et al. (2010) Belczynski K., Bulik T., Fryer C. L., Ruiter A., Valsecchi F., Vink J. S., Hurley J. R., 2010, ApJ, 714, 1217
  * Belczynski et al. (2016) Belczynski K., et al., 2016, A&A, 594, A97
  * Benacquista & Downing (2013) Benacquista M. J., Downing J. M. B., 2013, Living Reviews in Relativity, 16, 4
  * Bettwieser & Sugimoto (1984) Bettwieser E., Sugimoto D., 1984, MNRAS, 208, 493
  * Blanchet (2014) Blanchet L., 2014, Living Reviews in Relativity, 17, 2
  * Blanchet & Iyer (2003) Blanchet L., Iyer B. R., 2003, Classical and Quantum Gravity, 20, 755
  * Brown & Gnedin (2021) Brown G., Gnedin O. Y., 2021, MNRAS, 508, 5935
  * Bulirsch & Stoer (1966) Bulirsch R., Stoer J., 1966, Numerische Mathematik, 8, 1
  * Chatterjee et al. (2013) Chatterjee S., Rasio F. A., Sills A., Glebbeek E., 2013, ApJ, 777, 106
  * Chin (1997) Chin S. A., 1997, Physics Letters A, 226, 344
  * Chin (2007a) Chin S. A., 2007a, arXiv e-prints, p. arXiv:0704.3273
  * Chin (2007b) Chin S. A., 2007b, Phys. Rev. E, 75, 036701
  * Chin & Chen (2005) Chin S. A., Chen C. R., 2005, Celestial Mechanics and Dynamical Astronomy, 91, 301
  * Correia et al. (2016) Correia A. C. M., Boué G., Laskar J., 2016, Celestial Mechanics and Dynamical Astronomy, 126, 189
  * Danby (1992) Danby J. M. A., 1992, Fundamentals of celestial mechanics. Willmann-Bell, Richmond, Va., U.S.
  * Dehnen (2014) Dehnen W., 2014, Computational Astrophysics and Cosmology, 1, 1
  * Dehnen & Hernandez (2017) Dehnen W., Hernandez D. M., 2017, MNRAS, 465, 1201
  * Dehnen & Read (2011) Dehnen W., Read J. I., 2011, European Physical Journal Plus, 126, 55
  * Deuflhard (1983) Deuflhard P., 1983, Numerische Mathematik, 41, 399–422
  * Duchêne & Kraus (2013) Duchêne G., Kraus A., 2013, ARA&A, 51, 269
  * Eggleton & Kiseleva (1995) Eggleton P., Kiseleva L., 1995, ApJ, 455, 640
  * Einstein et al. (1938) Einstein A., Infeld L., Hoffmann B., 1938, Annals of Mathematics, 39, 65
  * Fryer et al. (2012) Fryer C. L., Belczynski K., Wiktorowicz G., Dominik M., Kalogera V., Holz D. E., 2012, ApJ, 749, 91
  * Gessner & Janka (2018) Gessner A., Janka H.-T., 2018, ApJ, 865, 61
  * Goldman & Kaper (1996) Goldman D., Kaper T. J., 1996, SIAM Journal on Numerical Analysis, 33, 349
  * Goodwin (2010) Goodwin S. P., 2010, Philosophical Transactions of the Royal Society of London Series A, 368, 851
  * Gragg (1965) Gragg W. B., 1965, SIAM Journal on Numerical Analysis, 2, 384
  * Grudić & Hopkins (2020) Grudić M. Y., Hopkins P. F., 2020, MNRAS, 495, 4306
  * Grudić et al. (2021) Grudić M. Y., Guszejnov D., Hopkins P. F., Offner S. S. R., Faucher-Giguère C.-A., 2021, MNRAS, 506, 2199
  * Hamers & Portegies Zwart (2016) Hamers A. S., Portegies Zwart S. F., 2016, MNRAS, 459, 2827
  * Hamers et al. (2021) Hamers A. S., Rantala A., Neunteufel P., Preece H., Vynatheya P., 2021, MNRAS, 502, 4479
  * Heggie (1993) Heggie D. C., 1993, arXiv e-prints, pp astro–ph/9312018
  * Heggie & Hut (2003) Heggie D., Hut P., 2003, The Gravitational Million-Body Problem: A Multidisciplinary Approach to Star Cluster Dynamics. Cambridge University Press
  * Heggie & Mathieu (1986) Heggie D. C., Mathieu R. D., 1986, in Hut P., McMillan S. L. W., eds, , Vol. 267, The Use of Supercomputers in Stellar Dynamics. Springer-Verlag, Berlin Heidelberg New York, p. 233, doi:10.1007/BFb0116419
  * Heggie et al. (2006) Heggie D. C., Trenti M., Hut P., 2006, MNRAS, 368, 677
  * Hellström & Mikkola (2010) Hellström C., Mikkola S., 2010, Celestial Mechanics and Dynamical Astronomy, 106, 143
  * Hirai et al. (2021) Hirai Y., Fujii M. S., Saitoh T. R., 2021, PASJ, 73, 1036
  * Hopcroft & Tarjan (1973) Hopcroft J., Tarjan R., 1973, Commun. ACM, 16, 372–378
  * Hu et al. (2014) Hu C.-Y., Naab T., Walch S., Moster B. P., Oser L., 2014, MNRAS, 443, 1173
  * Hulse & Taylor (1975) Hulse R. A., Taylor J. H., 1975, ApJ, 195, L51
  * Hurley et al. (2000) Hurley J. R., Pols O. R., Tout C. A., 2000, MNRAS, 315, 543
  * Hurley et al. (2002) Hurley J. R., Tout C. A., Pols O. R., 2002, MNRAS, 329, 897
  * Hurley et al. (2007) Hurley J. R., Aarseth S. J., Shara M. M., 2007, ApJ, 665, 707
  * Izzard et al. (2018) Izzard R. G., Preece H., Jofre P., Halabi G. M., Masseron T., Tout C. A., 2018, MNRAS, 473, 2984
  * Jänes et al. (2014) Jänes J., Pelupessy I., Portegies Zwart S., 2014, A&A, 570, A20
  * Kamlah et al. (2022) Kamlah A. W. H., et al., 2022, arXiv e-prints, p. arXiv:2205.04470
  * Kidder (1995) Kidder L. E., 1995, Phys. Rev. D, 52, 821
  * Kiminki & Kobulnicky (2012) Kiminki D. C., Kobulnicky H. A., 2012, ApJ, 751, 4
  * Kobulnicky et al. (2014) Kobulnicky H. A., et al., 2014, ApJS, 213, 34
  * Kochanek (1992) Kochanek C. S., 1992, ApJ, 385, 604
  * Konstantinidis & Kokkotas (2010) Konstantinidis S., Kokkotas K. D., 2010, A&A, 522, A70
  * Korol et al. (2017) Korol V., Rossi E. M., Groot P. J., Nelemans G., Toonen S., Brown A. G. A., 2017, MNRAS, 470, 1894
  * Kouwenhoven et al. (2008) Kouwenhoven M. B. N., Brown A. G. A., Goodwin S. P., Portegies Zwart S. F., Kaper L., 2008, Astronomische Nachrichten, 329, 984
  * Kremer et al. (2020) Kremer K., et al., 2020, ApJS, 247, 48
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Küpper et al. (2011) Küpper A. H. W., Maschberger T., Kroupa P., Baumgardt H., 2011, MNRAS, 417, 2300
  * Kustaanheimo & Stiefel (1965) Kustaanheimo P., Stiefel E., 1965, J. Reine Angew. Math, 218, 204
  * Lada & Lada (2003) Lada C. J., Lada E. A., 2003, ARA&A, 41, 57
  * Lahén et al. (2020) Lahén N., Naab T., Johansson P. H., Elmegreen B., Hu C.-Y., Walch S., 2020, ApJ, 904, 71
  * Li et al. (2022) Li H., Vogelsberger M., Bryan G. L., Marinacci F., Sales L. V., Torrey P., 2022, MNRAS, 514, 265
  * Lim & Rodriguez (2020) Lim H., Rodriguez C. L., 2020, Phys. Rev. D, 102, 064033
  * Longmore et al. (2014) Longmore S. N., et al., 2014, in Beuther H., Klessen R. S., Dullemond C. P., Henning T., eds, Protostars and Planets VI. University of Arizona Press, p. 291 (arXiv:1401.4175), doi:10.2458/azu_uapress_9780816531240-ch013
  * Lynden-Bell & Eggleton (1980) Lynden-Bell D., Eggleton P. P., 1980, MNRAS, 191, 483
  * Lynden-Bell & Wood (1968) Lynden-Bell D., Wood R., 1968, MNRAS, 138, 495
  * Mackey et al. (2008) Mackey A. D., Wilkinson M. I., Davies M. B., Gilmore G. F., 2008, MNRAS, 386, 65
  * Maggiore (2007) Maggiore M., 2007, Gravitational Waves. Vol. 1: Theory and Experiments. Oxford Master Series in Physics, Oxford University Press
  * Mannerkoski et al. (2019) Mannerkoski M., Johansson P. H., Pihajoki P., Rantala A., Naab T., 2019, ApJ, 887, 35
  * Mannerkoski et al. (2021) Mannerkoski M., Johansson P. H., Rantala A., Naab T., Liao S., 2021, ApJ, 912, L20
  * Mannerkoski et al. (2022) Mannerkoski M., Johansson P. H., Rantala A., Naab T., Liao S., Rawlings A., 2022, ApJ, 929, 167
  * Marchal (1990) Marchal C., 1990, The three-body problem. Elsevier
  * Mardling & Aarseth (2001) Mardling R. A., Aarseth S. J., 2001, MNRAS, 321, 398
  * McMillan (1986) McMillan S. L. W., 1986, in Hut P., McMillan S. L. W., eds, , Vol. 267, The Use of Supercomputers in Stellar Dynamics. Springer New York, p. 156, doi:10.1007/BFb0116406
  * Memmesheimer et al. (2004) Memmesheimer R.-M., Gopakumar A., Schäfer G., 2004, Phys. Rev. D, 70, 104011
  * Mikkola (1987) Mikkola S., 1987, Celestial Mechanics, 40, 329
  * Mikkola (2020) Mikkola S., 2020, Gravitational Few-Body Dynamics: A Numerical Approach.. Cambridge University Press
  * Mikkola & Aarseth (1993) Mikkola S., Aarseth S. J., 1993, Celestial Mechanics and Dynamical Astronomy, 57, 439
  * Mikkola & Aarseth (1996) Mikkola S., Aarseth S. J., 1996, Celestial Mechanics and Dynamical Astronomy, 64, 197
  * Mikkola & Aarseth (2002) Mikkola S., Aarseth S., 2002, Celestial Mechanics and Dynamical Astronomy, 84, 343
  * Mikkola & Merritt (2006) Mikkola S., Merritt D., 2006, MNRAS, 372, 219
  * Mikkola & Merritt (2008) Mikkola S., Merritt D., 2008, AJ, 135, 2398
  * Mikkola & Tanikawa (1999) Mikkola S., Tanikawa K., 1999, MNRAS, 310, 745
  * Moe & Di Stefano (2017) Moe M., Di Stefano R., 2017, ApJS, 230, 15
  * Mukherjee et al. (2021) Mukherjee D., Zhu Q., Trac H., Rodriguez C. L., 2021, ApJ, 916, 9
  * Mukherjee et al. (2022) Mukherjee D., Zhu Q., Ogiya G., Rodriguez C. L., Trac H., 2022, arXiv e-prints, p. arXiv:2205.12289
  * Naoz (2016) Naoz S., 2016, ARA&A, 54, 441
  * Naoz et al. (2013) Naoz S., Farr W. M., Lithwick Y., Rasio F. A., Teyssandier J., 2013, MNRAS, 431, 2155
  * Neumayer et al. (2020) Neumayer N., Seth A., Böker T., 2020, A&ARv, 28, 4
  * Pelupessy et al. (2012) Pelupessy F. I., Jänes J., Portegies Zwart S., 2012, New Astron., 17, 711
  * Peters (1964) Peters P. C., 1964, Phys. Rev., 136, B1224
  * Peters & Mathews (1963) Peters P. C., Mathews J., 1963, Physical Review, 131, 435
  * Piotto et al. (2004) Piotto G., et al., 2004, ApJ, 604, L109
  * Plummer (1911) Plummer H. C., 1911, MNRAS, 71, 460
  * Podsiadlowski et al. (2004) Podsiadlowski P., Langer N., Poelarends A. J. T., Rappaport S., Heger A., Pfahl E., 2004, ApJ, 612, 1044
  * Poisson & Will (2014) Poisson E., Will C. M., 2014, Gravity. Cambridge University Press
  * Portegies Zwart et al. (2010) Portegies Zwart S. F., McMillan S. L. W., Gieles M., 2010, ARA&A, 48, 431
  * Press et al. (2007) Press W., Teukolsky S., Vetterling W., Flannery B., 2007, Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press
  * Preto & Tremaine (1999) Preto M., Tremaine S., 1999, AJ, 118, 2532
  * Raghavan et al. (2010) Raghavan D., et al., 2010, ApJS, 190, 1
  * Rantala et al. (2017) Rantala A., Pihajoki P., Johansson P. H., Naab T., Lahén N., Sawala T., 2017, ApJ, 840, 53
  * Rantala et al. (2018) Rantala A., Johansson P. H., Naab T., Thomas J., Frigo M., 2018, ApJ, 864, 113
  * Rantala et al. (2020) Rantala A., Pihajoki P., Mannerkoski M., Johansson P. H., Naab T., 2020, MNRAS, 492, 4131
  * Rantala et al. (2021) Rantala A., Naab T., Springel V., 2021, MNRAS, 502, 5546
  * Rein & Spiegel (2015) Rein H., Spiegel D. S., 2015, MNRAS, 446, 1424
  * Rein & Tamayo (2015) Rein H., Tamayo D., 2015, MNRAS, 452, 376
  * Rieder et al. (2022) Rieder S., Dobbs C., Bending T., Liow K. Y., Wurster J., 2022, MNRAS, 509, 6155
  * Rizzuto et al. (2021) Rizzuto F. P., et al., 2021, MNRAS, 501, 5257
  * Rizzuto et al. (2022) Rizzuto F. P., Naab T., Spurzem R., Arca-Sedda M., Giersz M., Ostriker J. P., Banerjee S., 2022, MNRAS, 512, 884
  * Samsing et al. (2018) Samsing J., Leigh N. W. C., Trani A. A., 2018, MNRAS, 481, 5436
  * Sana et al. (2012) Sana H., et al., 2012, Science, 337, 444
  * Sandage (1953) Sandage A. R., 1953, AJ, 58, 61
  * Sheng (1989) Sheng Q., 1989, IMA Journal of Numerical Analysis, 9, 199
  * Spera et al. (2015) Spera M., Mapelli M., Bressan A., 2015, MNRAS, 451, 4086
  * Spitzer (1987) Spitzer L., 1987, Dynamical evolution of globular clusters. Princeton University Press
  * Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New Astron., 6, 79
  * Springel et al. (2021) Springel V., Pakmor R., Zier O., Reinecke M., 2021, MNRAS, 506, 2871
  * Spurzem (1999) Spurzem R., 1999, Journal of Computational and Applied Mathematics, 109, 407
  * Stone et al. (2017) Stone N. C., Küpper A. H. W., Ostriker J. P., 2017, MNRAS, 467, 4180
  * Sugimoto & Bettwieser (1983) Sugimoto D., Bettwieser E., 1983, MNRAS, 204, 19P
  * Suzuki (1991) Suzuki M., 1991, Journal of Mathematical Physics, 32, 400
  * Takahashi & Imada (1984) Takahashi M., Imada M., 1984, Journal of the Physical Society of Japan, 53, 3765
  * Thorne & Hartle (1985) Thorne K. S., Hartle J. B., 1985, Phys. Rev. D, 31, 1815
  * Tokovinin (2014a) Tokovinin A., 2014a, AJ, 147, 86
  * Tokovinin (2014b) Tokovinin A., 2014b, AJ, 147, 87
  * Toonen et al. (2016) Toonen S., Hamers A., Portegies Zwart S., 2016, Computational Astrophysics and Cosmology, 3, 6
  * Trani et al. (2019) Trani A. A., Fujii M. S., Spera M., 2019, ApJ, 875, 42
  * Valtonen & Karttunen (2006) Valtonen M., Karttunen H., 2006, The Three-Body Problem. Cambridge University Press
  * Vynatheya et al. (2022) Vynatheya P., Hamers A. S., Mardling R. A., Bellinger E. P., 2022, MNRAS, 516, 4146
  * Wall et al. (2020) Wall J. E., Mac Low M.-M., McMillan S. L. W., Klessen R. S., Portegies Zwart S., Pellegrino A., 2020, ApJ, 904, 192
  * Wang et al. (2015) Wang L., Spurzem R., Aarseth S., Nitadori K., Berczik P., Kouwenhoven M. B. N., Naab T., 2015, MNRAS, 450, 4070
  * Wang et al. (2016) Wang L., et al., 2016, MNRAS, 458, 1450
  * Wang et al. (2020a) Wang L., Nitadori K., Makino J., 2020a, MNRAS, 493, 3398
  * Wang et al. (2020b) Wang L., Iwasawa M., Nitadori K., Makino J., 2020b, MNRAS, 497, 536
  * Wang et al. (2021) Wang Y.-H., Leigh N. W. C., Liu B., Perna R., 2021, MNRAS, 505, 1053
  * Wetzstein et al. (2009) Wetzstein M., Nelson A. F., Naab T., Burkert A., 2009, ApJS, 184, 298
  * Will (2014) Will C. M., 2014, Phys. Rev. D, 89, 044043
  * Wisdom & Hernandez (2015) Wisdom J., Hernandez D. M., 2015, MNRAS, 453, 3015
  * Xu & Wu (2010) Xu J., Wu X., 2010, Research in Astronomy and Astrophysics, 10, 173
  * Yoshida (1990) Yoshida H., 1990, Physics Letters A, 150, 262
  * Zlochower & Lousto (2015) Zlochower Y., Lousto C. O., 2015, Phys. Rev. D, 92, 024022

##  Appendix A Star cluster initial conditions

###  A.1 A novel fast initial conditions generator

We setup the star cluster initial conditions (IC) for the purposes of this
study using our new initial conditions generator we term SCIC (for star
cluster initial conditions). Our IC generator code is functionally very
similar to the widely used McLuster code (Küpper et al., 2011). The main
difference is the speed of the two codes: generating a million-body star
cluster IC takes only a few seconds using SCIC. The reason for the speed of
our new code is the efficient MPI parallelization allowing the use of several
supercomputer nodes, a feature which McLuster as a shared-memory single-node
code lacks. Besides the standard $\mathcal{O}(N^{2})$ loop parallelization
also time-consuming $\mathcal{O}(N)$ loops such as single and binary stellar
evolution are parallelised in SCIC. The code is written in standard C language
with a number of Fortran libraries from the literature included for stellar
evolution. In order to summarised the basic features of the SCIC code we have
collected the main fixed and user-given parameters into Table 2. Further
details of the features of the code are given below.

###  A.2 Star cluster models

We sample the masses of individual stars using the common Kroupa (2001)
initial mass function. In the code we typically use the number of stars as the
user-given input parameter instead of the total cluster mass. This choice is
made to avoid the iterative search for the correct particle number to match
the given cluster mass when stellar evolution is used.

The spherically symmetric mass density profile of the star clusters is set
according to the Plummer (1911) profile

$\rho(r)=\frac{3M_{\mathrm{p}}}{4\pi
a_{\mathrm{p}}}\left(1+\frac{r^{2}}{a_{\mathrm{p}}^{2}}\right)^{-5/2}$ (74)

in which $M_{\mathrm{p}}$ and $a_{\mathrm{p}}$ are the mass and scale radius
of the model, respectively. The positions of single stars and binary center-
of-masses are sampled from the corresponding cumulative mass profile using the
standard inversion sampling technique as the inverse function $r(M)$ is easy
to calculate. Two random angles from a unit sphere are used to get the
Cartesian components of the position vector. The scale radius $a_{\mathrm{p}}$
is related to the half-mass radius of the model as
$r_{\mathrm{h}}=\left(2^{2/3}-1\right)^{-1/2}a_{\mathrm{p}}$. In SCIC the
half-mass radius is either user-given or obtained from the observed mass-
radius relations of young star clusters and nuclear star clusters. For lower-
mass star clusters we use the relation

$\frac{r_{\mathrm{h}}}{1\mathrm{pc}}=2.55\left(\frac{M}{10^{4}M_{\mathrm{\odot}}}\right)^{0.24}$
(75)

of Brown & Gnedin (2021). For more massive nuclear star clusters we devise a
simple ad-hoc relation

$\begin{split}\log_{\mathrm{10}}\left(\frac{R_{\mathrm{e}}}{1\mathrm{pc}}\right)&=\frac{1}{2},&M<10^{6.5}M_{\mathrm{\odot}}\\\
\log_{\mathrm{10}}\left(\frac{R_{\mathrm{e}}}{1\mathrm{pc}}\right)&=\frac{1}{2}\log_{\mathrm{10}}\left(\frac{M}{1M_{\mathrm{\odot}}}\right)-\frac{11}{4},&M\geq
10^{6.5}M_{\mathrm{\odot}}\end{split}$ (76)

with $r_{\mathrm{h}}=((1/2)^{-2/3}-1)^{-1}R_{\mathrm{e}}$ for the Plummer
(1911) model. The formula phenomenologically captures the observed flat radius
distribution at smaller masses and the power-law behaviour above the threshold
mass of $M\sim 10^{6.5}M_{\mathrm{\odot}}$ (see e.g. Fig. 7 of Neumayer et al.
2020).

The space velocities of single stars and binary center-of-masses are obtained
using the Plummer (1911) potential

$\phi(r)=-\frac{GM}{a_{\mathrm{p}}}\left(1+\frac{r^{2}}{a_{\mathrm{p}}^{2}}\right)^{-1/2}$
(77)

and the energy distribution function

$f(\mathcal{E})=\begin{cases}\frac{24\sqrt{2}a_{\mathrm{p}}^{2}}{7\pi^{3}G^{5}M_{\mathrm{p}}^{5}}(\mathcal{E})^{7/2}&\hskip
28.45274pt\mathcal{E}\geq 0\\\ 0&\hskip 28.45274pt\mathcal{E}<0\end{cases}$
(78)

using $\mathcal{E}=-E=-\frac{1}{2}m\|\boldsymbol{v}\|^{2}-\phi(r)$ together
with the von Neumann rejection sampling method (see e.g. the Art of
Computational Science online material222www.artcompsci.org by P. Hut & J.
Makino and Heggie & Hut 2003 for extensive practical details of the
implementation). As with the positions the isotropic components of the
velocity vectors are obtained by sampling two angles uniformly from the unit
sphere.

Table 2: The main initial conditions generator SCIC parameters, their symbols and their typical values as described in the text. Parameter | Symbol | Details
---|---|---
particle number | $N$ | $\lesssim 10^{8}$
IMF | $\xi(m)$ | Kroupa (2001)
half-mass radius | $r_{\mathrm{h}}$ | user-given or observed
binary fraction | $f_{\mathrm{bin}}$ | 0 $\leq f_{\mathrm{bin}}\leq 1.0$
semi-major axis distribution | $f(a)$ | flat in log space
min, max semi-major axis |  $a_{\mathrm{min}}$, $a_{\mathrm{max}}$ | $a_{\mathrm{min}}\leq a_{\mathrm{max}}$
eccentricity distribution | $f(e)$ | thermal
binary pairing mass limit | $m_{\mathrm{pair}}$ | $5M_{\mathrm{\odot}}$
mass ratio distribution | $f(q)$ | const., $0.1\leq q\leq 1.0$
stellar population age | $t_{\star}$ |  $<$ age of the Universe
stellar population metallicity | $Z_{\star}$ | $10^{-4}<Z\star<0.03$
artificial binary breaking limit | $r_{\mathrm{break}}$ |  $\sim$ a few times $a_{\mathrm{max}}$

###  A.3 Binary stars

Defining the number of single stars and binary systems in a star cluster as
$N_{\mathrm{s}}$ and $N_{\mathrm{b}}$ respectively the total number of stars
is $N=N_{\mathrm{s}}+2N_{\mathrm{b}}$. The number of binary systems can also
be expressed using the binary fraction $f_{\mathrm{bin}}$ of the cluster as
$N_{\mathrm{b}}=1/2f_{\mathrm{bin}}N$. We use the binary fraction as an user-
given input parameter defining the binary content of the initial conditions.

The initial binary component masses are paired above
$m_{\mathrm{pair}}=5M_{\mathrm{\odot}}$ using a flat mass ratio distribution
$q=m_{\mathrm{2}}/m_{\mathrm{1}}\propto\text{const}$ with $0.1\leq q<1.0$
(Kiminki & Kobulnicky 2012; Sana et al. 2012; Kobulnicky et al. 2014). In the
practical implementation beginning from the most massive star we assign
secondary star masses by finding the actual closest mass to the sampled random
secondary mass from the remaining stars until no primary stars above
$5M_{\mathrm{\odot}}$ remain.

For binary stars the initial semi-major axis distribution is flat in
logarithmic space between the user given minimum and maximum values of
$a_{\mathrm{min}}$ and $a_{\mathrm{max}}$. For reasonable cluster models
$a_{\mathrm{max}}\ll r_{\mathrm{h}}$. The binary eccentricity distribution is
assumed to be thermal initially, i.e $f(e)=2e$.

The age and metallicity of the stellar population are user-given initial
parameters. At the moment our code supports a single stellar population per
star cluster. In the case of an evolved stellar population the properties of
the single (e.g. type, mass, radius) and binary stars (e.g. types, masses,
radii, semi-major axes, eccentricities) of the given age and metallicity are
computed using the widely-used SSE (Hurley et al., 2000) and BSE (Hurley et
al., 2002) stellar evolution libraries. The SCIC calls the libraries in
parallel (a single task per binary) through a simple C-to-Fortran code
interface. Due to the parallel approach SCIC is faster than McLuster in
stellar evolution as the latter code only has a serial library interface.

We allow the binary stars to merge during binary stellar evolution using the
BSE library merging criteria (Hurley et al., 2002). In addition to the BSE
merger criteria we merge all systems which orbit closer than the sum of the
radii (or sum of ISCO radii) of the components of the binary. We also check
the systems with a single compact object after stellar evolution for tidal
disruption events just as in the BIFROST code described in Section 3.14. The
merger remnant properties are obtained using the BSE package. Depending on the
types and masses of the merging objects there may or may not to be a remnant.
Finally we check the binary systems containing a single or two compact objects
for gravitational-wave driven inspirals within the next $0.01$ Myr again as in
Section 3.14. In our test simulations with BIFROST initial conditions
generated using McLuster usually lead in up to a few tens of tidal disruption
mergers within the first few time-steps immediately in the beginning of the
BIFROST simulation. For initial conditions generated with SCIC this somewhat
spurious burst of merger events does not occur due to a more detailed set of
merger criteria in the IC generator code.

After stellar evolution the relative positions and velocities of the binary
star components are added to the position and velocity of the center-of-mass
of the system. While the binary semi-major axes and eccentricities are given
by the stellar evolution libraries we assign random values to the other
orbital elements (inclination $i$, longitude of the ascending node $\Omega$,
argument of periapsis $\omega$ and mean anomaly $M$). For obtaining the
eccentric anomaly required for the positions and velocities we use the same
Kepler solver as in the BIFROST code described in Section 3.11.

###  A.4 Massive black holes

We have included an option to add a central black hole in the star cluster
initial conditions in our code. The user-given input parameters for the
massive black hole particle are the mass and the magnitude of spin vector $s$.
Typically we conservatively set $s$ to 10% of the maximal black hole spin. For
a single central black hole the initial orientation of the spin vector does
not matter as the cluster model is spherically symmetric anyway so we decide
to align it with the z-axis of the coordinate system.

###  A.5 Avoiding artificial binaries and multiplets

Even when setting the binary fraction $f_{\mathrm{bin}}=0\%$ a star cluster
initial condition typically contains a few binary systems formed by random
chance which we refer to here as artificial binaries. The number
$N_{\mathrm{art}}$ of these systems increases with increasing cluster particle
number $N$ and decreasing half-mass radius $r_{\mathrm{h}}$ which both
increase the stellar number density of the cluster.

A binary system which survives in a star cluster for a long time needs to be
sufficiently tight (e.g. Heggie & Hut 2003), typically on milliparsec-scale
separations or less. This separation scale is considerably smaller than the
half-mass radius of typical star clusters (of the order of a few pc) or mean
inter-stellar distance in the cluster centers (of the order of $0.01$ pc to
$0.1$ pc). Still, the sampling of the initial star particle positions
independently of each other will inevitable create close-by particle pairs. An
artificial binary is formed if the relative velocity of the two particles is
small enough for the binary to be bound. The artificial binaries typically
have high eccentricities originating from their required small apocenter
velocities.

Artificial binaries formed in the initial conditions generator are not a
problem in itself but rather complicate the comparison of simulation runs
using same star cluster model properties with different random seeds,
especially when $f_{\mathrm{bin}}$ is low and $N$ is small. As artificial
binaries are most very eccentric they are thus computationally expensive to
integrate compared to typical binary systems. This may have an effect on for
example timing and scaling tests of a simulation code with a set of initial
conditions with a low or zero $f_{\mathrm{bin}}$ as initial conditions with
tight eccentric binaries are the slowest to integrate in wall-clock time. We
emphasised that for fair and controlled timing and scaling tests initial
conditions containing artificial binaries should be avoided. Thus, a robust
way beyond simply removing artificial binaries by hand should be devised.

We also note that in the case of $f_{\mathrm{bin}}>0$ (typically short-lived)
artificial multiplets of three or more stars will also occur. As these
artificial multiple systems are rarely stable they will quickly dissolve in a
strong few-body interaction (e.g. (Valtonen & Karttunen, 2006)). This may lead
to spurious formation of high-velocity stars or artificially enhanced merger
rates. We however note in addition that including physical multiple systems
(as motivated by observations) is an interesting prospect for future work.

A simple method to avoid artificial binaries and multiplets in the initial
conditions generation is to impose a minimum distance $r_{\mathrm{break}}\sim$
a few times $a_{\mathrm{max}}$ between particles which are not members of the
same binary system. We search for the shortest inter-particle distances
between the $N_{\mathrm{s}}+N_{\mathrm{b}}$ single stars and binary systems.
If an artificial binary is detected (i.e. $a<r_{\mathrm{break}}$) the
positions and velocities of the particles in question are re-sampled just as
they were sampled for the first time: the resulting new realisation of the
cluster model is as realistic a cluster model as the original one, only the
random numbers differ. The process may need to be repeated up to $\sim 10$
times to get rid of all the artificial binaries. The recipe also holds for the
binary center-of-masses as for single stars as the binaries treated as point
particles before assigning the positions and velocities of the binary
components in the final phase of the IC generation.

However, differences in the partially re-sampled random realisations begin to
unfortunately arise when either $a_{\mathrm{max}}$ is large or the cluster is
very dense: re-sampled particles from the centre of the cluster will end up
more probably to larger radii and the density profile of the cluster starts to
deviate from the desired one and the $f_{\mathrm{bin}}$ begins to develop a
radial gradient. A solution for this issue is to only re-sample the angular
coordinates (both for positions and velocities) leaving the original radial
position $\|\boldsymbol{r}_{\mathrm{i}}\|$ and space velocity
$\|\boldsymbol{r}_{\mathrm{i}}\|$ unchanged. A practical algorithm in our SCIC
code first performs the angular re-sampling method for the artificial binaries
using up to $10$ iterations before switching to the original re-sampling if
artificial systems remain. After another $5$-$10$ iteration rounds typically
only few or zero artificial binaries remain in the cluster model.


# The First Law for the Kerr-NUT Spacetime

Ernesto Frodden<EMAIL_ADDRESS>Departamento de Física, Universidad de
Santiago de Chile, Avenida Ecuador 3493, Santiago, Chile.  Diego Hidalgo
<EMAIL_ADDRESS>Centro de Estudios Científicos (CECs), Av. Arturo Prat
514, Valdivia, Chile
Departamento de Física, Universidad de Concepción, Casilla 160-C, Concepción,
Chile,  Instituto de Ciencias Físicas y Matemáticas, Universidad Austral de
Chile, Casilla 567, Valdivia, Chile.

###### Abstract

We show that the first law for the Kerr-NUT is straightforwardly established
with the surface charge method. The entropy is explicitly found as a charge
and its value is not proportional to the horizon area. We conclude that there
are unavoidable contributions from the Misner strings to the charges, still
the mass and angular momentum gets standard values. However, there are no
independent charges associated with the Misner strings.

Surface Charges, Misner Strings

##  I Introduction

Taub-NUT (TN) spacetime is a vacuum solution of General Relativity Taub ; NUT
. For a long time, the physical interpretation of TN spacetime has been
elusive. This spacetime contains unusual properties that make it attractive
for study and constitute a typical counter-example for weakly established
statements in the field (think about the area law in black hole physics or
standardized asymptotic behavior for spacetimes). The characterizing feature
of TN spacetime is that the additional parameter, the NUT parameter $n$,
causes string-like defects Misner:1963fr , nowadays called Misner strings.
These strings are interpreted as the gravitational analogue of the Dirac
strings appearing for magnetic monopoles in the electromagnetic context. Thus,
the NUT property is further interpreted as a gravitomagnetic charge, i.e. as a
charge dual to the mass in the same vein of electromagnetic duality.

There are at least two remarkable interpretations of the TN solution. The
first one, presented by Misner Misner:1963fr , argued that the string-like
defect could be made invisible if the time coordinate is made periodic. The
second interpretation, employed by Bonnor bonnor1969 , keeps the strings as a
certain “massless source of angular momentum.”

Most of the studies on TN spacetimes is done in the Euclidean regime
Chamblin:1998pz ; Ghezelbash:2008zz ; Emparan:1999pm ; Mann:1999pc ;
Mann:1999bt ; Garfinkle:2000ms ; Johnson:2014xza ; Johnson:2014pwa ;
Clarkson:2002uj ; Astefanesei:2004kn ; Kalamakis:2020aaj ; Ciambelli:2020qny ;
Arratia:2020hoy ; Hawking:1998jf ; Hawking:1998ct  because gravitational
instantons Hawking:1976jb  might play a role in the path integral formulation
of gravity. However, the Lorentzian version of TN spacetime was critically
revisited Clement:2015cxa ; Clement:2015aka . The authors considered the
geodesics motion in the presence of Misner strings. The main result is that TN
spacetime turns out to be geodesically complete: i.e Free-falling observers do
not see strings as they cross them. Guided by this last analysis, renewed
interest in TN spacetime emerged in the community as there is a window for the
TN solutions to have physical relevance Kubiznak:2019yiu ; Gonzalez:2017sfq ;
Cano:2021qzp ; BallonBordo:2019vrn ; Bordo:2019slw ; Durka:2019ajz ;
Bordo:2019rhu ; Bordo:2019tyh ; Zhang:2021pvx .

In this work, we pursue a consistent formulation of thermodynamics for the
Kerr-Nut spacetime, also known as rotating TN spacetime in the Lorentzian
regime. This solution is also referred to as the Kerr NUT solution because it
is a generalization of the Kerr black hole solution including the NUT
parameter. The final hope is to settle down a thermodynamics interpretation
for the TN spacetime to increase the knowledge about this maybe-physical
solution and may help to have a new arena where quantum gravity ideas might be
explored. As we will see, the spacetime entropy for this solution is different
to the usual black holes.

Recently, Kubiznak:2019yiu , it has been claimed that the first law for
Lorentzian rotating TN spacetimes, Kerr-NUT, can be formulated assuming: 1) An
entropy given by the one-four of the area of the horizon (area law), and 2)
the existence of independently NUT parameters that control the effect of
Misner strings at the first law.

In this work, we review how the Misner strings enter in the formulation of the
first law. Using the surface charge method Frodden:2019ylc . The TN charges:
Mass and angular momentum, are derived with a suitable spacetime integration.
With these conserved quantities, the variation of the entropy is naturally
obtained as a charge too. This procedure, plus integrability considerations,
lead us to restrict the parameter space and find an entropy that differs from
the area law, in accordance with Hawking:1998jf .

Our paper is organized as follows: In Section II, we briefly review basic
facts about the Kerr-NUT. The surface charge method improved to deal with
string-like defects is presented in Section III. Section IV is devoted to
constructing the first law of thermodynamic for TN solution and find the
appropriate entropy that fits on it. As a consistency check of the results, in
Subsection IV.1, we show that all exposed quantities and the first law can be
framed in the grand canonical ensemble. Summary and final discussion are
presented in Section V. For completeness we include an appendix exhibiting
intermediate steps of the computations.

##  II Kerr-NUT Solution

The Kerr-NUT, or rotating TN spacetime, generalizes the Kerr metric, with an
extra NUT parameter $n$ as

$\displaystyle
ds^{2}=-\frac{\Delta}{\Sigma}\left(dt+(2n\cos\theta+2Cn-a\sin^{2}\theta)d\phi\right)^{2}+\frac{\Sigma}{\Delta}dr^{2}$
$\displaystyle+\,\frac{\sin^{2}\theta}{\Sigma}\left(adt-(r^{2}+a^{2}+n^{2}-2anC)d\phi\right)^{2}+\Sigma\,d\theta^{2}\,,$
(1)

with

$\displaystyle\Delta$ $\displaystyle=$ $\displaystyle
r^{2}+a^{2}-2mr-n^{2}\,,$ (2) $\displaystyle\Sigma$ $\displaystyle=$
$\displaystyle r^{2}+(n+a\cos\theta)^{2}\,,$ (3)

where $m$ is the integration constant associated with the mass, $n$ the NUT
parameter, and $a$ the rotation parameter. The solution with $C=0$ is a
particular Demianski-Newman vacuum spacetime metriccite . The parameter $C=\pm
1,0$ was first introduced in Manko:2006bf . The solution presents Misner
strings as $(\nabla t)^{2}$ becomes ill-defined along the axis Misner:1963fr .
The discrete parameter $C$ specifies different configurations for the strings
Manko:2006bf : For $C=+1$ one string is at the north pole axis while the south
pole axis is regular, for $C=-1$ is the opposite situation, and for $C=0$
there are two strings symmetrically located.

The solution has a Killing horizon BallonBordo:2019vrn ; podolskybook  located
at the largest root of $\Delta=0$ and it is generated by the Killing vector
$\xi=\partial_{t}+\Omega\,\partial_{\phi}$, with $\Omega$ the angular velocity
of a ZAMO (zero angular momentum observer) at the horizon as seen by a static
observer at infinity.

##  III Surface charges for spacetimes with strings

To characterize TN spacetime we can compute charges associated with
isometries. General Relativity is a gauge theory and a naive use of the First
Noether Theorem produces misleading charges. An improved method to extract
charges is required. We use the standard surface charge method Lee:1990nz ;
Iyer:1994ys ; Wald:1999wa ; for isometries it is equivalent to the symplectic
method explained in Barnich:2001jy  (both methods were revisited in
Frodden:2019ylc ).

An exact symmetry (here equivalent to a Killing symmetry) that is generated by
a Killing vector field $\xi$ produces an on-shell conservation law equation:
$d\boldsymbol{k}_{\xi}\approx 0$. The integration of this equation over a
spacetime volume $\Sigma$, with two disconnected boundaries
$\partial\Sigma=S\cup S^{\prime}$, defines a conserved surface charge
variation as

$\not{\delta}Q_{\xi}=\oint_{S}\boldsymbol{k}_{\xi}\,,$ (4)

the conservation is understood as the fact that for any other closed two-
surface $S^{\prime}$, obtained from a continuous deformation of $S$; the
previous expression produces the same $\not{\delta}Q_{\xi}$. In particular,
this is true for deformations in time and radial directions. In’ (4) we use
bold notation to indicate spacetime differential forms, namely the surface
charge density is the two-form
$\boldsymbol{k}_{\xi}=\frac{1}{4}k^{\mu\nu}_{\xi}\varepsilon_{\mu\nu\rho\sigma}dx^{\rho}\wedge
dx^{\sigma}$. See Frodden:2019ylc  for more details and the explicit formula
of $k^{\mu\nu}_{\xi}$ in terms of a general metric or tetrad variables.

Here we stress two aspects of (4). First, the slashed delta symbol in
$\not{\delta}Q_{\xi}$ means we compute a one-form in phase space which is not
necessarily exact. Exactness in phase space, namely
$\delta(\not{\delta}Q_{\xi})=\delta^{2}Q_{\xi}=0$, is an extra condition
required for the existence of the charge $Q_{\xi}$, in the usual parlance this
is the integrability condition. A second aspect is that integration is over
any closed smooth surface $S$ where the conservation law holds (not
necessarily at infinity). As we check in a moment, Misner strings violate this
requirement, and the integration region needs to be corrected accordingly.

As in Kerr spacetime, for the Kerr-NUT, there are also two exact symmetries
present: Time translations generated by $\partial_{t}$ and axial rotations
generated by $\partial_{\phi}$. A direct computation of (4) fails as it
produces charges dependent on the surface $S$ (e.g. dependent on the radius of
a given sphere): Conservation is broken! This is a problem the reader may be
tempted to solve by taking a sphere at infinity. But, as explained before,
conservation here means precisely that one can compute charges not necessarily
at the asymptotic region. The solution to the problem comes by excluding the
region with the Misner strings. More precisely, we integrate the conservation
law equation in a volume without the strings. As a consequence, the result
does not decouple into the usual two integrals, one over each closed surface
$S^{\prime}$ and $S$ (which is what makes of (4) a seed to compute a charge).
Instead, extra contributions to (4) emerge.

Figure 1: Integration region: A three-volume, $\Sigma$, delimited by two
spheres and two cones. The blue lines stand for the Misner strings
symmetrically positioned (case $C=0$).

Let us consider the surface charge density $\boldsymbol{k}_{\xi}$, with $\xi$
a Killing vector field, either $\partial_{t}$ or $\partial_{\phi}$. We
integrate the three-form $d\boldsymbol{k}_{\xi}$ on a volume
$\Sigma_{\epsilon}$ limited by two spherical shells, say $S_{1}$ and $S_{2}$,
with radius $r_{1}<r_{2}$, and two tubes $T_{S}$ and $T_{N}$ (actually cones)
excluding the strings from the volume, see Figure 1. The south cone $T_{S}$
has a small opening angle, $\theta=\pi-\epsilon$, towards the south pole, and
the north cone $T_{N}$ has a small opening angle, $\theta=0+\epsilon$, towards
the north pole. We integrate $d\boldsymbol{k}_{\xi}\approx 0$ over
$\Sigma_{\epsilon}$, thus

$0=\int_{\Sigma_{\epsilon}}\\!\\!d\boldsymbol{k}_{\xi}=\int_{\partial\Sigma_{\epsilon}}\\!\\!\boldsymbol{k}_{\xi}=\int_{S_{1}}\\!\\!\boldsymbol{k}_{\xi}-\int_{S_{2}}\\!\\!\boldsymbol{k}_{\xi}+\int_{T_{S}}\\!\\!\\!\boldsymbol{k}_{\xi}+\int_{T_{N}}\\!\\!\\!\boldsymbol{k}_{\xi}\,.$
(5)

When there are no strings, the third and fourth term vanishes in the shrinking
limit $\epsilon\\!\to\\!0$ and we recover the usual expression for the charge
(4). However, for TN spacetime, strings contribute in the limit of a vanishing
$\epsilon$.

In general, the radial integration over each tube decouples into two terms
associated with the circles found as the intersection of the tube with the
shells

$\displaystyle\int_{T_{S}}\boldsymbol{k}_{\xi}$ $\displaystyle=$
$\displaystyle
K_{\xi}\big{|}_{r=r_{1};\,\theta=\pi-\epsilon}-K_{\xi}\big{|}_{r=r_{2};\,\theta=\pi-\epsilon}\,,$
(6a) $\displaystyle\int_{T_{N}}\boldsymbol{k}_{\xi}$ $\displaystyle=$
$\displaystyle
K_{\xi}\big{|}_{r=r_{1};\,\theta=\epsilon}-K_{\xi}\big{|}_{r=r_{2};\,\theta=\epsilon}\,.$
(6b)

Then, we take the limit $\epsilon\to 0$ and reorder (5) with terms associated
with each shell to get an expression for the surface charge that is defined
exclusively at a sphere $S=S_{1}$. This quantity is truly conserved, i.e. its
value does not depend on the specific surface $S$

$\not{\delta}Q_{\xi}=\oint_{S}\boldsymbol{k}_{\xi}+K_{\xi}\big{|}_{south}+K_{\xi}\big{|}_{north}\,,$
(7)

where we refer to $\theta=\pi$ as south, and $\theta=0$ as north. Although we
formally use a specific parametrization to refer to geometrical objects, the
result does not depend on them. Before shrinking the latter around the
strings, we can deform either the shells or the tubes and still have the same
resulting surface charge. In particular, (7) does not depend on the specific
radius we compute it.

Therefore, we have shown how the $K$ contributions solve the problem of a
naive computation with (4) when string-like defects cross the surface $S$. For
spacetimes without strings the $K$ terms automatically vanish. Still, we note
that we can introduce artificial string-like defects by changing coordinates
even for regular spacetimes. This can break down the conservation of (4)
unless terms similar to the $K$ are taken into account. An example of this is
the defect appearing in the axis for some versions of the Kerr-(a)dS metric
solution which can be avoided by a $\phi$ escalation.

We emphasize that to get (7), we rely on the correct integration of the local
conservation law. This requires identifying the tubes to exclude the defects:
The regions where $d\boldsymbol{k}_{\xi}\approx 0$ does not hold.

###  III.1 Mass and Angular Momentum

The two Killing symmetries produce two independent varied surface charges. We
check that our corrected formula delivers the expected result.

The solution space is spanned by three parameters $m$, $a$, and $n$. For time
translations exact symmetry, we expect the mass/energy, thus we rename $\delta
M\equiv\not{\delta}Q_{\partial_{t}}$, and (7) produces (for details check the
appendix)

$\delta M\,=\,\delta m\,,$ (8)

which is of course integrable as $\delta^{2}M=0$. We set the integration
constant equal to zero to obtain the finite charge $M=m$, which means that for
Kerr-NUT spacetime the parameter $m$ is still identified with the total
mass/energy of spacetime.

For the axial symmetry generator we expect the angular momentum, thus we
rename $\delta J\equiv\not{\delta}Q_{-\partial_{\phi}}$, and we get

$\delta J\,=\,(a-3Cn)\delta m+m\delta a-3Cm\delta n\,.$ (9)

which again is an integrable expression, $\delta^{2}J=0$, and produces the
angular momentum $J=(a-3Cn)m$. Note that, as expected, the TN spacetime with
$a=0$ still carries angular momentum but only when strings are not
symmetrically distributed along the rotation axis.

Let us conclude here with two remarks. First, the charges do not present
divergences (e.g. $r$-divergences when evaluated in asymptotic regions). While
each integral in (7) is $r$-dependent, therefore susceptible to divergences,
it is only the sum of the three that produces $r$-independent charges. Then,
the $K$-es are essentials and the Misner strings contribute to the charges
though them. The second remark is that no asymptotic region is required. We do
not demand a particular location for the sphere $S$. In that sense, the
charges computed above are quasi-local.

##  IV Entropy and The First Law

In analogy with black holes thermodynamics, we now treat the Kerr-NUT
spacetime as a thermodynamic system. We constructively present the first law
to determine the entropy of this spacetime following the black hole
prescription, i.e. we introduce a temperature and an angular momentum related
to the surface $r=r_{+}$, the larger root of $\Delta=0$. For later
convenience, we focus on the case with $C=0$.

Let us introduce a temperature (proportional to the surface gravity at $r_{+}$
redshifted to infinity)

$T=\frac{\sqrt{m^{2}-a^{2}+n^{2}}}{4\pi\left(m\left(m+\sqrt{m^{2}-a^{2}+n^{2}}\right)+n^{2}\right)}.\,$
(10)

The angular velocity, of the surface $r=r_{+}$ as measured from infinity, is

$\Omega=\left.\frac{g_{t\phi}}{g_{\phi\phi}}\right|_{r_{+}}=-\frac{a}{2\left(m\left(m+\sqrt{m^{2}-a^{2}+n^{2}}\right)+n^{2}\right)}\,.$
(11)

The expectation is that in the first law the temperature multiplies the
entropy variation, we call it $\not{\delta}S$, and the angular velocity
multiplies the angular momentum variation (9). Then, we can isolate
$\not{\delta}S$ from the first law to get

$\not{\delta}S=\frac{1}{T}\left(\delta M-\Omega\delta
J\right)=\frac{2\pi}{\sqrt{m^{2}-a^{2}+n^{2}}}\times\\\
\left(\big{(}2m(m+\sqrt{m^{2}-a^{2}+n^{2}})-a^{2}+2n^{2}\big{)}\delta
m-ma\delta a\right)\,.$ (12)

Note that for $n=0$ the previous expression is integrable, and as expected, it
is the variation of the usual one-fourth of the area black hole entropy of
Kerr spacetime. Even for non-rotating TN spacetime, $n\neq 0$, the expression
is not integrable anymore, $\delta^{2}S\neq 0$, which is a problem to have a
well-defined entropy charge.

Still the entropy variation is conserved in the sense that it can be directly
computed from (7) with the combination of Killing vectors
$\xi=\frac{1}{T}\left(\partial_{t}+\Omega\,\partial_{\phi}\right)$. That means
it is a proper varied surface charge as it does not depend on the surface one
chooses to compute it. Therefore, this expression should be a good seed to
compute an entropy as a Noether charge as we already did with the mass or the
angular momentum (we keep the Noether charge name, understanding that surface
charge method to compute charges generalizes the Noether procedure
Frodden:2019ylc ).

In the literature, it has been explored the possibility of adding extra terms
to the first law associated with the Misner strings. At this level, this
strategy is equivalent to the splitting of (12) in different terms and
providing interpretation to each one. Here we choose another path. We use the
surface charge method as a guide and keep the possibility of having an entropy
that is indeed conserved: A Noether charge. Therefore the expression in (12)
must not be split; otherwise, we can not write the entropy variation as (7)
with an arbitrary closed surface.

The alternative path: We have checked the suggestive alternative of splitting
(12) into surface and tubes terms associated to the Misner strings as in (7).
As explained before, separately surfaces and tubes terms are not conserved
quantities and neither they are integrable. One may be tempted to fix the
closed surface at $r=r_{+}$. In this case, the surface integral in (7) indeed
becomes the variation of one-fourth of the area. However, we are not computing
a conserved charge if we have to fix the surface where it is defined.
Furthermore, insisting in this approach, the terms associated with the tubes
are not integrable, and another temperatures should be introduced to achieve
their integrability. This path would be similar to the one presented in
Bordo:2019rhu  (still, we have a different angular momentum stemming from the
other methods used to compute charges).

A way to have, simultaneously, the entropy as conserved and integrable is to
reduce the parameter space. Let us assume that $m$, $a$, and $n$ are dependent
quantities and solve this dependence such that the integrability condition is
satisfied. Consider $m=m(n,a)$. Because of the way these parameters enter in
the metric components, it is easy to see that all the tree parameters scale
with the same factor under scaling of the metric. Then, the allowed relation
among the parameters is simpler: $m=n\,h(a/n)$, with $h(a/n)$ a function to be
determined by the integrability condition $\delta(\not{\delta}S)=0$. This
produces a differential equation solved simply by $h(a/n)=\alpha$, with
$\alpha$ an integration constant. Thus

$m=\alpha\,n\,,$ (13)

ensures integrability of the entropy, and it reduces the parameter space of
our solution from three to two parameters. Replacing this in the equation for
the varied charge (12), we can integrate it to obtain an entropy

$S=2\pi\alpha\,n^{2}\left(\alpha+\sqrt{1+\alpha^{2}-\frac{a^{2}}{n^{2}}}\right)\,,$
(14)

where we have set to zero the integration constant. Our method does not fix
$\alpha$ and in what follows we keep it free.111We can fix $\alpha$ if we
require the entropy to be compatible with the non-rotating TN entropy obtained
through other methods. For instance, Hawking and Hunter in Hawking:1998jf
found for the Euclidean TN-bolt counterpart $S_{HH}=\pi n^{2}$, this fixes
$\alpha=\frac{1}{2\sqrt{2}}$. This entropy fulfills a standard first law

$\delta M=T\delta S+\Omega\delta J\,,$ (15)

where $M=\alpha\,n$ and $J=\alpha\,a\,n$. In contrast with the usual black
hole solutions, the entropy is not one-fourth of the area, first noted in
Hawking:1998jf  for non-rotating TN and using Euclidean methods. The first law
presented is minimal because it does not require additional terms not
justified in terms of charges. It may be seen as non-standard as we had
reduced the parameter space to achieve integrability. But as explained before,
this is required to have $M$, $J$, and $S$ as well-defined charges to compose
the first law. Therefore, we conclude that there is a consistent thermodynamic
description only for a subset of parameter space of the Kerr-NUT spacetime. If
TN spacetimes will ever have a physical meaning this is a prediction of our
approach.

###  IV.1 Grand Canonical Ensemble

Our result is compatible with the Euclidean action to describe the
thermodynamics. To show this, we use the grand canonical ensemble function,
$G$, already computed in the Appendix of Bordo:2019rhu ; Emparan:1999pm  as
the flat limit of a regularized Euclidean (a)dS action divided by $\beta=1/T$.
The result is simply $G=m/2$. We express it as a function of the intensive
thermodynamic variables at play to get

$G(T,\Omega)=\frac{1}{8\pi
T\left(1+\alpha^{\prime}\sqrt{1+\displaystyle\frac{\Omega^{2}}{4\pi^{2}T^{2}}}\right)}\,.$
(16)

With $\alpha^{\prime}=\sqrt{1+1/\alpha^{2}}$. The expression is consequence of
a simple algebraic manipulation of $T(a,n)$, $\Omega(a,n)$, and $m=\alpha\,n$.
Note its simplicity. Now, the first law for the Kerr-NUT spacetimes can be
reproduced in the standard way. As expected, the entropy (14), the angular
momentum, and the mass/energy are: $S=-\frac{\partial G}{\partial T}$,
$J=\frac{\partial G}{\partial\Omega}$ and $M=G-J\Omega+ST$, respectively.

##  V Conclusions and Outlook

In this note, we reported a procedure to use the covariant surface charge
method for spacetimes with string-like defects and immediately applied it to
the Kerr-NUT spacetime. The critical step was to note that the equation,
$d{\boldsymbol{k}}_{\xi}\approx 0$, is not well behaved over the Misner
strings; therefore, the three-volume where it is integrated can not contain
them. The procedure is general and can be easily applied in particular for
accelerating black holes.

For the Kerr-NUT spacetime, besides the explicit computation of the natural
charges, mass and angular momentum, the thermodynamic picture was completed by
writing down the first law and computing the entropy as an integrated charge.
We achieve the integrability of the entropy without introducing additional
terms in the first law. As a general rule, the introduction of new terms in
the first law should be done with caution, we try to avoid it. Instead, we
impose a reduction in the parameter space of the solution. Thus, the principle
that charges appearing in the first law ($M$, $J$, and $S$) are Noether
charges, all at the same footing, is respected. The entropy found is not
proportional to the area, which is something already expected for TN
spacetimes.

Let us finish with a comment on our attempt to extend the results to the (a)dS
Kerr-NUT family. The strategy to compute mass and angular momentum charges is
equally successful there, however we could not solve the differential equation
for the integrability of the entropy. Then, we leave for a future project to
extend the thermodynamic picture for that spacetime family or in full
generality for the Plebanski-Demianski family Griffiths:2005qp .

###### Acknowledgements.

The authors would like thank to Remigiusz Durka for a valuable discussion on
the topic. The Centro de Estudios Científicos (CECs) is funded by the Chilean
Government through the Centers of Excellence Base Financing Program of ANID.

##  Appendix A Charge Computations

The surface charge for GR is given by Lee:1990nz ; Iyer:1994ys ; Wald:1999wa ;
Barnich:2001jy ; Frodden:2019ylc

$\displaystyle k_{\xi}^{\mu\nu}$ $\displaystyle=$
$\displaystyle\sqrt{-g}\kappa\left(\xi^{[\nu}\nabla_{\sigma}\delta
g^{\mu]\sigma}-\xi^{[\nu}\nabla^{\mu]}\delta g+\xi_{\sigma}\nabla^{[\mu}\delta
g^{\nu]\sigma}\right.$ (17)
$\displaystyle~{}~{}~{}~{}~{}~{}~{}~{}~{}\left.-\frac{1}{2}\delta
g\nabla^{[\nu}\xi^{\mu]}+\delta
g^{\sigma[\nu}\nabla_{\sigma}\xi^{\mu]}\right)\,,$

where $\kappa=c^{4}/(8\pi G)$, $c$ the speed of light, $G$ the Newton’s
constant, $g$ stands for the determinant of the metric $g_{\mu\nu}$, and
$\delta g\equiv g_{\mu\nu}\delta g^{\mu\nu}$.222We work in units where
$c=G=1$. Notice that the last three terms of the equation corresponds to the
variation of the Komar’s integrand commonly used to compute charges.

The integral of $k_{\xi}^{\mu\nu}$ over a closed surface splits in four
regions: Two concentric spheres with holes and two cones (see Fig. 1). Over a
sphere at constant $r$ the required component is $k^{tr}_{\xi}$ while for the
cones it is $k^{\theta t}_{\xi}$. The explicit expressions for them in the
Kerr-Taub-NUT spacetime are quite large, thus, for the sake of space with do
not present them here (still it is a straightforward computation using a
software). To give a taste we exhibit just the coefficient accompanying the
$\delta m$ term of $k_{\partial_{t}}^{tr}\equiv k^{tr}_{(t)}$ associated to a
Killing vector field $\xi=\partial_{t}$. That is, we decompose

$k^{tr}_{(t)}=\left[k^{tr}_{(t)}\right]_{\delta m}\delta
m+\left[k^{tr}_{(t)}\right]_{\delta a}\delta
a+\left[k^{tr}_{(t)}\right]_{\delta n}\delta n\,,$

and the first term is explicitly

$\left[k^{tr}_{(t)}\right]_{\delta m}=\\\
\frac{\sin\theta}{8\pi}\left(1+\frac{\left(r^{2}+a^{2}-2anC+n^{2}\right)(r^{2}-a^{2}\cos^{2}\theta)}{\left(r^{2}+a\cos\theta(a\cos\theta+2n)+n^{2}\right)^{2}}\right)\,.$
(18)

Note that it is a function of the coordinates. Here we can already identify
the first term that will contribute to half of the final charge $m$. We can
also check why, for this particular charge, the asymptotic analysis works as
in the limit $r\to\infty$ the second term in the parenthesis reduces to 1.

Before continuing with the calculation, let us define some functions that will
allow us to compactly write our next expressions

$\displaystyle h_{\pm}$ $\displaystyle=$ $\displaystyle(a\pm n)^{2}+r^{2}\,,$
(19) $\displaystyle f_{\delta a_{\pm}}$ $\displaystyle=$ $\displaystyle
n(-m(a\pm n)^{2}\pm 2n(a\pm n)r+mr^{2})\,,$ (20) $\displaystyle f_{\delta
n_{\pm}}$ $\displaystyle=$ $\displaystyle m(a\pm n)^{2}(2a\pm n)\mp n(a\pm
n)(3a\pm n)r$ $\displaystyle+\,m(2a\pm 3n)r^{2}\mp 3nr^{3}\,.$ (21)

For the Killing vector field $\xi=\partial_{t}$, the terms from the cones
contributing to the charges in (6) are

$\displaystyle\int_{T_{S}}\boldsymbol{k}_{(t)}$ $\displaystyle=$
$\displaystyle\int_{r_{1}}^{r_{2}}dr\int_{0}^{2\pi}d\phi\,k^{\theta
t}_{(t)}\bigg{|}_{\theta=\pi-\epsilon}=K_{(t)}\bigg{|}_{\begin{subarray}{c}r=r_{1},\\\
\theta=\pi-\epsilon\end{subarray}}-K_{(t)}\bigg{|}_{\begin{subarray}{c}r=r_{2},\\\
\theta=\pi-\epsilon\end{subarray}}\,,$ (22)
$\displaystyle\int_{T_{N}}\boldsymbol{k}_{(t)}$ $\displaystyle=$
$\displaystyle\int_{r_{1}}^{r_{2}}dr\int_{0}^{2\pi}d\phi\,k^{\theta
t}_{(t)}\bigg{|}_{\theta=\epsilon}=K_{(t)}\bigg{|}_{\begin{subarray}{c}r=r_{1},\\\
\theta=\epsilon\end{subarray}}-K_{(t)}\bigg{|}_{\begin{subarray}{c}r=r_{2},\\\
\theta=\epsilon\end{subarray}}\,,$ (23)

where $\epsilon$ is the small opening angle for the cones. We take
$\epsilon\to 0$ after integrating on each surface. Note that we have included
an explicit intermediate term, and note that for that term $\theta$ names a
component of the tensorial density $k^{\mu\nu}_{(t)}$ but it is also a
coordinate where the spacetime function $k^{\mu\nu}_{(t)}$ can be evaluated.
We concentrate in the region with $r_{2}\to r$, then the two ending terms of
the last equations are, explicitly, given by

$\displaystyle K_{(t)}\big{|}_{north}$ $\displaystyle=$
$\displaystyle\lim_{\epsilon\to 0}K_{(t)}\big{|}_{\begin{subarray}{c}r,\\\
\theta=\epsilon\end{subarray}}$ $\displaystyle=$
$\displaystyle\frac{C+1}{2h_{+}^{2}}\left(n(a+n)h_{+}\delta m+f_{\delta
a_{+}}\delta a+f_{\delta n_{+}}\delta n\right)\,,$ $\displaystyle
K_{(t)}\big{|}_{south}$ $\displaystyle=$ $\displaystyle\lim_{\epsilon\to
0}K_{(t)}\big{|}_{\begin{subarray}{c}r,\\\ \theta=\pi-\epsilon\end{subarray}}$
$\displaystyle=$ $\displaystyle\frac{C-1}{2h_{-}^{2}}\left(n(a-n)h_{-}\delta
m+f_{\delta a_{-}}\delta a+f_{\delta n_{-}}\delta n\right)\,.$

For the sphere contribution we have

$\displaystyle\oint_{S}\boldsymbol{k}_{(t)}$ $\displaystyle=$
$\displaystyle\lim_{\epsilon\to
0}\left[\int_{\epsilon}^{\pi-\epsilon}d\theta\int_{0}^{2\pi}d\phi\
k^{tr}_{(t)}\right]$ $\displaystyle=$
$\displaystyle\left(1-\frac{(C-1)n(a-n)}{2h_{-}}-\frac{(C+1)n(a+n)}{2h_{+}}\right)\delta
m$ $\displaystyle-\left(\frac{(C-1)f_{\delta
a_{-}}}{2h_{-}^{2}}+\frac{(C+1)f_{\delta a_{+}}}{2h_{+}^{2}}\right)\delta a$
$\displaystyle-\left(\frac{(C-1)f_{\delta
n_{-}}}{2h_{-}^{2}}+\frac{(C+1)f_{\delta n_{+}}}{2h_{+}^{2}}\right)\delta
n\,.$

Finally, the variation of the charge reads (see Eq. (7))

$\delta
M\equiv\not{\delta}Q_{(t)}=\oint_{S}\boldsymbol{k}_{(t)}+K_{(t)}\big{|}_{south}+K_{(t)}\big{|}_{north}=\delta
m\,,$ (27)

which happens to be trivially integrable. We see that most of the
contributions to the charge nicely cancel such that the final result does not
depend on coordinates (here the radius of the sphere). Note how complicated
the $r-$dependency of each term may be. Still, the whole charge is independent
as the method ensures conservation by virtue of
$\partial_{\mu}k^{\mu\nu}_{(t)}\approx 0$.
The computation of the (varied) angular momentum $\delta
J\equiv\not{\delta}Q_{-\partial_{\phi}}=\delta\left((a-3Cn)m\right)$
associated with the Killing vector field $\xi=-\partial_{\phi}$ can be carried
out following the same steps as the derivation of the mass.

## References

  * (1) A. H. Taub, Empty space-times admitting a three parameter group of motions, Annals of Mathematics 53 472-490 (1951).
  * (2) E. Newman, L. Tamburino, and T. Unti, Empty-space generalization of the Schwarzschild metric, Journal of Mathematical Physics 4 (1963), no. 7 915–923.
  * (3) C. W. Misner, “The Flatter regions of Newman, Unti and Tamburino’s generalized Schwarzschild space,” J. Math. Phys. 4, 924-938 (1963) doi:10.1063/1.1704019
  * (4) Bonnor, W. B. 1969 Proc. Camb. Phil. Soc. 66 145
  * (5) A. Chamblin, R. Emparan, C. V. Johnson and R. C. Myers, “Large N phases, gravitational instantons and the nuts and bolts of AdS holography,” Phys. Rev. D 59, 064010 (1999) doi:10.1103/PhysRevD.59.064010 [arXiv:hep-th/9808177 [hep-th]]. Copy to ClipboardDownload
  * (6) A. M. Ghezelbash, R. B. Mann and R. D. Sorkin, Can. J. Phys. 86, 579-582 (2008) doi:10.1139/p07-184
  * (7) R. Emparan, C. V. Johnson and R. C. Myers, “Surface terms as counterterms in the AdS / CFT correspondence,” Phys. Rev. D 60, 104001 (1999) doi:10.1103/PhysRevD.60.104001 [arXiv:hep-th/9903238 [hep-th]].
  * (8) R. B. Mann, “Misner string entropy,” Phys. Rev. D 60, 104047 (1999) doi:10.1103/PhysRevD.60.104047 [arXiv:hep-th/9903229 [hep-th]].
  * (9) R. B. Mann, “Entropy of rotating Misner string space-times,” Phys. Rev. D 61, 084013 (2000) doi:10.1103/PhysRevD.61.084013 [arXiv:hep-th/9904148 [hep-th]].
  * (10) D. Garfinkle and R. B. Mann, “Generalized entropy and Noether charge,” Class. Quant. Grav. 17, 3317-3324 (2000) doi:10.1088/0264-9381/17/16/314 [arXiv:gr-qc/0004056 [gr-qc]].
  * (11) C. V. Johnson, “Thermodynamic Volumes for AdS-Taub-NUT and AdS-Taub-Bolt,” Class. Quant. Grav. 31, no.23, 235003 (2014) doi:10.1088/0264-9381/31/23/235003 [arXiv:1405.5941 [hep-th]].
  * (12) C. V. Johnson, “The Extended Thermodynamic Phase Structure of Taub-NUT and Taub-Bolt,” Class. Quant. Grav. 31, 225005 (2014) doi:10.1088/0264-9381/31/22/225005 [arXiv:1406.4533 [hep-th]].
  * (13) R. Clarkson, L. Fatibene and R. B. Mann, “Thermodynamics of (d+1)-dimensional NUT charged AdS space-times,” Nucl. Phys. B 652, 348-382 (2003) doi:10.1016/S0550-3213(02)01143-4 [arXiv:hep-th/0210280 [hep-th]].
  * (14) D. Astefanesei, R. B. Mann and E. Radu, “Nut charged space-times and closed timelike curves on the boundary,” JHEP 01, 049 (2005) doi:10.1088/1126-6708/2005/01/049 [arXiv:hep-th/0407110 [hep-th]].
  * (15) G. Kalamakis, R. G. Leigh and A. C. Petkou, “Aspects of holography of Taub-NUT- AdS4 spacetimes,” Phys. Rev. D 103, no.12, 126012 (2021) doi:10.1103/PhysRevD.103.126012 [arXiv:2009.08022 [hep-th]].
  * (16) L. Ciambelli, C. Corral, J. Figueroa, G. Giribet and R. Olea, “Topological Terms and the Misner String Entropy,” Phys. Rev. D 103, no.2, 024052 (2021) doi:10.1103/PhysRevD.103.024052 [arXiv:2011.11044 [hep-th]].
  * (17) E. Arratia, C. Corral, J. Figueroa and L. Sanhueza, “Hairy Taub-NUT/bolt-AdS solutions in Horndeski theory,” Phys. Rev. D 103, no.6, 064068 (2021) doi:10.1103/PhysRevD.103.064068 [arXiv:2010.02460 [hep-th]].
  * (18) S. W. Hawking and C. J. Hunter, “Gravitational entropy and global structure,” Phys. Rev. D 59, 044025 (1999) doi:10.1103/PhysRevD.59.044025 [arXiv:hep-th/9808085 [hep-th]].
  * (19) S. W. Hawking, C. J. Hunter and D. N. Page, “Nut charge, anti-de Sitter space and entropy,” Phys. Rev. D 59, 044033 (1999) doi:10.1103/PhysRevD.59.044033 [arXiv:hep-th/9809035 [hep-th]].
  * (20) S. W. Hawking, “Gravitational Instantons,” Phys. Lett. A 60, 81 (1977) doi:10.1016/0375-9601(77)90386-3
  * (21) G. Clément, D. Gal’tsov and M. Guenouche, “Rehabilitating space-times with NUTs,” Phys. Lett. B 750, 591-594 (2015) doi:10.1016/j.physletb.2015.09.074 [arXiv:1508.07622 [hep-th]].
  * (22) G. Clément, D. Gal’tsov and M. Guenouche, “NUT wormholes,” Phys. Rev. D 93, no.2, 024048 (2016) doi:10.1103/PhysRevD.93.024048 [arXiv:1509.07854 [hep-th]].
  * (23) R. A. Hennigar, D. Kubizňák and R. B. Mann, “Thermodynamics of Lorentzian Taub-NUT spacetimes,” Phys. Rev. D 100, no.6, 064055 (2019) doi:10.1103/PhysRevD.100.064055 [arXiv:1903.08668 [hep-th]].
  * (24) P. A. Cano and D. Pereñiguez, [arXiv:2101.10652 [hep-th]].
  * (25) H. Gonzalez, D. Grumiller, W. Merbis and R. Wutte, EPJ Web Conf. 168, 01009 (2018) doi:10.1051/epjconf/201816801009 [arXiv:1709.09667 [hep-th]].
  * (26) A. Ballon Bordo, F. Gray, R. A. Hennigar and D. Kubizňák, “The First Law for Rotating NUTs,” Phys. Lett. B 798, 134972 (2019) doi:10.1016/j.physletb.2019.134972 [arXiv:1905.06350 [hep-th]].
  * (27) A. B. Bordo, F. Gray and D. Kubizňák, “Thermodynamics and Phase Transitions of NUTty Dyons,” JHEP 07, 119 (2019) doi:10.1007/JHEP07(2019)119 [arXiv:1904.00030 [hep-th]].
  * (28) R. Durka, “The first law of black hole thermodynamics for Taub-NUT spacetime,” [arXiv:1908.04238 [gr-qc]].
  * (29) A. B. Bordo, F. Gray, R. A. Hennigar and D. Kubizňák, “Misner Gravitational Charges and Variable String Strengths,” Class. Quant. Grav. 36, no.19, 194001 (2019) doi:10.1088/1361-6382/ab3d4d [arXiv:1905.03785 [hep-th]].
  * (30) A. Ballon Bordo, F. Gray, R. A. Hennigar and D. Kubizňák, “The First Law for Rotating NUTs,” Phys. Lett. B 798, 134972 (2019) doi:10.1016/j.physletb.2019.134972 [arXiv:1905.06350 [hep-th]].
  * (31) M. Zhang and J. Jiang, “NUT charges and black hole shadows,” Phys. Lett. B 816, 136213 (2021) doi:10.1016/j.physletb.2021.136213 [arXiv:2103.11416 [gr-qc]].
  * (32) J. Lee and R. M. Wald, J. Math. Phys. 31, 725-743 (1990) doi:10.1063/1.528801
  * (33) V. Iyer and R. M. Wald, Phys. Rev. D 50, 846-864 (1994) doi:10.1103/PhysRevD.50.846 [arXiv:gr-qc/9403028 [gr-qc]].
  * (34) R. M. Wald and A. Zoupas, Phys. Rev. D 61, 084027 (2000) doi:10.1103/PhysRevD.61.084027 [arXiv:gr-qc/9911095 [gr-qc]].
  * (35) G. Barnich and F. Brandt, Nucl. Phys. B 633, 3-82 (2002) doi:10.1016/S0550-3213(02)00251-1 [arXiv:hep-th/0111246 [hep-th]].
  * (36) E. Frodden and D. Hidalgo, “Surface Charges Toolkit for Gravity,” Int. J. Mod. Phys. D 29, no.06, 2050040 (2020) doi:10.1142/S0218271820500406 [arXiv:1911.07264 [hep-th]].
  * (37) Demianski M and Newman E T 1966 Bull. Acad. Polon. Sci. Ser. Math. Astron. Phys. 14 653\.
  * (38) V. S. Manko, J. Martin and E. Ruiz, Class. Quant. Grav. 23, 4473-4484 (2006) doi:10.1088/0264-9381/23/13/011 [arXiv:gr-qc/0603002 [gr-qc]].
  * (39) J. B. Griffiths and J. Podolsky, “Exact Space-Times in Einstein’s General Relativity,” doi:10.1017/CBO9780511635397.
  * (40) J. B. Griffiths and J. Podolsky, Int. J. Mod. Phys. D 15, 335-370 (2006) doi:10.1142/S0218271806007742 [arXiv:gr-qc/0511091 [gr-qc]].