NUTs and bolts beyond Lovelock

We construct a plethora of new Euclidean AdS-Taub-NUT and bolt solutions of several four- and six-dimensional higher-curvature theories of gravity with various base spaces $\mathcal{B}$. In $D=4$, we consider Einsteinian cubic gravity, for which we construct solutions with $\mathcal{B}=\mathbb{S}^2,\mathbb{T}^2$. These represent the first generalizations of the Einstein gravity Taub-NUT/bolt solutions for any higher-curvature theory in four dimensions. In $D=6$, we show that no new solutions are allowed for any Generalized quasi-topological gravity at cubic order. They exist however when we consider quartic Quasi-topological and Generalized quasi-topological terms, for which we construct new solutions with $\mathcal{B}=\mathbb{CP}^2,\mathbb{S}^2\times\mathbb{S}^2,\mathbb{S}^2\times\mathbb{T}^2,\mathbb{T}^2\times\mathbb{T}^2$. In all cases, the solutions are characterized by a single metric function, and they reduce to the corresponding ones in Einstein gravity when the higher-curvature couplings are set to zero. While the explicit profiles must be constructed numerically (except for a few cases), we obtain fully analytic expressions for the thermodynamic properties of all solutions. The new solutions present important differences with respect to Einstein gravity, including regular bolts for arbitrary values of the NUT charge, critical points, and re-entrant phase transitions.


JHEP10(2018)095
Higher curvature theories of gravity are proving to be increasingly useful in providing us with knowledge connected with fundamental questions in gravitational physics. Quadratic curvature theories are renormalizable [1], and it has long been realized that corrections to the Einstein-Hilbert action that go like powers in the curvature generically arise as low energy corrections from a UV complete theory of gravity, e.g. string theory [2]. Furthermore, in the context of the AdS/CFT correspondence conjecture, higher curvature corrections correspond to 1/N corrections in the large N limit of the dual CFT, allowing investigation of a much broader class of CFTs [3][4][5][6][7][8]. Imposing the condition that higher-curvature gravity be ghost-free (even if only on a constant-curvature background) severely limits the number of sensible theories available, since they generically contain such excitations. The most well-known class of higher curvature theories that are ghost free on maximally symmetric backgrounds is Lovelock gravity [9], in which the dimensionally extended Euler densities are included in the gravitational action. In D ≥ 2k + 1, the k th order Euler density is non-trivial; in this sense Lovelock gravity is a natural extension of Einstein gravity, and is the unique higher curvature theory maintaining second order field equations for the metric.
However other ghost-free higher curvature theories exist, and one class of particular interest has recently been identified [10][11][12][13][14][15][16][17]. The Lagrangian of this new family can be written schematically as L(g ab , R cdef ) = 1 16πG where L is some length scale, µ n are independent dimensionless couplings, and R (n) are certain linear combinations of order-n densities constructed from contractions of the metric and the Riemann tensor. The theories are characterized by the following properties: i) they have second-order equations of motion when linearized around any maximally symmetric spacetime, i.e., just like Einstein gravity, they only propagate a massless and traceless graviton on such backgrounds; ii) they possess a continuous and well-defined Einstein gravity limit corresponding to µ n → 0 for all µ n ; iii) they admit generalizations of the Schwarzschild-(A)dS black hole -so they reduce to it in the Einstein gravity limitcharacterized by a single function, where dΣ 2 (D−2) is the metric of the horizon cross sections (not necessarily spherical); iv) the function f (r) is determined from (at most) a second-order differential equation which, for a fixed set of µ n , admits a unique black-hole solution completely characterized by its ADM energy and which, at least in the spherically symmetric case, describes the exterior field of matter distributions with that symmetry [14]; v) the thermodynamic properties of such black holes can be obtained from a system of algebraic equations with no free parameters.

JHEP10(2018)095
The above class of theories can be subdivided if one considers more restrictive criteria. In particular, replacing i) by the requirement that the full non-linear equations of the theory are second order selects the Lovelock family [9,18,19]. Keeping i) as it is, but replacing instead iv) by the requirement that f (r) is determined by an algebraic equation, selects a more general class of theories, known as Quasi-topological gravities [20][21][22][23] (which of course include Lovelock as particular cases). More generally, the family of higher-curvature gravities satisfying i)-v) is larger than the Quasi-topological one. The missing theories possess black holes whose metric function is determined by second-order differential equations, and the full set has been coined "Generalized Quasi-topological Gravity" (GQTG) in [13].
One intriguing feature of these new theories, in contradistinction to those belonging to the Quasi-topological subset (except for Einstein gravity itself) is that some of them are nontrivial in D = 4. The simplest possible case of that kind, and the first to be identified, corresponds to a single additional cubic term and goes by the name of four-dimensional Einsteinian cubic gravity 1 (ECG) [10], whose action is given in (2.1) below. Many examples of GQTG theories in general dimensions have now been constructed, and their respective black hole solutions studied and characterized [11][12][13][14][15][16][17][24][25][26][27][28][29].
A different class of exact static solutions of Einstein gravity is given by the Taub-NUT family. The Euclidean section of the corresponding metrics can be written as which, in even dimensions, can be understood as U(1) fibrations over (D − 2)-dimensional Kähler-Einstein base spaces B with metric g B . In (1.3), τ is a periodic coordinate parametrizing the S 1 , and J = dA B is the Kähler form on B. The non-triviality of the fibration is controlled by the presence of a non-zero parameter n, customarily called "NUT charge". Depending on the dimension of the set of fixed points of the U(1) isometry -namely those for which V B (r) = 0 -the solution is said to be a "NUT" or a "bolt". Taub-bolt solutions are characterized by (D − 2)-dimensional fixed-point sets, whereas smaller dimensionalities give rise to Taub-NUT solutions. It has been known for some time that NUT-charged solutions exist in Lovelock gravity [30][31][32][33]. A broad understanding of their thermodynamics remains an ongoing subject of investigation [34][35][36][37][38][39][40][41][42][43][44][45][46] since it was realized that their contribution to the entropy does not obey the area law, even in Einstein gravity [47,48]. A recent review of Taub-NUT spacetimes and their symmetries has appeared [49].
On general grounds, one expects two independent functions to be required to describe Taub-NUT solutions in general higher-curvature gravities. The relevant observation for us is that both for Einstein gravity and Gauss-Bonnet, all Taub-NUT solutions are characterized by a single function for each choice of base space. This is analogous to the situation encountered for static black-hole solutions. One is then naturally led to wonder whether 1 The D-dimensional version of ECG was originally obtained as the most general cubic theory defined in a dimension-independent way -i.e., so that the relative coefficients of the cubic invariants involved do not depend on D -possessing second-order linearized equations on general maximally symmetric backgrounds [10]. However, it is only for D = 4 that ECG additionally satisfies properties ii)-v) [11,12].

JHEP10(2018)095
the rest of GQTG theories also admit generalizations of the Einstein gravity Taub-NUT solutions characterized by a single function, V B (r), just like they admit generalizations of the Schwarzschild black hole with that property. The answer turns out to be yes and, as we show here, a plethora of new Taub-NUT and Taub-bolt solutions of the form (1.3) can be constructed in various dimensions and for different base spaces.
Each choice of D, base space and Taub-NUT class has its peculiarities, but some aspects of our construction can be explained in general. First of all, and in a similar fashion to what occurs with black holes, inserting ansatz (1.3) in the equations of motion of the corresponding GQTG theory we will observe that, whenever the corresponding theory admits Taub-NUT solutions of that form, 2 they reduce to a single third-order equation for V B (r). Interestingly, this equation always admits a simple integrable factor which allows us to integrate it once. Hence, in each case we are left with a single equation of the form (1.4) where C is an integration constant related to the ADM energy of the solution. Also, in analogy with the black-hole case, one of the integration constants of this second-order differential equation will always be fixed by imposing the solutions to be locally asymptotically AdS. With regards to the second, recall that the defining properties characterizing V B (r) are, respectively, where β τ is the period of τ , and the bolt location satifies r b > n. While the first condition determines whether we are considering a NUT or a bolt, the second ensures that the solutions are smooth at r = n and r = r b , respectively and, when possible, it will fix the other integration constant in (1.4) for our solutions. For concreteness, we shall restrict ourselves to four-and six-dimensional theories. In D = 4, we will focus on the simplest possible modification to the Einstein-Hilbert action, namely, the ECG term. For this, we will construct solutions with base spaces B = S 2 and T 2 . These are, to the best of our knowledge, the first higher-curvature generalizations of the Einstein gravity Taub-NUT solutions in four dimensions. Turning to D = 6, we find that no non-trivial solutions of the form (1.3) can be constructed at cubic order in the GQTG family. They will exist, however, when quartic invariants are included, and we will restrict ourselves to that case. For those, we will construct solutions with B = CP 2 , S 2 ×S 2 , S 2 ×T 2 and T 2 × T 2 .
Although we will not be able to solve (1.4) analytically for V B (r) in general (except for the critical theories), the thermodynamic properties of the solutions will be accesible in a fully analytic fashion, again similar to what happened for the black hole solutions constructed in [11][12][13][14][15][16]. The relation between the ADM energy of the solutions, the NUT charge and r b (when present) will be accessible in each case from the asymptotic and near JHEP10(2018)095 r = n or r = r b expansions. On the other hand, in order to compute the free energy of the solutions, we will make use of the method introduced in [27]. According to this, given some higher-curvature gravity with Lagrangian density L(g ef , R abcd ) whose linearized equations on pure AdS match those of Einstein gravity (up to a normalization of Newton's constant), the Euclidean on-shell action of any asymptotically AdS solution can be computed using the formula 3 where Ω D−2 ≡ 2π (D−1)/2 /Γ((D − 1)/2) is the area of the unit sphere S D−2 ,L is the AdS radius, and a * is the charge appearing in the universal contribution to the entanglement entropy across a spherical entangling surface S D−3 in the dual CFT. This quantity is related, for any higher-curvature theory of gravity, to the on-shell Lagrangian of the theory on pure AdS through [27,[55][56][57][58] a * = − π (D−1)/2LD (D − 1)Γ D−1 2 L| AdS . (1.7) In [27], it was also argued that the same counterterms required to produce finite on-shell actions for Einstein gravity solutions can also be used for higher-curvature gravities of this class if we weight them by the same overall coefficient. With minor modifications in the D = 6 case -see discussion below eq. (3.38) -associated with the fact that the solutions are only locally asymptotically AdS, eq. (1.6) satisfactorily removes all divergent terms in the corresponding on-shell actions, and yields thermodynamic masses that agree with the ADM ones in all cases. The structure of our paper is simple. In sections 2 and 3 we construct Taub-NUT/bolt solutions of D = 4 Einsteinian cubic gravity and D = 6 Quartic Generalized quasitopological gravities, respectively. In each case, we compute the relevant thermodynamic quantities of the solutions, with special emphasis on the most standard cases B = S 2 and B = CP 2 , for which we study the corresponding phase spaces finding interesting new phenomena. Subsection 2.3 is somewhat different from the rest. It is devoted to the critical limit of Einsteinian cubic gravity, for which the solutions can be constructed analytically. We conclude in section 4. In appendix A, we repeat the D = 4 analysis in section 2 and construct solutions for Einsteinian cubic gravity plus an additional quartic density of the Generalized quasi-topological class. In appendix B, we present a detailed calculation of the on-shell action for the D = 4 NUT solution with B = S 2 which should be illustrative of the method utilized for the other cases. Some details regarding our numerical computations can be found in appendix C.

JHEP10(2018)095
2 Four dimensions: Einsteinian cubic gravity The first theory we will consider is four-dimensional Einsteinian cubic gravity with a negative cosmological constant [10].

Its Euclidean action reads
where the cubic density P is defined as The explicit form of the field equations of (2.1) can be found e.g., in [11]. The theory admits pure AdS 4 solutions of radiusL related to the action scale L byL Throughout the paper we will assume 0 ≤ µ ≤ 4/27, for which a unique branch of stable AdS vacua reducing to the Einstein gravity one as µ → 0 exists. In general, stable vacua exist for µ < 0 as well, but these are eliminated by the requirement that black holes have positive energy [11,12]. On the other hand, values of µ larger than 4/27 always give rise to unstable vacua. The "critical" limit of the theory [17], corresponding to µ = 4/27, warrants special attention. For that value of the coupling, the effective Newton constant diverges, and a number of simplifications take place, including the existence of analytic black hole solutions -as well as various exotic results from the point of view of a putative CFT dual [27].
Let us consider a metric ansatz with NUT charge n of the form (1.3) where, initially, we choose the base spaces B = S 2 , T 2 and H 2 , although we shall only construct explicit solutions for the first two. The base-space metrics and 1-forms appearing in (1.3) can then be written, respectively, as We stress again that the most general ansatz for a Taub-NUT metric in a general highercurvature gravity should involve an additional function -for example, g τ τ = V B (r)N B (r) 2 instead. It is a remarkable and highly nontrivial property of ECG that, when evaluated on (1.3) with the above choice of base spaces, its field equations reduce to a single differential JHEP10(2018)095 equation for V B (r). This is given by (omitting the 'B' subscript to reduce the clutter) where k = +1, 0, −1 for S 2 , T 2 and H 2 , respectively. Despite its challenging appearance, the above equation has two remarkable properties. First, it is of third order, instead of fourth, which is what one would have naively expected. Second, it allows for an integrable factor: after multiplying by (1 − n 2 /r 2 ), the equation becomes a total derivative and it can be integrated once. By doing so, we are left with a second-order differential equation of the form (1.4), namely where C is an integration constant which will be related to the energy of the solution. We now require the metric (1.3) to be locally asymptotically AdS; as a consequence we must demand V (r) → f ∞ r 2 L 2 + O(r) as r → +∞. Performing a 1/r expansion we find The effective Newton constant of the theory is given by so, at least in the spherical case, we can identify the integration constant in (2.6) with the ADM mass of the solution as C = GM . In an abuse of notation, we will use this definition for all base spaces. Now, note that since (2.6) is a second-order differential equation, it possesses a two-parameter family of solutions, of which (2.7) corresponds to a particular JHEP10(2018)095 one. In order to find the remaining asymptotic solutions, let us write V (r) = V p (r)+ r 2 L 2 g(r) and expand linearly in g. Taking into account only the leading terms when r → +∞, we find that g satisfies the following equation Leaving aside the limiting values µ = 0, 4/27, the general solution is given by where AiryAi[x] and AiryBi[x] are the Airy functions of the first and second kind, respectively. When GM µ > 0, the solution involving AiryBi grows exponentially, while the one with AiryAi decays. Therefore, we must set B = 0 in order for the solutions to be locally asymptotically AdS. Hence, we learn that the asymptotic boundary condition is fixing one of the integration constants in (2.6). The remaining one will be fixed by the corresponding regularity conditions in the bulk, as we will show in the following sections. When GM µ < 0, the solutions (2.10) have an oscillatory character and they are all singular at infinity (except the trivial one, g = 0). To remove this behaviour we would need to set both A and B to zero, which would fully specify the solution. This would leave us with no integration constants to impose regularity in the bulk. This behaviour is very similar to that found for the static black hole solutions of the theory [11,12] -see also [16], and leads us to choose µ ≥ 0, so that the solutions with GM µ > 0 have positive energy.

Einstein gravity
In the following subsections we will consider the base spaces S 2 and T 2 independently, and we will construct new Taub-NUT and bolt solutions for them for general values of µ. It is illustrative however to start analyzing the Einstein gravity case, for which the analysis can be performed at the same time for all base spaces. Indeed, if we set µ = 0, (2.6) can be easily solved for V B (r). Imposing the NUT condition V B (r = n) = 0 first, one is left with where we already fixed the integration constant as The regularity condition (1.5) imposes which means that τ cannot be a compact coordinate for B = T 2 or, in other words, the solution is extremal, in the sense that the temperature T ≡ 1/β τ is forced to vanish. Similarly, for B = H 2 , one finds that the period of τ would need to be negative. This means JHEP10(2018)095 that V H 2 (r) actually becomes negative for values of r greater than n, which is forbidden by assumption. Hence, no regular Taub-NUT solution exists in that case for Einstein gravity. If we impose the bolt condition V B (r = r b ) = 0 instead, we find where the integration constant was fixed as The regularity condition (1.5) fixes now the bolt radius as a function of n and β τ , namely (2. 16) In order for each solution to be allowed, it must be such that r b > n. Furthermore, the quantity inside the square root must be positive, which restricts the allowed values of n for which the corresponding solutions exists.
On general grounds, in order to remove the so-called Misner string [59], an additional condition must be imposed on β τ both for NUT and bolt solutions when B = S 2 . As we explain in the next subsection, this reads β τ = 8πn. It is a remarkable (and peculiar) fact that in Einstein gravity eq. (1.5) automatically implements this condition in the case of the NUT solution. In general, both conditions must be imposed separately.

B = S 2
Let us now turn on the Einsteinian cubic gravity coupling. We begin by assuming the base space to be the one-dimensional complex projective space CP 1 or, equivalently, the two-dimensional round sphere, B = S 2 . Then, the metric (1.3) reads This metric has "wire singularities" at θ = 0, π, for which it becomes noninvertible. As shown by Misner [59], it can nevertheless be made regular everywhere using two coordinate patches. The idea is to define new coordinates τ ± = τ ± 2nφ covering the θ ≥ π/2 and θ ≤ π/2 regions respectively. In the overlap region, τ + = τ − + 4nφ, and since β φ = 2π, one is forced to impose the periods of τ ± to be β τ ± = 8πn. For clarity reasons, in what follows we will work with the metric (2.17) in a single patch, but taking into account that the period of τ is related to the NUT charge through β τ = 8πn. Observe that this condition is a consequence of choosing B = S 2 and does not depend on the theory. When combined with the general regularity condition (1.5), this gives rise to the conditions V S 2 (r = n) = 1/(2n) and V S 2 (r = r b ) = 1/(2n) respectively for NUTs and bolts. The function V S 2 in (2.17) is determined from (2.6) with k = 1. Using the asymptotic expansion (2.7), we see that when r → +∞ the metric induced on a constant-r hypersurface JHEP10(2018)095 is given by where we have introduced the angle coordinate ψ = τ /(2n), whose period is 4π. When 4f ∞ n 2 = L 2 , the previous metric is the one of a round S 3 . For any other value of n, it is the metric of a squashed sphere, and it is customary [60][61][62][63][64] to rewrite the NUT charge in terms of a 'squashing parameter' α as 4f ∞ n 2 /L 2 = 1/(1 + α). In order to specify the solution, we need to choose a boundary condition at some finite r = r b . Depending on whether we choose r b = n, or r b > n, we will be considering Taub-NUT or Taub-bolt solutions.

Taub-NUT solutions
As we have explained, the Euclidean Taub-NUT metric is characterized by the conditions V S 2 (r = n) = 0 and V S 2 (r = n) = 1/(2n). Let us then expand V S 2 (r) around r = n as for some a i . Plugging this expansion into (2.6), we observe that the O (r−n) and O (r−n) 2 equations are automatically satisfied, whereas the O(1) one gives rise to the following relation between the mass and the NUT charge, Observe that this reduces to the Einstein gravity expression (2.12) for µ = 0. The following term in the expansion gives a relation between a 3 and a 2 , which we can use to write the former as a function of the latter, a 3 (a 2 ). Similarly, the following term allows us to obtain a 4 (a 2 ), and so on. Hence, as in the black hole case [11,12], the full series is determined by a single free parameter a 2 . This parameter must be chosen in a way such that B = 0 in (2.10), which ensures that the solution is locally asymptotically AdS. In practice, the shooting method can be used to identify a 2 for each value of µ, so that the near r = n expansion yields a good approximation to the exact solution that connects with the asymptotic expansion (2.7). There is a unique a 2 for each µ that does the job, corresponding to a unique Taub-NUT solution in each case. We plot the metric function V S 2 (r) for different values of µ in figure 1. These solutions generalize the Einstein gravity Taub-NUT solution (the red curve in figure 1), whose metric function is given by (2.11) with k = 1. As we can see, the qualitative behaviour of V S 2 (r) is very similar to that of Einstein gravity for nonvanishing values of the ECG coupling. One peculiarity of (2.20) for nonvanishing µ is that the mass becomes negative for small values of n. In particular, the mass is non-negative only when (2.21) The existence of a finite lower bound for n is a new feature, which does not occur for Einstein gravity. Indeed, in that case (2.21) becomes L 2 /(4f ∞ ) ≥ n 2 ≥ 0. For general values of the gravitational coupling, we cannot expect the solution to exist whenever n lies outside the interval in (2.21), because the asymptotic behaviour for negative masses is pathological. Indeed, as we explained in the discussion below eq. (2.10), negative mass solutions would be highly oscillating at infinity, and hence they are not asymptotically AdS. Solutions with zero mass occur when either the upper or the lower bounds are saturated. In terms of µ, the M = 0 condition reads In the case of Einstein gravity, the possibilities are n 2 = L 2 /4 and n = 0, for which the solution reduces to pure Euclidean AdS 4 foliated by round S 3 slices and S 1 × S 2 slices, respectively 4 [65]. For any nonvanishing value of µ, (2.22) is satisfied identically for n 2 = L 2 /(4f ∞ ), which can be straightforwardly checked using (2.3). In this case, the solution also reduces to pure AdS 4 , the metric factor being simply given by V S 2 (r) = f ∞ r 2 /L 2 −1/4. Besides this solution, there exists another one obtained by choosing n to saturate the lower bound in (2.21), and which is analogous to the n = 0 one in Einstein gravity. Interestingly, this solution no longer reduces to pure Euclidean AdS 4 for 0 < µ < 4/27 but, rather, it has a nontrivial profile. In the critical limit, µ = 4/27, the range allowed by (2.21) collapses to a single possible value, corresponding to n 2 = L 2 /6. In that case, the solution does correspond to pure AdS 4 . For other values of n, the critical solution can also be accessed

JHEP10(2018)095
analytically (see section 2.3), and the result for the metric function reads This solution has a vanishing mass parameter M = 0, namely, it only exists if we fix the integration constant C to zero in (2.6). Note that this solution has V cr S 2 (n) = 3n/L 2 ; it has a conical singularity at r = n in all cases but one, corresponding to the value n 2 = L 2 /6, for which it becomes pure AdS 4 , as mentioned above.
Let us now evaluate the on-shell action of the solutions. In order to do so, we make use of the generalized action (1.6) where, for ECG, the charge a * is given by a * = (1 + 3µf 2 ∞ )L 2 /(4G). Using this, the full ECG action takes the form (2.24) where R stands for the Ricci scalar of the induced metric on the boundary. The last two terms in the second line are the standard counterterms in D = 4 which, as explained in the Introduction, also appear weighted by a * without further modification according to the prescription in [27]. A detailed evaluation of all the terms in the above expression for our new Taub-NUT solutions, which can be performed fully analytically using the near r = n and asymptotic expansions, is presented in appendix B. It yields the following remarkably simple finite answer This reduces to the free energy of the corresponding Einstein gravity Taub-NUT solution when µ = 0, as it should. The energy and entropy can be easily obtained now from E = ∂I E /∂β and S = βE − I E . Using this, we find that the energy precisely matches the result for the ADM mass obtained in (2.20), E = M , which is a highly nontrivial check of the calculation, whereas for the entropy we obtain which is not given by a simple area law due to contributions from the Misner string. It is also possible to consider the thermodynamics of these NUT charged solutions from the perspective of extended phase space thermodynamics. Within this framework, one introduces potentials conjugate to the cosmological constant -interpreted as a pressure P = −Λ/(8πG) -and any higher-curvature couplings that appear in the action [66,67]. These considerations are motivated by scaling arguments, since without these terms the Smarr relation fails to hold. In the case of Taub-NUT solutions in ECG, the extended first law reads where we have restored the dimensions to the ECG coupling constant. The new potentials read (2.28)

JHEP10(2018)095
Interestingly, the thermodynamic volume here is precisely the same as for Taub-NUT solutions in Einstein gravity [41]. The same conclusion holds for the thermodynamic volume of black holes in higher-curvature gravities that belong to the generalized quasi-topological class. With the thermodynamic quantities defined as above, the Smarr relation that follows directly from a scaling argument is found to hold: In figure 2 we plot the Euclidean action (2.25) for several values of µ. Dashed lines correspond to negative values of the mass, whereas solid lines correspond to solutions with M > 0. As we mentioned earlier, in principle we only expect solutions with positive mass M > 0 to exist. A numerical analysis seems to confirm this, since we were not able to construct any solution with M < 0. This also constrains the validity of the thermodynamic expressions (2.25) and (2.26) to the interval defined by (2.21). This interval becomes smaller as µ grows, and it reduces to a single point, n 2 = L 2 /6, in the critical limit. Interestingly, we observe that the free energy of the critical theory solutions (solid gray curve) acts as an envelope of all possible solutions with positive mass and arbitrary values of µ. Observe that this free energy cannot be obtained from (2.25) in the µ → 4/27 limit; the same applies to the mass, which cannot be obtained from (2.20). The correct result for the on-shell action associated to the critical solutions with metric function (2.23) reads however As we mentioned above, all these solutions except for the one with n 2 = L 2 /6 have conical singularities, so the result must be taken with care -e.g., the mass cannot derived from (2.30) using standard thermodynamic identities. It is a remarkable and somewhat striking fact that the free energy of this singular solution, as given by (2.30), precisely separates the free energies of negative-mass solutions from those corresponding to completely regular positive-mass solutions for general values of µ. The different nature of the critical solutions can also be seen from the fact that whenever µ = 4/27, the solutions with M = 0 correspond to values of n for which I E (n) is locally extremized, whereas the whole µ = 4/27 curve has M = 0.

Taub-bolt solutions
Let us now turn to bolt solutions. These are obtained by imposing V S 2 (r) to vanish for some r b > n, i.e., V S 2 (r b ) = 0, plus the regularity condition V S 2 (r b ) = 1/(2n). If we plug a Taylor expansion for V S 2 (r) around r = r b including these conditions in (2.6), the equations corresponding to the first nontrivial orders give rise to two equations involving the mass of the solution M , the bolt radius r b , and the NUT charge n. These read . The relation r b (n) has several remarkable differences with respect to the Einstein gravity case. Indeed, for µ = 0, (2.32) has two nontrivial roots, given by (2.16), namely (2.33) Since we want r b to be real and larger than n, this implies that n/L < (2 − If µ is greater than this quantity, there is a two-to-one relation between n and r b for all n. All this is shown in figure 3. For the set of parameters for which a unique bolt solution exists, the profile of V S 2 (r) can be accessed numerically following exactly the same logic as for the NUT solutions. We plot the resulting metric functions for some values of µ in figure 4. We can also compute the on-shell action for the bolt solutions analogously to the NUT case. The final result can be written as .

JHEP10(2018)095
Using the chain rule and the relation (2.32), one can show again that E ≡ ∂ β I E = M as given in (2.31), which is a consistency check of the calculation. In addition, the entropy, given by S = βE − I E , reads .
Just as in the case of the NUT solutions, we can also study the thermodynamics of the bolts in extended phase space. The extended first law has the same form as in the NUT case, but now the thermodynamic volume and coupling potential read and satisfy the Smarr formula that follows from scaling, which is of the same form as in the NUT case -see (2.29). Note that, once again, the basic formula for the thermodynamic volume of the bolts is unaltered by the higher-curvature terms [41]. However the thermodynamic volume is implicitly sensitive to the ECG coupling since the ECG term is important in determining the value of r b for a given n.
Although we cannot solve (2.32) exactly, we can study its behaviour in several limits. For example, let us consider the new branch of solutions for which r b is close to n in the limit µ 1. We can expand r b in powers of µ 1/2 . To second order, we get For Einstein gravity, we get r b = n, and the solution reduces to the NUT one. However, for any given nonvanishing µ, we have two inequivalent solutions: the NUT constructed in the previous subsection, and this one. In particular, as opposed to the Einstein gravity case, a bolt solution does exist for n 2 = L 2 /(4f ∞ ), which corresponds to a nonsquashed spherical boundary geometry. Expansions for the free energy and the mass of this branch of solutions can be easily obtained in the µ 1 limit, the results being Note that the mass is nonvanishing when the boundary geometry is that of a round S 3 , namely, GM (n 2 = L 2 /(4f ∞ )) = 3µL/4 + O(µ 3/2 ), so the free energy is not extremized in that case. Instead, the maximum is reached for , which is also the M = 0 value. Greater values of n would give rise to negative mass solutions, as illustrated in figure 5. Note also that I E (n 2 = L 2 /(4f ∞ )) is greater than the one for the NUT solution, so that in the region near n 2 = L 2 /(4f ∞ ), the NUT would dominate the corresponding holographic partition function. We can also study the behaviour near n = 0, for which there is also a single bolt solution for each nonvanishing µ. We get approximately and for the free energy, If we set µ = 0 in these expressions, we recover the small n expansions for Einstein gravity bolts corresponding to the (+) root in (2.33). Observe that in the critical limit, µ = 4/27, the leading term disappears, and the on-shell action is finite for n = 0. We summarize the different possibilities in figure 5, where we plot I E for ECG bolt solutions for several values of µ. As we can see, the result is very different from that of Einstein gravity. There are two cases that we can distinguish: if 0 < µ < 0.001126, the diagram contains three branches, since there are three different bolt solutions; for µ > 0.001126 there is a single (elephant-shaped) branch. At µ 0.001126, we expect to have a critical point which would represent a second-order phase transition if the bolt solution were dominant. In all cases, the solutions exist for much larger values of n than in Einstein gravity. However, there is an additional upper bound on n coming from imposing M > 0.

Free energy comparison
Finally, let us compare the Euclidean action of NUT and bolt solutions in order to determine which one dominates the partition function. In figure 6 we compare the Euclidean actions for several values of µ. The case for Einstein gravity is shown in the top-left panel, and we can see that the NUT solution dominates in all the region n > L/6 7 − 2 √ 10, where a first order NUT/bolt phase transition takes place. In particular, there are no bolt solutions near the undeformed 3-sphere n 2 = L 2 /4. When we switch µ on, there are some drastic changes. Specially, we recall that for positive values of µ there are no solutions with negative mass. We plot with dashed lines the would-be Euclidean action of these solutions, but they do not actually exist. This has the effect of inducing zeroth-order phase transitions in the points

JHEP10(2018)095
where some solution ceases to exist. Another new feature is the existence of bolt solutions near the round 3-sphere n 2 = L 2 /(4f ∞ ) ≡ n 2 0 for all values of µ > 0. In all the cases, we observe that for n = n 0 the NUT solution (corresponding to pure AdS) dominates, but for n > n 0 the NUT solution does not exist because it would have negative mass. However, for values of n slightly larger than n 0 , there is still a bolt solution of positive mass, and a zeroth-order phase transition from NUT to bolt must take place at n 0 . For larger values of n, the bolt solution also acquires a negative mass and there are no solutions. The behaviour is more interesting in the region n < n 0 . In all the cases the NUT solution dominates until certain value n = n min , where there is a transition to a bolt solution. When µ < 0.00569, the transition is of first-order, as shown in top-right and bottom-left panels in figure 6. When µ > 0.00569, the mass of the NUT solution vanishes before the value of the Euclidean action crosses that of the bolt solution, and a zeroth-order phase transition takes place, as shown in the bottom-right panel. After that phase transition the bolt solution exists and dominates for 0 < n < n min .
The appearance of zeroth-order transitions in Taub-NUT solutions is a new feature whose interpretation is not clear to us. This seems to be a problem that only appears in four dimensions, since, as we will see, in six dimensions there is no restriction on the mass of the solutions.

B = T 2
Let us now consider a toroidal base space, so that the metric ansatz (1.3) reads Here, the coordinates (η, ζ) parametrize a T 2 with periods which we choose to be equal, β η,ζ = l. We note that, unlike the spherical case, the periodicity of the variable τ , which we denote β τ ≡ 1/T , is not a priori fixed in terms of n [30,68]. The function V T 2 (r) satisfies (2.6) with k = 0. From the general asymptotic expansion (2.7), we can obtain the metric of constant-r hypersurfaces for r n. This reads If we define z = τ /(2n), η/L = −x, ζ/L = y, this can be rewritten as Remarkably, when n 2 = L 2 /(4f ∞ ) -i.e., for the same value n for which in the B = S 2 case the corresponding boundary metric becomes that of a round S 3 -this reduces to the so-called Nil geometry 5 [69]. The appearance of this kind of geometry should not come as a surprise, as T m -bundles over tori T n are always compact 2-step nilmanifolds (and vice versa) [70] -in our case above, m = 1 and n = 2. 6 On the other hand, it is also natural to defineτ = √ f ∞ τ , whose periodicity is βτ = √ f ∞ β τ ≡ 1/T . Then, (2.43) reduces to the standard metric on T 3 for n = 0, so we can also understand (2.43) as a sort of twisted three-torus metric.

Taub-NUT solutions
Let us start with the NUT solutions. Just like in the previous section, we assume that V T 2 (r = n) = 0, and we impose V T 2 (r = n) = 4πT in order to avoid a conical singularity

JHEP10(2018)095
at the NUT. Then, we can consider a Taylor expansion around r = n of the form Plugging it into (2.6) and solving order by order in (r −n), we obtain the following relations for the first terms (2.49) The first equation fixes the "mass" M in terms of n, L, µ and T , while the rest give us relations between the coefficients of the expansion and the temperature. We can try to solve these relations in two inequivalent ways. The first possibility, which is the only one available for Einstein gravity, consists in setting T = 0 -see equations (2.11), (2.12) and (2.13) with k = 0. This solves the last two equations and, in fact, completely determines the series expansion for any value of µ, i.e., we can obtain a 2 , a 3 , etc., from the subsequent equations. The series is convergent in a vicinity of r = n. However, note that in that case, GM = −4n 3 /L 2 < 0. Hence, according to the general discussion about the asymptotic behaviour, we expect this solution to be pathological at infinity unless some miraculous fine-tuning occurs. Unfortunately, this is not the case, and when we solve (2.6) starting from the near-horizon expansion with T = 0, the oscillatory character appears asymptotically. We are then led to conclude that regular extremal NUT solutions do not exist for any allowed value of µ.
The second possibility is setting a 2 = 2πT /n, which solves the second equation. The following equations can be used to determine the remaining coefficients, which turn out to have a nonperturbative dependence on µ, e.g., Observe that these do not possess a finite limit when µ → 0. As a matter of fact, we have failed to construct these solutions for any nonvanishing value of µ different from the critical limit value µ = 4/27 -see section 2.3 -so we strongly suspect that no regular NUT solution exists for B = T 2 for any 0 < µ < 4/27.

Taub-bolt solutions
Fortunately, the situation is different for bolt solutions. In that case, we impose the existence of some r b > n such that near r = r b , This fixes V T 2 (r b ) = 0 and V T 2 (r b ) = 4πT . Again, we plug this expansion in (2.6), and from the first two terms we get . (2.53) As usual, the first equation fixes M , while the second relates r b to n and T . It turns out that for µ ≥ 0, there is a unique solution for T for every n and r b > n. The solution can be written explicitly as For a given n, this is a one-to-one relation between every r b > n and T > 0. In particular, we have lim r b →n T = 0. However, in order to keep GM positive, r b is bounded from below, r b ≥ nγ(µ), for some constant γ(µ). In particular, for the limiting cases µ = 0 and µ = 4/27, we have, respectively, γ(0) = 3 + 2 √ 3 2.5425, and γ(4/27) 4.0171. In each case, for a given n, the radius r b and the "mass" M are fixed by the periodicity of the coordinate τ . The remaining coefficients in the expansion (2.51) are fully determined once we choose a 2 , which is the only free parameter. Once again, this is fixed by demanding the solution to have the correct asymptotic behaviour. In all cases we find that there is one and only one value of a 2 for which this happens, and so the solutions are completely determined by n and T . In figure 7 we show some of the metric functions corresponding to these solutions computed numerically.

JHEP10(2018)095
Let us now study the thermodynamic properties of the solutions. The Euclidean action can be once again evaluated using the generalized action (1.6), following the same steps as in the previous subsections. The result is where we used (2.53) to simplify the result. This on-shell action should be understood as a function of T and n, which appear implicitly through r b . In the case of Einstein gravity, for which the metric function can be obtained analytically -see (2.14) with k = 0 -the result for the on-shell action can be written explicitly as a function of T and n. Using (2.16) with k = 0, one finds (2.56) Now we must account for the fact that we have an extended thermodynamic phase space, since n is in this case a free variable. However, n cannot be the appropriate thermodynamic variable as it has units of length instead of energy. Hence, let us define θ ≡ 1/n, which has the right units. We are going to see that this is indeed the canonical thermodynamic variable, as the energy we get matches the ADM one. Then, associated with T and θ, we have two potentials: the usual entropy S, and a new potential Ψ. In terms of the free energy F ≡ T I E , these are given by Finally, the energy is defined as E = F + T S + θΨ, so that, by construction it satisfies the first law The energy is given by Now, this is a thermodynamic energy, but the energy of the solution should be computed using the ADM formula, which in this case tells us that E ADM = M l 2 /(4πL 2 ). Using the expression for M given in (2.52), we have checked that both energies actually coincide E ADM = E. Hence, the introduction of the variable θ ≡ 1/n is crucial for the first law of black hole mechanics to hold in this case. This picture goes through nicely when the ECG term is turned on. In that case, it is convenient to express the thermodynamic quantities in terms of the rescaled temperature In terms of this, we have the free energy F (T , θ; µ) = √ f ∞T I E , which can be obtained from (2.55), and the thermodynamic potentials S(T , θ; µ) and Ψ(T , θ; µ) defined as in (2.57) (but with respect toT instead of T ). We find that where the ADM energy is now given by It is interesting to study the isotherms on the diagram of Ψ-θ. These are shown in figure 8 for µ = 0 and µ = 0.12. In the case of Einstein gravity -the same happens for small values of µ -the isotherms are monotonous. However, when µ is large enough, the diagram changes drastically. In that case, the isotherms develop a maximum, and the limit Ψ → 0 corresponding to θ → +∞ is approached from above instead of from below. However, the phase space seems to be free of critical points. We have seen that to satisfy the quantum-statistical relation and first law, it is necessary to treat θ = 1/n as a thermodynamic variable. If we also wish to satisfy the Smarr formula, then once again we must consider both Λ and µ to be thermodynamic parameters. The basic construction is identical to the B = S 2 case, but now we include θ as well. A simple computation then reveals the following for the thermodynamic volume and coupling potential where in the second equation above we note that it is the un-normalized temperature that appears (i.e. T rather thanT ). With these definitions, the extended first law and Smarr formula hold, with the latter being identical to eq. (2.29) with the additional term 2θΨ added.

JHEP10(2018)095
Let us close this section by mentioning the possibility that the solutions considered in this subsection can be relevant holographically. In that context, and in analogy to the Taubbolt solutions with B = S 2 , we expect them to represent saddle points in the semiclassical partition function for boundary theories living on deformed tori with metric (2.43). While in the spherical case the boundary is only characterized by n -or, equivalently, the squashing parameter α -, the B = T 2 case is richer, given that n andT are independent parameters in that case.

Exact Taub-NUT solutions in the critical limit
In this subsection we study the Taub-NUT solutions of critical ECG, which can be constructed analytically. As we mentioned earlier, when µ = 4/27, the only AdS vacuum has a length scaleL 2 = 2L 2 /3, and the linearized equations on that background vanish identically [17]. The field equations simplify considerably, which has allowed for the construction of analytic black hole solutions [17]. 7 In the case of NUT-charged metrics, a similar simplification takes place, and we find the following family of exact Taub-NUT solutions, . (2.65) As we mentioned before, in the case of a spherical base space, B = S 2 , this solution has a conical singularity at r = n, except for n 2 = L 2 /6, in whose case the solution is simply globally Euclidean AdS 4 -also known as H 4 . Hence, only the cases B = T 2 , H 2 are of interest in Euclidean signature. The solutions (2.65) can be analytically continued to Lorentzian signature in different ways, giving rise to very interesting metrics. For example if we make the replacement n → in and τ = it we get the following metric .
This metric is regular everywhere and, in fact, we can allow r to take values in the whole real line. Hence, this solution usually represents a wormhole or wormbrane, depending on the topology, connecting two asymptotically AdS 4 regions. The cases k = n = 0 and k = −6n 2 /L 2 = −1 are special as they correspond to pure AdS 4 . Let us introduce a new radial coordinate r = n cosh ρ/( 2/3L) , so that the metric reads which has an explicit wormhole character. In the spherical case B = S 2 the solution reads 8 The existence of special points in the parameter space in which the equations simplify has been previously reported in [72] in the case of Lovelock gravity, where Lorentzian AdS wormholes similar to those in [17] were constructed. 8 Let us note that a very similar NUT-charged wormhole solution was reported in [73] in the context of Einstein-Skyrme theory.

JHEP10(2018)095
This solution has the problem that it suffers from closed time-like curves, because the time coordinate must be periodic t → t + 8πn. The B = T 2 , H 2 cases are free of them, because there is no periodicity condition on the time coordinate. In particular, after some rescalings we can write the T 2 solution as where the NUT charge has been absorbed in the period of the coordinates x, y. However, we can also allow x and y to be noncompact. Interestingly, there is an inequivalent Lorentzian solution that can be obtained by rotating the coordinates as (t, y) → (iy, it). This reads Going back to the general solution (2.67), we can consider the following transformation: t → z, ρ → it, L → iL. Here we are changing the sign of L 2 , which amounts to changing the sign of the cosmological constant in the ECG action (2.1). Hence, the corresponding metric is a solution of the critical theory with a positive cosmological constant. The general solution reads These represent bouncing cosmologies with different topologies for the spatial sections connecting two asymptotically (NUT charged) de Sitter spaces for t → ±∞. The only exception is the case k = 1, n 2 = L 2 /6, which is actually de Sitter space foliated by S 3 spheres. Particularly relevant for cosmology is the flat case k = 0, which after rescaling of the coordinates can be written as The transverse geometry is again a Nil space. Interestingly, this solution represents a homogeneous but nonisotropic bouncing cosmology. Homogeneity follows from the fact that Nil space is a coset space and it possesses the isometries (x, y, z) → (x + a, y + b, z + c − a √ 6y/L), for arbitrary (a, b, c). Let us also mention in passing that this solution seems to be disconnected from the isotropic and homogeneous bouncing solution found in [17], since we do not recover it in any limit.

Six dimensions: quartic generalized quasi-topological gravities
We now move on to consider theories in six dimensions. In this case the generalized quasitopological gravity class [13][14][15] includes additional densities beyond those in D = 4. In JHEP10(2018)095 particular, the Gauss-Bonnet term X 4 is no longer topological. Analytic generalizations of the Einstein gravity static black-hole [74,75] and Taub-NUT/bolt solutions [30,31] have been constructed in the presence of this contribution. In particular, the Taub-NUT solutions of Gauss-Bonnet [30,31] and 3rd-order Lovelock gravity [33] were, prior to this paper, and to the best of our knowledge, the only known examples of solutions of that class for any higher-curvature gravity theory.
In principle, the six-dimensional GQTG family includes two nontrivial terms at cubic order, corresponding to the usual quasi-topological gravity density, plus an additional one. As observed in [20,21,76], including the quasi-topological gravity density in D ≥ 6 is equivalent, from the point of view of static black-hole solutions, to including the cubic Lovelock interaction. In D = 6, this is a topological term, and therefore no new nontrivial black holes exist in that case for quasi-topological gravity. They do exist, however, when the additional GQTG term is included [13]. Hence, following the same reasoning as for ECG in D = 4, one would have expected that new Taub-NUT solutions of the form (1.3) should also exist for GQTG. Remarkably, we find that this is not the case, and that no cubic theory admits nontrivial generalizations of the Einstein gravity or Gauss-Bonnet Taub-NUT solutions characterized by a single function in six dimensions (independent of the base manifold considered). While there do exist cubic theories that satisfy the necessary constraints to admit solutions of the form (1.3), it turns out that for any such theory the field equations on these spaces vanish identically. This is analogous to the C (i) D terms discussed in [15] and quasi-topological gravities of order n in D = 2n [76] for the case of static, spherically symmetric metrics. Mathematically, this arises due to the proportionality of certain cubic invariants on the class of metrics considered here.
Happily, at quartic order in curvature there exist non-trivial options of both quasitopological and generalized quasi-topological type. 9 We will not present a detailed classification of such theories here -see [22] and [15] -but limit ourselves to some brief remarks. Beginning from a general action containing all 26 possible quartic invariants [77], we constrain the action by imposing the conditions listed in appendix A of [15]. This selects theories admitting black hole solutions of the form (1.2) and which, as a consequence [14], do not propagate ghosts on maximally symmetric backgrounds. Next, for each of the possible four dimensional base manifolds listed below, we generate additional constraints to ensure the corresponding theory admits solutions of the form (1.3). Surprisingly, demanding the constraints to be simultaneously satisfied for all four base manifolds B = CP 2 , S 2 × S 2 , S 2 × T 2 and T 2 × T 2 results in a family of theories that yield trivial field equations. When one relaxes this condition, considering only a subset of the base manifolds, then nontrivial options exist.
The non-trivial theories can be classified into two groups: quasi-topological and generalized quasi-topological. For the base manifold S 2 × S 2 , the only nontrivial theories are of 9 Recall that the distinction between both classes comes from considering the theories for static spherically symmetric spacetimes. While both admit solutions of the form (1.2), the field equations for the quasi-topological theories reduce to algebraic polynomial equations for f (r), while for those of the generalized quasi-topological type, the metric function satisfies a non-linear second order differential equation in each case.

JHEP10(2018)095
the generalized quasi-topological type. The constraints can be satisfied simultaneously for B = CP 2 , S 2 × S 2 and T 2 × T 2 , resulting in a three parameter family of non-trivial theories, each making the same contribution to the field equations for a given base manifold. We can deduce all of the relevant physics by considering only one member of this class, which we denote as S below.
Excluding the base manifold S 2 ×S 2 , then quasi-topological options exist for all remaining base manifolds. The constraints can be satisfied simultaneously, resulting in a three parameter family of non-trivial quasi-topological theories. Again, each theory makes the same contribution to the field equations for a given base manifold, and we need only consider a single member of this family, which we denote as Z below. Thus, in six dimensions, we will consider the following two quartic Lagrangian densities: The generalized quasi-topological term S is an appropriate choice for all base manifolds besides T 2 × S 2 , while the quasi-topological term Z is an appropriate choice for all base manifolds besides S 2 × S 2 . The complete action we consider is then where we have allowed for the possible contribution of the Gauss-Bonnet term. In this case, the AdS 6 vacua of the theory are characterized by being solutions to h(f ∞ ) = 0, where a definition that will turn out to be useful later on. The field equations can be obtained from the reduced action or computed directly. As anticipated, when we insert the single-function Taub where C B is an integration constant. The explicit form of the various terms appearing in the field equation is the following. The Einstein gravity contributions to the field equation can be expressed in the form where κ is defined by (3.7) Next are the Gauss-Bonnet contributions, which for the various base spaces read The quartic contributions to the field equations are, of course, more complicated. The ones due to the generalized quasi-topological term read

JHEP10(2018)095
where E B is a base-dependent contribution, which takes the explicit form On the other hand, the quartic quasi-topological contributions yield where now the base-dependent factors E It should be emphasized that while for static and spherically symmetric solutions the quartic quasi-topological term yields algebraic field equations [22], these become non-linear second-order differential equations for Taub-NUT metrics. This is an interesting difference with respect to the Gauss-Bonnet case, for which the equations determining the metric function are algebraic for both kinds of solutions.

Einstein gravity
Just like in the four-dimensional case, in the following subsections we will be studying the different base spaces independently. As before, it is illuminating to start with a quick study JHEP10(2018)095 of the situation for Einstein gravity, for which the analysis can be performed at the same time for all base spaces. Indeed, if we set λ GB = ξ = ζ = 0, (3.5) can be easily solved for V B (r). Imposing the NUT condition V B (r = n) = 0 first, one is left with where we set the integration constant The regularity condition (1.5) imposes Hence, we find β τ = 12πn for B = CP 2 and B = S 2 × S 2 , β τ = 24πn for B = S 2 × T 2 , and β τ = ∞ for B = T 2 × T 2 , which forbids the existence of regular solutions with a compact S 1 in that case. If we impose the bolt condition V B (r = r b ) = 0 instead, we find where in this case we related C B to n and r b through Finally, the regularity condition (1.5) produces the following relation between r b , n and the period of τ , (3.25) Just like in D = 4, we must require the quantity inside the square root to be positive and, of course, r b > 0, which in each case restricts the values of n for which solutions exist. Besides the regularity condition (1.5), additional constraints on β τ arise both for NUTs and bolts when demanding the absence of Misner string singularities -see e.g., discussion in [78]. For example, for B = CP 2 , we must demand β τ = 12πn. Just like in D = 4, for the Einstein gravity NUT this condition is automatically implemented by (1.5). This is not the case in general, and the conditions must be imposed separately.

B = CP 2
Let us now turn on again the higher-curvature couplings in (3.3). The first base space we consider is CP 2 . For this, we can write

JHEP10(2018)095
where the coordinate ranges 10 are 0 ≤ ξ 1,2 ≤ π/2 and 0 ≤ ψ 1,2 ≤ 2π. Now, we consider the metric asymptotically, making the rescalings τ → 6nψ and r → r/ √ 6. This gives at large r. For the specific case of 6n 2 f ∞ /L 2 = 1 , this boundary metric is just that of a round S 5 provided that the coordinate ψ has period 2π to ensure regularity. In all other cases, it is the metric of a squashed sphere [64] and, in analogy with the D = 4 case, it is customary to parametrize such squashing with the parameter α, defined in terms of n through 6n 2 f ∞ /L 2 = 1/(1 + α) .
We begin our study of Taub-NUT/bolt solutions with this base space by considering the asymptotic behaviour of the metric. The asymptotic solution for V CP 2 (r) consists of a particular and homogeneous solution. The particular solution is found by expanding V CP 2 (r) in a 1/r series and solving the field equations to determine the constants order by order. The result is where h (f ∞ ) denotes the derivative of h(f ∞ ) -see (3.4) -with respect to f ∞ . To obtain the form of the homogeneous equation, we again write V (r) = V p (r) + g(r), and work to linear order in g(r). While both S and Z contribute in the same way to the particular solution, the contributions differ in the homogeneous equation. The resulting equation, in the limit of large r, takes the form: where a(r) and b(r) are the leading terms in this expansion, taking the explicit forms and we have definedξ to simplify the presentation of the terms above. We recognize that the contributions in parenthesis in a(r) and b(r) are directly related to the squashing parameter and vanish when the base is a round sphere; in that case, the solution reduces to just pure AdS.

JHEP10(2018)095
In the limit of large r, the homogeneous equation can be solved in terms of special functions. First, when ξ = 0 the homogeneous solution reads while if ξ = 0 it takes the form: where I ν (x) and K ν (x) are the modified Bessel functions of the first and second kinds, respectively. The explicit form of these solutions is not as important as what their asymptotic behaviour tells us: because c(r) > 0 by virtue of demanding the graviton is not a ghost, when a(r) < 0, in each case the homogeneous solution consists of a super-exponentially growing and super-exponentially decaying part. When a(r) > 0, both of the above solutions oscillate more and more rapidly near infinity and are ultimately pathological. Therefore, to ensure that the solutions are physically reasonable, we must demand that a(r) < 0, while also requiring h (f ∞ ) < 0. The general solutions to these constraints with λ GB = 0 are a bit messy, and so we quote the result explicitly only in the case λ GB = 0. In that case it is a straight-forward matter to show that these conditions are satisfied -independent of n -provided that Interestingly, in contrast to the four dimensional case, here the mass parameter M does not enter into the constraints, with the result that there is no pathology associated with the negative mass solutions (see below). Furthermore, note that for ζ non-zero, simply demanding ξ < 0 is not enough since one must also require that h (f ∞ ) < 0 -this is the origin of the more complicated constraint in that case. It can be shown that, for λ GB = 0, ξ = 27/256 − ζ corresponds to the critical limit of the theory, which has f ∞ = 4/3. A consequence of these bounds on the coupling is that, when one considers a theory that contains only a single one of the quartic terms, then it is not possible to reach the critical limit of the theory at physical coupling. This situation is similar to what happens to cubic GQTG for spherically symmetric black hole solutions in D ≥ 6.

Taub-NUT solutions
We now consider NUT solutions where V CP 2 (r = n) = 0. Further restrictions on V CP 2 (r) arise due to regularity of the metric. Recall from the discussion above that the boundary is a squashed S 5 . Regularity of this boundary metric requires that ψ := τ /(6n) has period 2π, which in turn means τ ∼ τ + 12πn. A further constraint is imposed on the derivative of V CP 2 (r) near the NUT where the absence of conical singularities at a zero of V CP 2 requires that τ is periodic with period β τ given by β τ = 4π/V CP 2 (r = n). Consistency of these two regularity conditions fixes β τ = 12πn and so we therefore have the following series JHEP10(2018)095 expansion near the NUT: Substituting this expression into the field equations, and expanding in (r − n), we find 16 3 where we have conveniently redefined the integration constant C CP 2 = −4GM/(9π), where M will correspond to the ADM mass of the solution. The first condition (shown above explicitly) determines M in terms of the couplings and the NUT charge, and the next two relations are automatically satisfied. The next non-trivial relation is linear in a 3 , allowing one to solve for a 3 as a function of the free parameter a 2 . This trend continues to higher order in the field equations, and thus there is a single free parameter that is left unfixed by the field equations and regularity conditions. This is fully analogous to the D = 4 case. The near horizon solution can be joined to the asymptotic solution that was presented above by numerically integrating the field equations. The near horizon expansion is used as initial data, with the shooting method employed to determine the free parameter a 2 . A careful choice of this parameter is required to ensure the growing modes present in the asymptotic solution are not excited. Ensuring this, we find a unique a 2 for which the solution can be integrated, with the result for several values of the coupling shown in figure 9. For comparison, the Einstein gravity solution is shown in red and we see that the solutions to the higher curvature theories are qualitatively the same with the main difference being that they approach f ∞ r 2 /L 2 with f ∞ depending on the value of the couplings. We also note that, while the top left and bottom left plots depict solutions with positive mass, the top right and bottom right plots depict solutions with negative mass. The fact that the negative mass solutions can be constructed and possess no inherent pathology is in contrast with the four dimensional case, where the negative mass solutions possessed pathological asymptotic structure.
Let us now turn to the free energy of the NUTs and compute the regularized on-shell action for these solutions. With minor modifications, the prescription (1.6) introduced in [27] can be used to eliminate the divergent terms in the on-shell action. The Euclidean action, completed with the generalized boundary term and counterterms is given by The evaluation is facilitated via the asymptotic expansion presented above and the expansions near r = n in the NUT case or r = r b for the bolts. Near the boundary, the bulk action has several divergent components that are precisely canceled by the generalized boundary and counterterms. Note the addition of a new counterterm, nonproportional to a * , on the last line above. This appears because, strictly speaking, the spacetime is not asymptotically AdS -the boundary is not maximally symmetric except for the choice of NUT parameter that yields the undeformed five sphere. The additional counterterm was chosen since it vanishes identically when the boundary is maximally symmetric (and so could be dropped in those cases) but allows for the cancellation of the linear divergence in the case of B = CP 2 considered here. Just like the four-dimensional case discussed in detail in appendix B, eliminating the divergent terms also removes all possible constant terms coming from boundary contributions, leaving us with the bulk action evaluated at r = n, and nothing else. The final result is from which the total energy and entropy can be found to be, 40) and the first law dE = T dS is verified to hold.

JHEP10(2018)095
Similar to the discussion for the S 2 base in the case of ECG, here we can also enlarge the thermodynamic phase space and construct the extended first law. The expression is slightly more complicated, reading where we have again restored the dimensionality to the coupling constants. The potentials appearing in the extended first law read (3.42) The expression above for the thermodynamic volume holds also in Einstein gravity though there appear to be no previous computations of this quantity in the literature for higher dimensional Taub-NUT solutions. It is noteworthy that the thermodynamic volume here is positive, while the thermodynamic volume is negative in the D = 4 case. Of course, we find that the Smarr relation consistent with scaling is satisfied by the thermodynamic quantities defined above: Note that if we turn off the quartic couplings, then the result for the free energy reduces to that previously calculated for Einstein gravity and Gauss-Bonnet gravity in [34,39] up to an overall factor of 8/9, the same discrepancy noted in [64]. We have carefully revisited the calculations in [34,39] and have traced the discrepancy to the ratio of volumes of S 2 ×S 2 to CP 2 . In [34] it is claimed that the thermodynamic quantities for both base spaces are identical. However, we have found this to be true only up to an overall ratio of the volumes of the base spaces. For CP 2 normalized so that R ab = g ab the volume is Vol CP 2 = 18π 2 , while the volume of S 2 × S 2 is given by Vol S 2 × S 2 = (4π) 2 . The ratio of these volumes is precisely 8/9, which accounts for the observed discrepancy. 11 In figure 10 we show plots of the Euclidean on-shell action for the NUT solutions with λ GB = 0. As is clear from eq. (3.39), this depends on the higher curvature couplings only through the combination ξ + ζ. In each case, there is only a single branch and, from the figure, we see that its qualitative structure depends on whether ξ +ζ is positive or negative. For consistency with the plots presented earlier in the document, we have indicated regions of negative mass with dashed curves. However, unlike the four dimensional case, there is no pathology associated with the negative mass solutions for the quartic theories in six dimensions. The region of negative mass solutions shrinks and eventually vanishes as ξ + ζ → 27/256, which corresponds to the critical limit of the theory. 11 In [34] a different set of coordinates is used, but the metric of CP 2 is still normalized so that R ab = g ab .
Using the fact that the coordinates in [34] have ranges 0 ≤ u ≤ ∞, 0 ≤ θ ≤ π, 0 ≤ φ ≤ 2π and 0 ≤ ψ ≤ 4π we obtain the result presented here with the correct overall factor. In higher dimensions, Vol CP k = 2 k (k + 1) k π k /k! and the volume of k 2-spheres is Vol S 2 × · · · × S 2 = (4π) k . In higher dimensions, we find the ratio between thermodynamic quantities for the two bases is 2 k k!/(k + 1) k .

Taub-bolt solutions
We now consider Taub-bolt solutions which satisfy V CP 2 (r b ) = 0 for r b > n. In this section, we turn off the Gauss-Bonnet coupling to limit the size of the parameter space. Regularity demands that V CP 2 (r b ) = 1/(3n), and therefore we write the near horizon expansion as Substituting this expansion into the field equations and solving order by order in (r−r b ), we find the first two relations fix the integration constant C CP 2 and the relationship between r b and n: where, just as in the NUT case, we have set C CP 2 = −4GM/(9π).
Let us now examine the second relation above in more detail. When the higher curvature terms are turned off, the bolt radius is given by (3.25) with κ = −1, i.e.,

JHEP10(2018)095
Since r b must be real and larger than n, we must then have n < √ 15(2 − √ 2)L/30. As in the four-dimensional case, there is a maximum value of n for bolts in Einstein gravity. In particular, this means that there does not exist a bolt solution near the undeformed five sphere, for which n = L/ √ 6. Of course, the behaviour is different with higher curvature corrections, but there are some notable differences from what was observed in the four dimensional case.
Depending on the relative size of ξ and ζ, the behaviour of r b as a function of n can either resemble that of Einstein gravity (namely, there is a largest value of n for which a bolt exists) or resemble that observed in the four dimensional cubic case discussed earlier in this paper (bolts exist for arbitrarily large n). In fact, this classification is completely determined by the sign of the quantity ζ − 6ξ. If this quantity is positive, there exists a maximum value of n; if this quantity is negative, bolts exist for arbitrarily large n. 12 From the perspective of the phase structure of the bolts, the most interesting scenario occurs when there are three values of r b for a particular n -one would expect these cases could yield swallowtail type behaviour and critical phenomena. We can constrain the regions of parameter space where three bolts exist by searching for 'critical points'. More specifically, such a critical point would occur when ∂n/∂r b = ∂ 2 n/∂r b 2 = 0, while respecting eq. (3.46). These points will mark transitions in the maximum number of bolts for given couplings. We were unable to solve the resulting constraints analytically, but it is straightforward to do so numerically. This results in the breakdown of parameter space shown in figure 11. It is useful to understand the qualitative behaviour of the bolts in the various partitions of the parameter space shown in figure 11. We illustrate this in the top row of figure 12, which represents a 'vertical slice' through figure 11 for ξ = −4 × 10 −5 . The plot on the right is a zoomed-in copy of the left, and the decreasing opacity of the blue curves (left to right) denotes ζ becoming more negative, while the red curve corresponds to the Einstein gravity result when both couplings vanish. We see that when ξ and ζ are small (or, equivalently, when r b is large) the bolt radius reduces nicely to the Einstein gravity result. The interesting behaviour is observed for smaller bolt radius. The first curve corresponds to ζ = −20 × 10 −5 which is in the white region of figure 11 and above the orange line. We see that in this case, the behaviour is similar to Einstein gravity, with two possible values for the bolt radius. As ζ is further decreased, the structure of the curve remains similar but a small 'flattened' region begins to form, ultimately becoming vertical for ζ ≈ −23.16 × 10 −5 which corresponds to the point on the orange line of figure 11. Continuing to decrease ζ further, we see that a bump emerges, and as a result there are up to four values of r b for a given n. This behaviour continues until ζ < 6ξ, which corresponds to the blue shaded region of figure 11. At this point, the structure of the curve changes drastically, and there are bolts for arbitrarily large n. Further, in the region where ζ < 6ξ but remains above the region bounded by the black curve in figure 11, there are up to three bolts for a given  Here the orange and black curves denote lines of 'critical points', i.e., for a given NUT charge, three solutions for r b coalesce. The red dot represents the single point in the physical parameter space where there is a coalescence of four roots. Within the blue shaded region, ζ < 6ξ and there are bolts for arbitrarily large n. In the complement, ζ > 6ξ and there is a largest value of n for which bolts exist. On the black locus of critical points, the critical point is always physical (i.e. of lowest free energy) within the blue region, otherwise (for the dashed portion of the curve) the situation can depend on which branch of the cusp minimizes the free energy, and also on whether or not there are re-entrant phase transitions as described in the text. value of n. As ζ is further decreased we continue to see three bolts for a given n until we reach the black curve of figure 11, which corresponds to ζ ≈ 38.03 × 10 −5 . At this point, the three bolts coalesce, and for values of ζ smaller than this there is only ever a single bolt for a given n.
It is also possible for ζ to take on positive values, provided that ξ < 0 and ζ ≤ 27/256 − ξ. The bottom row of plots in figure 12 shows representative behaviour in this case. The qualitative shape of the curve is controlled by the ratio ξ/ζ. When ξ/ζ → 0 − , a peak forms at small n. The overall behaviour is similar to Einstein gravity: there is a maximum value of n beyond which bolts cannot exist. For n smaller than this value, there are two values of r b for any given n.
The above discussion highlights the general trend in this parameter space. The lines of 'critical points' mark the boundaries where there is a change in the maximum number of bolts for a given NUT charge. For a fixed ξ, the structure is (referring to figure 11): two bolts and Einstein-like structure in the white region above the orange line; up to four bolts in the white region below the orange line; up to three bolts in the blue shaded region above the black line, and one bolt in the blue shaded region below the black line. When ζ takes on positive values, the structure remains the same as in the white region above the orange line, but a peak forms at small n as ξ/ζ → 0 − . So far, our study of the bolt solutions has focused on the properties of the near horizon solutions. It is important to verify that these near horizon solutions can be joined smoothly on to the asymptotic solution (3.29) that was presented at the beginning of this section. This can be shown by numerically solving the field equations, with some relevant examples shown in figure 13. The left plot shows example bolt solutions for n/L = 7/100. In this regime, both Einstein gravity and the quartic theories admit bolt solutions, and the two can be compared. The solutions are qualitatively similar but, of course, the solutions to the quartic theories asymptote to f ∞ r 2 /L 2 with f ∞ = 1. In the right plot, we show examples for n/L = 2 -for this value of the NUT parameter, there are no bolt solutions in Einstein gravity.
Finally, turning to the on-shell action, it can be computed using the same prescription as in the NUT case, but now evaluating for the bolt at r = r b . In performing the calculation, we make use of the near horizon equation (3.45) to simplify the result. We find that, We can study the extended thermodynamics of these bolts in the same manner as the NUTs. The extended first law has the same form as (3.41) but for the bolts the potentials are given by

JHEP10(2018)095
and we recall that here we are working with λ GB = 0. These quantities also satisfy the Smarr relation that follows from scaling, which has the same form here as in (3.43). Again, the formula for the thermodynamic volume is unaltered from its form in Einstein gravity. Though, since r b implicitly depends on the higher-curvature couplings, the numerical value of the thermodynamic volume for fixed n, ξ and ζ will in general differ from the Einstein gravity value. Contrast this with the situation for the NUTs where the thermodynamic volume is completely insensitive to the theory of gravity, so long as the theory belongs to the generalized quasi-topological class. The Euclidean on-shell action exhibits rich structure for the bolt solutions. In understanding the behaviour, it is helpful to once again refer to figure 11. As it turns out, this figure partitions the parameter space into regions where the behaviour is qualitatively similar. Referring to figure 11, the most interesting changes in behaviour occur when the orange and black lines are crossed, which correspond to actual critical points in the thermal phase space marking the appearance/disappearance of swallowtail structures in the on-shell action. Also when transitioning from the white-shaded to blue-shaded region, the action switches from terminating at a cusp at some finite n to existing for all values of n.
In the white region, in all cases but Einstein gravity there will be a zeroth order phase transition between bolt solutions and NUT solutions at the value of n corresponding to the maximum value of r b . In Einstein gravity there is also a phase transition at this point, but in that case it is first order. As an example that highlights the salient points pertaining to the bolts, let us once again consider the ξ = −4 × Though the discussion will make reference to numerical values in only this particular case, the qualitative features are general.
Particularizing now the discussion to ξ = −4 × 10 −5 , for positive ζ through to ζ ≈ −23.1565×10 −5 (which corresponds to the orange line in figure 11), the behaviour is similar to Einstein gravity, with the on-shell action exhibiting two smooth branches that end at a cusp located at the maximum value of NUT charge. Precisely when ζ is chosen on the orange line shown in figure 11 (ζ ≈ −23.1565 × 10 −5 in this case), the upper branch of I E develops a cusp, corresponding to a critical point in the system. As ζ is further decreased, a swallowtail emerges from the cusp on the upper branch, as shown in the second plot of figure 14. Further decreasing ζ elongates the swallowtail, and eventually it intersects the lower branch of I E -for the particular case of ξ = −4 × 10 −5 , this intersection occurs for ζ ≈ −23.705 × 10 −5 . This intersection then gives rise to a region where a re-entrant phase transition occurs as n is increased, as shown in the center-left plot of figure 14. The two vertical black, dotted lines show the locations where these transitions occur. There is a first order phase transition from phase 1 to phase 2, followed by a zeroth order phase transition which returns the system back to the initial phase. It is in this sense that we have a re-entrant phase transition -a monotonous variation of the NUT charge gives rise to two phase transitions with the final and initial phases coinciding. Let us note that re-entrant phase transitions were first observed in nicotine/water mixtures in [80]. In the context of black hole physics, while somewhat exotic, they are well-established -see [81] JHEP10(2018)095 for an example in a rotating black hole spacetime, and [82,83] for black hole examples in higher curvature theories of gravity. We believe this is the first instance observed for NUT charged solutions. As ζ is further decreased, the swallowtail continues to elongate, and for ζ ≈ −23.753 × 10 −5 , the tip of the swallowtail extends past the cusp -this ends the region of parameter space for which re-entrant phase transitions occur.
There is a drastic change in structure at ζ = 6ξ. Corresponding to the boundary of the blue-shaded region in figure 11, this condition yields the largest ζ for which there is a maximum NUT parameter for which bolts exist. From the perspective of the onshell action, essentially what happens is, at this point, the swallowtail has now elongated "to infinity". Between ζ = 6ξ and ζ ≈ −38.026 × 10 −5 the action displays a swallowtail structure that is associated with a first-order phase transition. The swallowtail vanishes at a critical point when ζ ≈ −38.026×10 −5 (the black line in figure 11). For ζ −38.026×10 −5 , the on-shell action displays only a single branch for all values of n.
Lastly, let us make some remarks regarding the critical points that are present at some points in the parameter space. As mentioned above, the lines of critical points appearing in figure 11 are bonafide critical points in the thermodynamic parameter space. When in the region of the parameter space corresponding to the white region of figure 11, the action exhibits a cusp structure qualitatively similar to that shown in the top left plot of figure 14. We find that one of the critical points always occurs on the upper branch of this cusp (those corresponding to the orange curve in figure 11). These critical points will, therefore, not be realized since they do not comprise the dominant contribution to the partition function. The critical points that correspond to the points on the black curve shown in figure 11 belong to the lower branch of the cusp in the white region or are on the single physical branch in the blue shaded region. These critical points are physically realized.
At the critical point, certain physical quantities blow up in power law fashion. To get a sense of the critical exponents governing these divergences, we can study the behaviour of the specific heat, where F = T I E andα is the critical exponent governing this divergence. 13 Due to the complexity of the equations relating the bolt radius to the NUT parameter, it is difficult to perform an analytic study near the critical point. Instead, to make progress, we plot numerically and extract the slope of this line via a linear fit. As an example, we find in the case ξ = −4 × 10 −5 the following fit:

JHEP10(2018)095
which after some simple algebra yieldsα = 0.659 (3.54) which is consistent withα = 2/3 to within the numerical precision. The valueα = 2/3 for the critical exponent is often observed for the divergence of the specific heat at constant pressure in black hole systems -see, e.g., [85]. In this sense, it is not surprising to find that the same critical exponent governs the behaviour near the critical point for the bolts. A numerical survey of many critical points for different values of the couplings shows that they are all consistent with this result. The red dot shown in figure 11 represents a special point in the parameter space where two critical points merge. Because of this, one might hope to see novel critical exponents similar to how the coalescence of multiple critical points leads to non-mean field theory critical exponents for Lovelock black holes [86]. However, unfortunately, this is not the case here. The reason is that as the red dot is approached, there is one critical point on the upper branch of the cusp and one on the lower branch. When these critical points merge, they also meet at the cusp which acts as a phase boundary -no solutions exist beyond the tip of the cusp. To within the accuracy of our calculation, the critical exponent associated with each critical point as the cusp is approached remains consistent withα = 2/3.

B = S 2 × S 2
As mentioned at the beginning of this section, the case of B = S 2 × S 2 is somewhat special as it only admits theories of the generalized quasi-topological type. Our action is then (3.3) but now with ζ = 0.
The metric takes the form of eq. (1.3), with the following 1-form and base space metric, A 2 S 2 ×S 2 = 2 cos θ 1 dφ 1 + 2 cos θ 2 dφ 2 , dσ 2 S 2 ×S 2 = dθ 2 1 + sin 2 θ 1 dφ 2 1 + dθ 2 2 + sin 2 θ 2 dφ 2 2 . (3.55) The field equations are relatively complicated and have been listed in full detail above. The S 2 × S 2 case is somewhat interesting because there is the appearance of a logarithmic term in the integrated field equations. We first wish to consider the field equations asymptotically. The vacua of the theory are determined by the equation Again, the solution will consist of a homogenous part g(r) and a particular part V p (r), given by and −60 × 10 −5 (more to less opacity, respectively). Along the first curve, there is a critical point located at n/L ≈ 0.07815, while the other two curves are smooth. The structure of the on-shell action is qualitatively similar to these last two curves for all ζ −38.026 × 10 −5 .

JHEP10(2018)095
and obtained by performing a large-r series solution of the field equation. The homogeneous equation, at large r, is a second order differential equation for g(r), a(r)g (r) + b(r)g (r) + c(r)g(r) = 0 (3.58) with the coefficients given by, (3.59) The differential equation can be solved exactly in terms of modified Bessel functions, The asymptotic form of this solution will consist of a super-exponentially growing and a super-exponentially decaying mode provided that a(r) < 0. If a(r) > 0, the asymptotic behaviour is pathological. Therefore, we must demand that a(r) < 0 and use the asymptotic AdS boundary conditions to set the growing mode (C 1 ) to zero. Demanding a(r) < 0 is equivalent to demanding that the coupling ξ is negative, since the term in parentheses in the expression for a(r) is always positive.

Taub-NUT solutions
Taub-NUT solutions with base space S 2 × S 2 are generically pathological, even in Einstein gravity, due to a curvature singularity at the NUT r = n. The situation is actually worse here, since the logarithm in the field equations -see above -develops an essential singularity at r = n. Thus it seems that there are no well-behaved NUT solutions for B = S 2 × S 2 .

Taub-bolt solutions
Even in the case of Einstein gravity, Taub-NUT solutions are singular at the NUT when the base is either B = S 2 × S 2 or B = S 2 × T 2 . This pathology leads to problems when solving the field equations for higher-curvature gravities. These problems were first observed in the case of Gauss-Bonnet gravity [30] where it was found that NUT solutions with B = S 2 × S 2 do not exist for non-vanishing Gauss-Bonnet coupling. It was conjectured that this was due to the fact that the corresponding Einstein gravity solutions are singular [30].
In the higher-curvature gravities considered here, we are faced with a similar problem -the field equations develop an essential singularity at r = n, and therefore the NUT solutions do no exist. However, luckily, the situation is better in the case of bolt solutions. In this case, we will proceed in the same manner as before noting that the periodicity of Euclidean time must be τ ∼ τ + 12πn to ensure the absence of Dirac-Misner string singularities. This value is the same as the periodicity enforced by the field equations themselves in Einstein gravity. In a higher curvature theory, like the ones studied here, JHEP10(2018)095 the field equations may not naturally ensure regularity and we must enforce it by hand. This approach was used in [30] to successfully construct bolt solutions in Gauss-Bonnet gravity that limit smoothly to the Einstein gravity solutions as the Gauss-Bonnet coupling is turned off, and we will adopt the same approach here.
We begin by expanding the metric near the bolt The first two relationships are which correspond to the O(1) and O (r b − n) 1 terms in the field equations. At the next order, one determines a 3 in terms of a 2 , n, and r b . This pattern continues, and one is left with the unfixed parameter a 2 that should be, as usual, determined by demanding that the field equations have the correct asymptotic behaviour. The on-shell action can be computed in the same way as before and we find that the expression is given by The standard thermodynamic relationships are now where the second equation in (3.62) was used to simplify the results.

JHEP10(2018)095
Next we consider the metric (1.3) with B = S 2 × T 2 . The potential and base space metric read We take the action to be (3.3) but now with ξ = 0. The vacua of the theory are then determined as roots of the polynomial equation, At large distances, the solution consists of a particular and homogenous part, V (r) = V p (r) + g(r). The particular solution is given, as usual, by performing a series expansion of the equation in 1/r and matching the coefficients to give a consistent solution order by order. This yields for the particular solution. The homogenous equation, at large r, is a second order differential equation for g(r), a(r)g (r) + b(r)g (r) + c(r)g(r) = 0 (3.68) with the coefficients given by, (3.69) The large-r homogeneous solution has a known solution in terms of modified Bessel functions, The asymptotic form of this solution will consist of a super-exponentially growing and a super-exponentially decaying mode provided that a(r) < 0. If a(r) > 0, the asymptotic behaviour is pathological. Therefore, we must demand that a(r) < 0 and use the asymptotic AdS boundary conditions to set the growing mode (C 1 ) to zero. Demanding a(r) < 0 is equivalent to demanding that the coupling ζ is negative, since the term in parentheses in the expression for a(r) is always positive.

Taub-NUT solutions
First, we note that whenever the base is S 2 × T 2 , the NUT solutions will necessarily possess curvature singularities at r = n. With this in mind, let us begin by discussing the situation for Einstein and Gauss-Bonnet gravity, as subtleties arise already in these cases. In the pure Einstein gravity case, the solution to the field equations reads (3.71) where C S 2 ×T 2 is the integration constant in the field equation (3.5). If we demand that this solution permits NUTs, then we must have V (r = n) = 0 and this fixes the integration constant uniquely as which yields (3.20) with κ = 0. Alternatively, we could have begun by expanding V (r) as a series near r = n as Then we would find that, near the NUT, the field equation takes the form, From either perspective, it is easy to check that, at the location of the NUT we have V (n) = 4πT = 1/(6n). Next we turn on the Gauss-Bonnet term. The addition of this term gives rise to two exact (but messy) solutions for V (r). One of these solutions limits to the Einstein gravity result as λ GB → 0 (this solution is given in [30] -see eq. (6) in that paper), and the other blows up in that limit. If we expand the field equations near r = n using the same ansatz as before, we now find the first few terms to be This expansion of the field equations is the same as (3.74) in the limit λ GB → 0. However, we can see that the solution of these equations for the temperature does not limit to the Einstein gravity result. While the result for the integration constant is the same as before, the Gauss-Bonnet term now contributes at a lower order in the field equations, bringing a new term into the expansion that is not present in the case of pure Einstein gravity -see the O((r − n) 2 ) term above. There are three possibilities for a solution at that order: (1) λ GB = 0, in which case we recover the Einstein result quoted above, (2) T = 0, which when extended to the full non-perturbative solution corresponds to the Gauss-Bonnet branch that limits to the Einstein branch (the one discussed in [30]), or (3) 4πT = 1/(3n), which

JHEP10(2018)095
corresponds to the Gauss-Bonnet branch that is singular in the limit λ GB → 0. We conclude that none of these possibilities for the temperature actually limits to the temperature of the Einstein gravity solution, even though the full non-perturbative solution corresponding to the T = 0 branch does limit to the Einstein gravity solution.
What we have arrived at here is an order of limits problem: performing the λ GB → 0 before the r → n limit gives a different result than first performing the r → n limit followed by λ GB → 0. The limit of the temperature expression is not continuous: for any non-zero λ GB it should be T = 0, but when λ GB is precisely zero, it "jumps" to 4πT = 1/(6n). The origin of this incompatibility of limits would seem to be linked to the fact that the space is actually singular at r = n for B = S 2 × T 2 .
Having reviewed this structure, let us now consider the case with both ζ and λ GB nonvanishing. We expand the metric function as above and demand the field equations are satisfied order by order. The first terms in the expansion of the field equations are and we see here that the new higher curvature terms contribute at even lower order than the Gauss-Bonnet term. Again, we have two (non-trivial) possibilities for a solution here. The first possibility is to have T = 0. This case leads to a non-perturbative solution with the value of a 2 determined as a solution to the following equation: with a 3 , a 4 and so on given directly (and uniquely) as functions of the couplings and a 2 . From this equation, it is clear that one of the three possible roots for a 2 will limit to the extremal Gauss-Bonnet solution when ζ → 0. The second possibility is for 78) and the rest of the constants are determined uniquely in terms of T , the couplings, and n. In this solution, T is left as an arbitrary parameter and the solution limits to neither the Einstein result nor the Gauss-Bonnet case as ζ → 0. Based on our intuition from the Gauss-Bonnet situation analyzed above, the most reasonable conclusion would seem to be that the NUT solutions on B = S 2 × T 2 should be regarded as extremal solutions. That is, taking T = 0 seems to be the most reasonable of the various options discussed above.

Taub-bolt solutions
For simplicity, we will at this point set λ GB = 0. Then, we expand the metric function as

JHEP10(2018)095
and demand the field equations are satisfied order by order. There will be no curvature singularities unless r b = n. Since the NUT solutions for B = S 2 × T 2 are somewhat pathological, and should probably be regarded as zero temperature or extremal solutions, there is not a natural periodicity enforced on the bolts via regularity of the NUTs. If we wish the temperature to match the bolts from Einstein gravity, then we would require that β = 24πn. Here, similar to what was done in [30] for this base, we will keep the temperature explicitly present in the following analysis.
The first few components of the field equations are given by The first equation determines the mass parameter, while the second gives a relation between the bolt radius and the NUT parameter. At the next order in the field equations, a 3 is determined by a 2 , the bolt radius and the NUT parameter. This pattern continues in the usual way, and once a 2 is determined by demanding the appropriate asymptotics the full near horizon solution will be determined. We can compute the free energy from the on-shell Euclidean action as before. Denoting the compactification length of the T 2 as l, we find the following: from which we obtain for the energy and entropy.

B = T 2 × T 2
Let us finally consider the case of T 2 × T 2 . In that case, the 1-form and the metric of the base space are given by Each pair of coordinates (η 1 , ζ 1 ) and (η 2 , ζ 2 ) parametrize a T 2 , and for simplicity we can assume that η 1 , ζ 1 have both periodicity l 1 and η 2 , ζ 2 have periodicity l 2 . The coordinate τ

JHEP10(2018)095
is also compact and its period β τ is a free parameter, since there are no regularity conditions in this case. The function V (r) is determined as usual by eq. (3.5). Both densities Z and S can be introduced in this case, and for completeness we will also consider the Gauss-Bonnet term. As usual, let us start by determining the asymptotic behaviour of the solution, which we decompose as a particular solution plus a homogeneous one as V (r) = V p (r)+r 2 g(r)/L 2 . The particular solution is found by performing a 1/r expansion, which reads On the other hand, at linear order g(r) satisfies the equation where the functions a(r), b(r) and c(r) take the following form asymptotically The solution is qualitatively different depending on whether ξ = 0 or ξ = 0. In the latter case, taking into account only the leading terms in a, b c, we get the following solution asymptotically: On the other hand, when ξ = 0, the asymptotic solution reads instead (3.90) In these expressions, J k and Y k are Bessel functions of the first and second kinds respectively. Let us note that when the argument of the Bessel functions is real, their asymptotic behaviour is oscillatory, while for imaginary arguments they behave as real growing and decaying exponentials, which is the kind of behaviour we require in order to impose an appropriate boundary condition asymptotically. Noting that h (f ∞ ) < 0 for the physical vacuum, we find the condition ξ + δ ξ,0 ζ ≤ 0, where δ ξ,0 = 1 if ξ = 0, and 0 otherwise. Then, by demanding that the solution has the right asymptotic behaviour, we fix one of the integration constants, as usual. From (3.84), we see that the boundary metric when r → ∞ takes the form (3.91)

JHEP10(2018)095
We also can determine the ADM energy, which is proportional to the coefficient of the 1/r 3 term in (3.84), Let us now explore the regularity conditions which as usual will fix the remaining integration constant of the solution.

Taub-NUT solutions
We start with NUT solutions which, as usual, are determined by the conditions V T 2 ×T 2 (n) = 0 and V T 2 ×T 2 (n) = 4πT . In the case of Einstein-Gauss-Bonnet gravity the solutions are forced to be extremal, i.e., T = 0 -see (3.22) for the Einstein case. Here we find that this is also the only possibility if we want the solutions to reduce to those of EGB. Setting T = 0 and plugging an expansion of V T 2 ×T 2 (r) around r = n into the field equations, we obtain and an infinite number of equations that fix the coefficients a i as functions of the couplings λ GB , ξ, ζ. For example, the first two of them read Note that in this case there are not free parameters and the full series is univocally determined once we choose one of the roots of the first equation, which should be the one that reduces to the Einstein gravity one when λ GB , ξ, ζ → 0. Therefore, at least near r = n we have solved the equation by using the series expansion. However, when ξ, ζ = 0 one also has to make sure that the asymptotic behaviour is the correct one. In this case it is unclear whether the solution constructed from the near-horizon expansion satisfies this property.
Were this not the case, it would imply that no NUT solutions exist when we include the quartic corrections. In any case, assuming the solution exists globally, it is illustrative to compute the free energy as F = T I E , the result being On the other hand, the ADM energy (3.92) reads As we can see, n acts as a thermodynamical variable and it has to be taken into account when we compute the energy E from the free energy F . This is not a coincidence and, as we show below, the same observation is valid for bolt solutions.

Taub-bolt solutions
Let us turn now to bolt solutions, which are more interesting. As usual, we assume V T 2 ×T 2 (r b ) = 0 and V T 2 ×T 2 (r b ) = 4πT for certain r b > n. Then, we can Taylor expand the solution as 3.98) and the equations of motion fix the value of the integration constant C T 2 ×T 2 , the relation between r b , T and n, and all coefficients of the expansion a i>2 in terms of a 2 . As usual, this constant is determined by the boundary condition at infinity and can be found, along with the full V T 2 ×T 2 (r), using numerical methods. We will focus on the thermodynamic properties, which can be obtained analytically. The bolt radius r b is determined implicitly by n and T through the equation (3.99) while the integration constant C T 2 ×T 2 is given by (3.100) from which we can obtain the ADM energy of the solution using (3.92). Observe that while the Gauss-Bonnet term does not modify (3.99) or (3.100) with respect to the corresponding Einstein gravity expressions, the quartic theories produce important modifications.
The free energy can be computed from the Euclidean action analogously to the rest of cases, and it reads (3.101) The analysis now follows the same lines as in section 2.2: since we have an extended phase space, we must introduce the variable θ = 1/n, and the free energy should be interpreted as a function of the physical temperatureT = T / √ f ∞ and of θ: F = F (T , θ). From this we obtain the entropy S and the potential Ψ as defined in (2.57) (but with respect toT ): To ensure consistency of the Smarr relation, we consider Λ and the couplings as thermodynamic parameters and find that the following potentials satisfy the extended first law , With these thermodynamic potentials, the Smarr formula that follows from scaling holds -this is of the same form as eq. (3.43), but now we must include an additional 4θΨ term. Even though it is known the first law should hold in general theories, it is remarkable that this can be explicitly checked in these very non-trivial examples.

Concluding remarks
In this paper we have constructed new Taub-NUT and Bolt solutions for several highercurvature gravities for various base spaces in D = 4 and D = 6. In particular, the solutions constructed in section 2 for Einsteinian cubic gravity are the first examples of four-dimensional higher-curvature generalizations of the Einstein gravity Taub solutions. In all cases, the solutions generalize the Einstein gravity (and Gauss-Bonnet) ones, and reduce to them as the higher-curvature couplings are set to zero. Also in all cases, and analogous to the new classes of black holes constructed in [11][12][13][14][15][16][17], the solutions are always characterized by a single base-space-dependent metric function V B (r). Even though we cannot compute this function analytically in most cases, the thermodynamic properties of the solutions can be accessed in a fully analytic form, as we have shown. As we have seen, turning on the higher-curvature couplings notably modifies the structure and thermodynamic properties of the solutions with respect to the Einstein gravity case. In particular, they typically modify the multiplicity of solutions as the NUT charge is varied -e.g., new bolt solutions exist for values of n that are forbidden in Einstein gravity -and drastically modify the phase spaces -e.g., for NUT solutions with B = CP 1,2 , the free energy generically diverges as n → 0, instead of going to zero as is the case in Einstein gravity. Remarkably, the phase space of D = 6 solutions with B = CP 2 present re-entrant phase transitions, these being the first examples of this kind observed for NUT-charged solutions. It would be interesting to better understand this phase structure from the perspective of extended phase space thermodynamics.
The cases studied here (see also appendix A) are just the simplest possible within the Generalized quasi-topological family. We expect additional solutions of the form (1.3) to exist in D = 4, 6 for even higher curvature theories of this class, and similarly for D ≥ 8.

JHEP10(2018)095
It would be interesting to construct them. In particular, in D = 4 one could consider the invariants presented in [16]. It is possible that a closed form expression -valid for arbitrarily high curvature terms -for the equation determining the metric function V B (r), analogous to the one found in [16] for the black hole solutions, can be found. This would allow for a characterization of the solutions for infinitely many theories.
One of the motivations for this work, and an obvious application of our results, can be found in the holographic context. There, Taub solutions with CP D−2 2 base spaces generically dominate the semiclassical partition function for holographic theories on a particularly interesting class of squashed spheres [60][61][62][63][64]. The fact that their thermodynamic properties can be computed analytically makes our solutions particularly appealing from a holographic perspective. In particular, they can be used to study the properties of squashed-sphere partition functions for a class of theories much broader than the one available so far. As it turns out, our results here can be used to identify new universal properties, presumably valid for general CFTs (holographic or not) [87].
Although we have only considered asymptotically AdS solutions here, the corresponding asymptotically flat counterparts can be easily obtained by taking L → ∞ while keeping the dimensionful higher-curvature couplings (such as µL 4 ≡μ) finite. Asymptotically flat Taub-NUT solutions are the main constituents of Kaluza-Klein monopoles, which arise when constructing lower-dimensional solutions from compactifications of higherdimensional theories [88][89][90]. It would be interesting to explore the new solutions from this perspective.

JHEP10(2018)095
solutions with B = S 2 constructed in section 2.1 are modified by the introduction of this term. Let us then consider the Euclidean action is a particular GQTG density. Let us start by determining the AdS vacua of (A.1). As usual, we write the relation between the action scale L and the AdS radiusL asL 2 = L 2 /f ∞ . Then, the possible values of f ∞ are determined by the positive roots of the polynomial For a given vacuum, the effective gravitational constant can be computed as G eff = −G/h (f ∞ ). Hence, in order to get a positive energy graviton, we must demand h (f ∞ ) < 0, the critical case corresponding to h (f ∞ ) = 0. Just like for ECG, there is an additional constraint coming from imposing the existence of positive-energy solutions. This reads µ + 2f ∞ ξ ≥ 0 and, interestingly, it is equivalent to h (f ∞ ) ≥ 0 (assuming f ∞ > 0). Therefore, we need to identify solutions to ( All these conditions bound the space of parameters (µ, ξ), and we can write the allowed set as If the parameters belong to this set, there exists at least one AdS vacuum satisfying all the aforementioned constraints with f ∞ = 1/α, G eff = G/β. Remarkably, we do not find any other allowed vacuum, so in this region of the parameter space the vacuum exists and it is unique. In figure 15 we show the region defined by (A.4). It is convenient to divide it into three different zones. Zone A is the one with ξ ≥ 0, and in this case the allowed AdS vacua is the second largest real root of h. The largest root has h (f ∞ ) > 0 so it is not allowed. Zone B 1 corresponds to ξ < 0. There a third root appears, which becomes the largest one. This one has h (f ∞ ) < 0 but h (f ∞ ) < 0, so it is not suitable. At this point the physical vacuum is the third largest root of h. If ξ is negative enough, the two roots larger than the physical one disappear. They coalesce for (µ, ξ) = (3α 2 − 4α 3 , 3α 4 − 2α 3 ), 0 ≤ α ≤ 1/2, in which case there appears a special critical point. This line is the one which separates zones B 1 and B 2 in figure 15. Below the line, in zone B 2 , the physical vacuum is the largest root of h, which has interesting consequences for the Taub-NUT solutions, as we will see. There is a one-parameter family of critical theories, i.e., for which h (f ∞ ) = 0. We can use f ∞ to parametrize the value of the couplings in that case, namely Of course, if we impose ξ cr to be zero, we recover critical ECG, for which f ∞ = 3/2 and µ cr = 4/27.
When evaluated on the ansatz (1.3) for B = S 2 , T 2 , H 2 , we find again that the field equations of (A.1) reduce to a single equation for the function V B . As before, we find that this equation allows for an integrable factor (1 − n 2 /r 2 ), and we can write it as in (1.4), namely V 2n 2 r −2r + 2 kL 2 n 2 +r 2 −3n 4 −6n 2 r 2 +r 4 L 2 r (A.6) where k = +1, 0, −1 for B = S 2 , T 2 and H 2 , respectively.

JHEP10(2018)095
Let us start by determining the asymptotic behaviour in this case. As usual, we can separate the solution as the sum of a particular solution plus a homogeneous one. The particular solution can be obtained by performing a 1/r expansion, which yields where h (f ∞ ) = −1 + 3µf 2 ∞ + 4ξf 3 ∞ < 0, according to the unitarity constraint. From this asymptotic expansion, and using the fact that G eff = −G/h (f ∞ ) [87], we see that for a spherical base space, C = GM , where M is the ADM mass [91,92], or more appropriately, the Abbott-Deser energy [93][94][95][96]. For the remaining topologies, C is also proportional to the total energy, but the proportionality constant is different. If we now consider V (r) = V p (r)+ r 2 L 2 g(r) and expand linearly in g, we obtain the following differential equation keeping only the leading terms when r → ∞, 14 Just like for ECG, the solution is again given in terms of Airy functions, and the analysis is analogous. If Ch (f ∞ ) > 0 there is a growing mode and a decaying one, so by eliminating the former we obtain an asymptotically AdS solution. If Ch (f ∞ ) < 0 all solutions except the trivial one are pathological at infinity. Then, in order to ensure the existence of solutions of positive mass, C > 0, we demand that h (f ∞ ) > 0, which is the constraint anticipated before.
From this point on, we focus on the case B = S 2 . Then, the Taub-NUT metric takes the form (2.17), where V S 2 (r) satisfies (A.6) with k = 1. As usual, the period of τ is fixed to β τ = 8πn, which removes the Dirac-Misner string.

Taub-NUT solutions
Assuming V S 2 (n) = 0 and the regularity condition V S 2 (n) = 1/(2n), we can write an expansion around r = n as If we introduce this expansion in (A.6), we obtain a series of relations that must be satisfied order by order in (r − n). From the first one we read the mass of the solution, which is given by (A.11) 14 For example, we are neglecting a term g /r 3 against g /r.

JHEP10(2018)095
Naturally, this generalizes the ECG result (2.20) and reduces to it for ξ = 0. Also analogously to the ECG case, the following term in the expansion gives a relation between a 3 and a 2 from where we obtain a 3 (a 2 ), the next fixes a 4 (a 2 ), and so on. Therefore, once again, the complete series is determined by a single free parameter that must be chosen so that the condition B = 0 in (A.9) is met.
Let us now compute the Euclidean on-shell action of the solutions. For that, we use the generalized action (1.6), where the charge a * is given in this case by a * = (1 + 3µf 2 ∞ + 2ξf 3 ∞ )L 2 /(4G). Then, we can write the full action as (A.12) The evaluation of all terms in (A.12) is analogous to the one performed in detail for ECG in appendix B. We observe that the divergent terms coming from the various contribution cancel, and we are left with the following finite answer which generalizes the ECG result (2.25). Taking into account that β = 8πn, we can obtain the energy and the entropy E = ∂I E /∂β, S = βE − I E . The first exactly coincides with the ADM mass in (A.11), E = M , whereas for the entropy we find which generalizes the ECG answer (2.26).

Taub-bolt solutions
Let us now assume that V S 2 vanishes for some r = r b > n. In order to avoid a conical singularity we demand again that V S 2 (r b ) = 1/(2n), so that V S 2 (r) should be Taylorexpanded as Plugging this expansion into (A.6), we find that the mass of the bolt is given by

JHEP10(2018)095
Just as for ECG, the remaining equations fix the coefficients a i>2 in terms of a 2 , which must be chosen so that the solution is asymptotically AdS, a condition that selects a unique value of a 2 . The roots of (A.17) behave in different ways depending on the values of the parameters. We can characterize several qualitative features depending on the region of the parameter space shown in figure 15. First, recall that in the case of Einstein gravity, µ = ξ = 0, there are two allowed roots when n/L < (2 − √ 3)/12 1/2 -see (2.33) -and no solutions otherwise. One of the roots goes to zero for n → 0 and the other one diverges. When µ = 0 or ξ = 0 there is no root going to 0 for n → 0. In fact, in this limit we can expand r b as , c 1 = 256c 3 0 − 48c 2 0 + µ 8c 2 0 (8c 0 − 1) , (A. 18) where we must demand c 0 > 0. The first equation gives us some information about the roots, depending on the region. If ξ > 0, there is a unique value of c 0 , so there is a single solution for n → 0. Indeed, we observe that there is a unique branch in the diagram (r b , n) if ξ > 0 and that there is a solution for every value of n, including large values. When −1/16 < ξ < 0, there are two different roots c 0 , so there are two different solutions for n → 0. We see that if ξ ∈ B 1 , then these solutions extend to every n, while for ξ ∈ B 2 , the solutions only exist for n smaller than certain value. Finally, if ξ < −1/16, we find that there are no bolt solutions. In figure 16 we summarize the different possibilities. In all possible cases in which solutions exist, the mass when n → 0 is given by which can be shown to be positive as long as the parameters lie in the allowed region. Then, the Euclidean on-shell action can be computed along the same lines as in the NUT case using (A.12). The final result reads
(B.4) Then, using the asymptotic expansion (2.7) we find the boundary contribution Adding up bulk and boundary contributions, we find

C Remarks on numerical methods
We have presented in our investigation a number of numerical solutions for the NUT and bolts. Here we provide some details on how these solutions were obtained. The differential equations solved here are in general stiff, which results in difficulties in the numerical scheme. All numerical solutions presented in this work were obtained using Mathematica, utilizing the ImplicitRungeKutta method of NDSolve. This method satisfies A-stability, making it a suitable method for stiff differential equations. High WorkingPrecision was used in the numerical solver, ranging between 20 and 50 on a case by case basis. Let us make some remarks on the details of the numerical scheme, focusing on the B = CP k bases. The metric function V B (r) was expanded near a NUT or a bolt as V CP k ( ) = (k + 1)n + a 2 2 , (C.1) where = (r−n) for the NUTs or = (r−r b ) for the bolts is taken to be some small, positive quantity -typically 10 −2 L−10 −3 L in this work. The parameter a 2 is not fixed by the near horizon solution, and must be determined via the shooting method. Specifically, for a given choice of a 2 , eq. (C.1) is used to generate initial data for the differential equation, namely V CP k ( ) and V CP k ( ). Finding a numerical solution then reduces to finding a sensible value of a 2 . A generic choice of a 2 will lead to the excitation of the growing modes that appear in the asymptotic expansion of the metric function. The correct choice of a 2 will result in a numerical solution that approaches the 1/r part of the asymptotic expansion at sufficiently large r. Regardless of the choice of a 2 , the numerical scheme will eventually breakdown because of the accumulation of errors due to finite working precision. It is useful to study the point at which the numerical solution fails as a function of the shooting parameter a 2 -an example of this is shown in the left plot of figure 17. This figure makes clear that there is a special value of a 2 that allows the solution to be integrated the furthest. It also appears that this is the unique value of a 2 that joins the numerical solution smoothly onto the asymptotic expansion -see figure 18. While with the proper choice of a 2 the solution can be visually seen to join onto the asymptotic expansion smoothly, it is nice to have quantitative confirmation of this. In the right plot of figure 17 we show a residual that measures how closely the numerical solution matches the asymptotic expansion in a region where they overlap. The residual shown was calculated according to Residual = r fail 0.9r fail L V numeric (r) − V 1/r (r) r 2 dr , (C. 2) where again r fail is the point at which the numerical solution breaks down. In performing this calculation, terms up to O(r −3 ) where included in the asymptotic expansion. The plot shows that the error blows up a 2 L 2 = 0.1181855186708097 is approached from the left, while it goes to zero when approached from the right. This confirms that the numerical solution is indeed becoming arbitrarily close to the asymptotic solution, and the asymptotic solution can be used to continue the solution to infinity. On the contrary, there are some regions in the parameter space (for example, when the mass is negative in D = 4) for which we argued that the solutions do not exist due to a bad asymptotic behaviour. In those cases we are not able to match the numerical solution with the asymptotic expansion, and this confirms that those solutions do not exist.
Several strategies may be used in order to improve the precision of the numerical methods. For example, instead of working with the function V (r) one may work with f (r) = L 2 V (r)/r 2 , which should approach the constant f ∞ at infinity. Also, more terms can be included in the expansion (C.1), so that one does not need to choose a very small (we recall that the full expansion depends only on a 2 ). Let us close by mentioning that the numerical problem is considerably more stiff in D = 6 than in D = 4. In the latter case we do not require to increase substantially the WorkingPrecision and the parameter can be chosen as small as 10 −3 L. The numerical integration in D = 6 is less stable and requires a larger value of and higher values of WorkingPrecision.
Open Access. This article is distributed under the terms of the Creative Commons Attribution License (CC-BY 4.0), which permits any use, distribution and reproduction in any medium, provided the original author(s) and source are credited.