Lorentzian Robin Universe

In this paper, we delve into the gravitational path integral of Gauss-Bonnet gravity in four spacetime dimensions, in the mini-superspace approximation. Our primary focus lies in investigating the transition amplitude between distinct boundary configurations. Of particular interest is the case of Robin boundary conditions, known to lead to a stable Universe in Einstein-Hilbert gravity, alongside Neumann boundary conditions. To ensure a consistent variational problem, we supplement the bulk action with suitable surface terms. This study leads us to compute the necessary surface terms required for Gauss-Bonnet gravity with the Robin boundary condition, which wasn't known earlier. Thereafter, we perform an exact computation of the transition amplitude. Through $\hbar\to0$ analysis, we discover that the Gauss-Bonnet gravity inherently favors the initial configuration, aligning with the Hartle-Hawking no-boundary proposal. Remarkably, as the Universe expands, it undergoes a transition from the Euclidean (imaginary time) to the Lorentzian signature (real time). To further reinforce our findings, we employ a saddle point analysis utilizing the Picard-Lefschetz methods. The saddle point analysis allows us to find the initial configurations which lead to Hartle-Hawking no-boundary Universe that agrees with the exact computations. Our study concludes that for positive Gauss-Bonnet coupling, initial configurations corresponding to the Hartle-Hawking no-boundary Universe gives dominant contribution in the gravitational path-integral.


I. INTRODUCTION
Transition amplitudes give us information on, whether a particular transition is allowed or not.A change of end points (also known as boundary) will affect the transition amplitude.Methods of quantum field theory (QFT) developed over the decades allow us to compute such transition amplitudes systematically, which can be matched with relevant experiments for verification.Infact, several experimental verifications of Standard Model predictions further strengthen our reliability on QFT and methods used to formulate it.path integral is one way to reliably study QFT and compute transition amplitudes systematically.
However, to properly define path integral in QFT one needs to address the issue of infinities that arise while computing transition amplitudes.These are taken care of by proper regularization and renormalization 1 .Besides these, one still has to specify an integration contour as the standard Lorentzian path integrals are not absolutely convergent and are highly oscillatory.In standard flat spacetime QFT (non-gravitational) this is achieved by Wick rotation: a process of going from real-time to imaginary time.This transforms the original flat spacetime Lorentzian path integral into a convergent Euclidean path integral.
These successes are hard to replicate when gravity is involved, where besides dealing with issues of divergences, renormalizability and gauge-fixing, one also has to address the issue of the contour of integration.For the gravitational system, this need not be standard Wick-rotation!Moreover, Euclidean gravitational path integral suffers from the conformal factor problem (the path integral over the conformal factor is unbounded from below [1]), motivating to study the gravitational path integral directly in the Lorentzian signature.A gravitational path integral with metric as the degree of freedom, on a manifold with boundaries can be written as where g µν is the metric and S[g µν ] is the corresponding action.The path integral is defined on manifold M with boundary ∂M.This abstract mathematical expression represents summation over all possible metrics with specified boundary conditions, where each geometry comes with a corresponding 'weight-factor'.
The gravitational action that we will focus on in this paper is given by following (see [2,3] for earlier works investigating the path integral of such gravitational theories) This is the Gauss-Bonnet gravity action where G is the Newton's gravitational constant, Λ is the cosmological constant term, α is the Gauss-Bonnet (GB) coupling and D is spacetime dimensionality.The mass dimensions of various couplings are: The action in eq. ( 2) falls in the class of lovelock gravity theories [4][5][6], and is a sub-class of higher-derivative gravity.In this, the dynamical evolution equation of the metric field remains second order in time always.Interestingly, the GB term also arises in the low-energy effective action of the heterotic string theory [7][8][9].It should also be emphasized that for the first time, the GB-coupling α has received observational constraints [10].These constraints arise from the analysis of the gravitational wave (GW) data of the event GW150914 which also offered the first observational confirmation of the area theorem [11].
Our interest in this paper is to study the path integral given in the eq.( 1) for the gravitational action specified by eq. ( 2), and to analyze the effects of the boundary conditions on the transition amplitude (see [12-14, 23, 24] for the role played by boundary terms).We start by considering a spatially homogeneous and isotropic metric in D spacetime dimensions.It is the FLRW metric in arbitrary spacetime dimension with dimensionality D. In polar coordinates {t p , r, θ, • • • } it is given by It has two unknown time-dependent functions: lapse N p (t p ) and scale-factor a(t p ), k = (0, ±1) is the curvature, and dΩ D−2 is the metric for the unit sphere in D − 2 spatial dimensions.This is the mini-superspace approximation of the metric.In this approximation, though, we don't have gravitational waves however, we do still retain diffeomorphism invariance of the time co-ordinate t p and the dynamical scale-factor a(t p ).This simple setting is enough to highlight the importance of boundary conditions and/or Gauss-Bonnet gravitational terms that may be becoming relevant at the boundary.The Feynman path integral with reduced degree of freedom becomes where beside the integration over the scale-factor a(t p ), lapse N p (t p ) and the Fermionic ghost C, we also have an integration over their corresponding conjugate momenta given by p, π, and P respectively.( ′ ) here denotes derivative with respect to t p , while the time t p co-ordinate is chosen to range from 0 ≤ t p ≤ 1. bd i and bd f refers to the field configuration at initial (t p = 0) and final (t p = 1) boundaries respectively.The Hamiltonian constraint H consists of two parts H = H grav [a, p] + H gh [N, π, C, P ] , (5) where H grav refers to the Hamiltonian for the gravitational action and the Batalin-Fradkin-Vilkovisky (BFV) [15] ghost Hamiltonian is denoted by H gh 2 .Although the degrees of freedom are quite reduced in mini-superspace approximation, the theory still retains time reparametrization invariance which need to be gauge-fixed (we choose proper-time gauge N ′ p = 0).For more elaborate discussion on BFV quantization process and ghost see [16][17][18].Most of the path integral in eq. ( 4) in the mini-superspace approximation can be performed exactly (the path integral over ghosts (C and P ) and the conjugate momenta (π and p)) leaving behind the following path integral The path integral Da(t p ) e iS[a,Np]/ℏ gives the transition amplitude for the Universe to evolve from one boundary configuration to another in the proper time N p , while the lapse integration implies that one needs to consider paths of every proper duration 0 < N p < ∞.This choice leads to causal evolution from one boundary configuration bd i to another bd f [19].The purpose of this paper is to study the path integral in eq. ( 6) for the gravitational action given in eq. ( 2) focusing mainly on the case of Robin boundary condition (RBC) imposed at the initial boundary.Imposing RBC at the initial time is like specifying a combination of field and its corresponding conjugate momenta.This translates into a oneparameter family of all allowed possible boundary conditions for a given specific value of the combination of field and its conjugate momenta.Motivation to study such configuration stems from the fact that imposition of Dirichlet BC at initial time leads to unsuppressed behavior of the gravitational fluctuations in the no-boundary proposal of the Universe [20][21][22].This forces one to investigate the situation in the case of Neumann and Robin BC at the initial boundary, and also tells that Dirichlet BC is perhaps not the right boundary condition for gravity [23,24].
Neumann boundary condition (NBC) has been studied in [2,3,[28][29][30] using Picard-Lefschetz methods, where it was seen that gravitational fluctuations are well-behaved in the no-boundary proposal of the Universe.Exact computations are done for the NBC case for the Gauss-Bonnet gravity, further showed some surprises mentioned in [3].For Einstein-Hilbert gravity, Robin boundary condition has been studied using the Picard-Lefschetz methods [28,31].These perturbative studies showed that fluctuations are well-behaved for the noboundary proposal of the Universe (see [32] for review on no-boundary proposal).
In this paper, we take a fresh look at the implications of imposing RBC at the initial time and study the path integral for the Gauss-Bonnet gravity.We study the Gauss-Bonnet gravity path integral using both perturbative and non-perturbative methods to gain a proper understanding of the effects of RBCs on the behavior of the transition amplitude.In the process, we also construct the suitable surface term for the Robin case needed for the Gauss-Bonnet gravity to have a consistent variational problem.We use time-slicing method of evaluating the gravitational path integral to compute the transition amplitude exactly.This is then compared with the results obtained via Picard-Lefschetz (PL) methods to gain a better understanding of the role played by saddle point geometries which could be complex.
PL methodology generalizes the notion of 'Wick-rotation' by adapting it to tackle highly oscillatory integral in a systematic manner 3 .By providing a framework that takes into account contributions from all the possibly relevant complex saddle points of the path integral, this framework uniquely determines contours of integration along which the oscillatory integrals (like the ones appearing in eq. ( 6)) becomes well-behaved.Such contours termed Lefschetz thimbles constitute the generalized Wick-rotated contour.In the context of Lorentzian quantum cosmology, these methods have been extensively used to analyze the nature of path integral [20][21][22][25][26][27] and the role played by boundaries [2,3,28,29,31]   PL-methods being inherently based on saddle-point approximation, fall short of analyzing situations dealing with degenerate case i.e. when saddles become degenerate.In such cases, one has to go beyond saddle-point approximation, as the approximation breaks down in such limiting cases.Non-perturbative and exact results whenever available become useful in investigating these such situations.This was particularly noticed in the study of gravitational path integral for the NBC case [3], where it was seen that our Universe undergoes a transition from an Euclidean to a Lorentzian phase as the size of the Universe increases.However, saddle-point analysis broke down at the turning point, a clear understanding of which came from the exact results.In the RBC case study, it is expected that a similar situation could arise.It is therefore best to approach the investigations of the gravitational path integral in the RBC case via both perturbative and non-perturbative manner.
The outline of the paper is following: section I gives introduction and motivation to the problem being studied.Section II studies path integral of a particle in potential in one-dimension with various boundary conditions.This not only helps us in understanding path integrals with non-trivial boundary conditions but also give rise to relations between 3 Some studies involving Wick-rotation in curved spacetime have been initiated in [33][34][35][36].However, more work needs to be done on this. 4 Earlier work using complex analysis methods to analyze Euclidean gravitational path integral was done in [37,38].Exploration of boundary effects in the Euclidean quantum cosmology was done in tunnelling proposal [39][40][41] and no-boundary proposal [37,38,42].Furthermore, the choice of a contour of integration via usage of complex analysis methods was also done in [43][44][45].
them which become useful later in the paper.Section III discusses the mini-superspace approximation and write the gravitational action in the mini-superspace approximation.Section IV studies the variational problem and computes the surface terms that are needed to have a consistent variational problem.Section V studies the transition amplitude for the NBC and RBC case, and compute the exact expression for it by exploiting the results of section II.Section VI studies the classical limit (ℏ → 0) of the exact transition amplitude which shows the boundary configuration giving a dominant contribution in the path integral.Section VII does the Picard-Lefschetz analysis of the contour integral.We conclude the paper in section VIII along with an outlook.

II. ONE-DIMENSIONAL QUANTUM MECHANICS
Before we proceed further with our study of gravitational path integral, we do a quick review of the path integral of one particle system to which our gravitational path integral reduces in the mini-superspace approximation.
We look at the path integral of one-particle system whose initial and final boundary conditions are specified.As the dynamical equation of motion involves two-derivatives of time, it needs only two boundary conditions.In the following, we will consider the cases where the final is always Dirichlet boundary condition and initial could be either Dirichlet or Neumann or Robin boundary condition.We will study a free-particle system and particle living in a linear potential which is of relevance to our studies in quantum cosmology.
Our purpose is to compute the following path integral where S tot [q] is the action for the one particle system m is a t-independent parameter, 'dot' denotes t-derivative, V (q) is the potential, S bd is the surface term added to have a consistent variational problem 5 .The 'bar' over G in eq. ( 7) is added to prevent confusion later on with gravitational transition amplitude.The path integral in eq. ( 8) can be computed by time-slicing method (first principles).The only subtlety arises at the end points around which the computation has to be done with care.To achieve this, we approach the problem in the Hamiltonian language as it is easier to incorporate the boundary conditions.The classical Hamiltonian for the above system is given by where p = m q.The quantum transition amplitude from one state to another in the Hamiltonian picture (which is equivalent to computing the above path integral) is given by, Bd f , t = 1 Bd i , t = 0 = Bd f e −iH(p,q)/ℏ Bd i .
Here we have lifted p and q to operators p and q respectively, which then obey the commutation relation q, p = iℏ .
It should be mentioned that the effect of surface terms S bd gets taken care in the Hamiltonian style of doing the computation from one boundary state to another.The Hamiltonian follows from the Lagrangian by Legendre transform.
A. Dirichlet boundary condition at t = 0 Let's first consider the case of a free particle where V (q) = 0.The quantum Hamiltonian is given by One is interested in computing a transition amplitude from one state to another.However, these states are given by specifying the initial and final position: q(t = 1) = q f and q(t = 0) = q i .The path integral in eq. ( 7) becomes This can be evaluated from first principles by the time-slicing method, which discretizes the problem.Breaking the time-length into M segments gives the time duration of each segment to be ϵ = 1/M .This means that q(t) :⇒ q 0 = q i , q 1 , q 2 , • • • , q k , •, q M −1 , q M = q f , with the end points q 0 = q i and q M = q f are fixed.This implies that we can write the transition amplitude given in eq. ( 13) in the following discretized form where Consider the following matrix element ⟨q k e −iH free ϵ/ℏ q k−1 ⟩.In this, one can plug completenessrelation in momentum space and perform the momentum integral where we performed the Gaussian integral in p k to obtain the last expression.Eventually, we are left with the following expression for Ḡfree The series of Gaussian q-integrals can be carried out one-by-one starting from q 1 till we reach integration over q M −1 which is the final integration that one has to perform.This will lead to the following expression for Ḡfree DBC (q f , t = 1; q i , t = 0) where we have used M ϵ = 1, when taking the limit M → ∞ and ϵ → 0. In the case interaction term V (q) is present, then eq. ( 17) gets modified into the following This is nothing but a discretized version of the path integral given in eq. ( 7) where qk = (q k − q k−1 )/ϵ.However, for a generic potential V (q) this can't be computed exactly.In the case for linear potential V (q) = λq (which will also appear in quantum cosmology, as we will later see), one can easily compute the above-discretized version of the path integral exactly.The end result after the q-integrations is the following, Notice the quadratic structure of the λ-dependent terms.The part coming from interaction vanishes for two values of λ: λ = 0 (free theory) and λ = −12m(q f + q i ) (interacting).

B. Neumann boundary condition at t = 0
The path integral for the Neumann boundary condition is easier to deal with when handled in the Hamiltonian framework.The quantity we are interested in computing is where at the initial time, we fix the momentum.We note that the initial momentum state can be written as the following dq 0 e ip i q 0 /ℏ q 0 ⟩ .
This allows one to express ḠNBC as a Fourier transform of another transition amplitude dq 0 e ip i q 0 /ℏ ⟨q f e −iH/ℏ q 0 ⟩ ḠDBC (q f ,t=1;q 0 ,t=0) This shows that ḠNBC (q f , t = 1; p i , t = 0) is related to ḠDBC (q f , t = 1; q 0 , t = 0) by a Fourier transform thereby implying that one can be obtained from another by transformation (if ḠNBC (q f , t = 1; p i , t = 0) is known then ḠDBC (q f , t = 1; q 0 , t = 0) can be obtained from it by an inverse Fourier transform).In the present situation as we have already worked out the transition amplitude ḠDBC (q f , t = 1; q 0 , t = 0), the above Fourier transform can be used to compute the expression for the Neumann transition amplitude.This is given by Ḡfree NBC (q f , t = 1; p i , t = 0) = e i(p i q f −p 2 i /2m)/ℏ .
Similarly, one can exploit the above technique to work out the path integral for a particle in the linear-potential with Neumann boundary condition at the initial time.Here if we plug eq. ( 19) on the RHS of eq. ( 23) and writing q i → q 0 , then on integration with respect to q 0 we get the ḠNBC (q f , p i ).
The interaction dependent term is quadratic in λ and vanishes for two values of λ: λ = 0 and λ = 3(p i − 2mq f ).

C. Robin boundary condition at t = 0
Now we are interested in computing the path integral for the case when Robin boundary condition (RBC) is specified at the initial time t = 0.An RBC at t = 0 means that a given linear combination of position and momentum is fixed at that time.
The boundary value problem we are interested in is the following Bd f : q(t = 1) = q f and Bd i : where q i ≡ q(t = 0) and p i ≡ p(t = 0).The path integral we are interested in computing is An interesting way to evaluate this is to do a canonical transformation and define a new momentum variable P and position variable Q as Under this transformation, the commutation relation is preserved, i.e., This also means that the Jacobian of transformation is unity.However, in terms of new variables, the Hamiltonian acquires a new form H(p, q) = H( P −β Q, Q) ≡ H β ( P , Q) (where subscript 'β' symbolizes β-dependence of the Hamiltonian).In terms of new variables, our transition amplitude acquires the following form In terms of new variables the Robin boundary condition problem transforms into an NBC problem.One can write the initial state as follows This means that for the Robin boundary condition the path integral becomes the following ḠRBC (q f , t = 1; The object of interest here is ⟨Q f e −iH β ( P , Q)/ℏ Q 0 ⟩ which is the DBC path integral for the transformed Hamiltonian H β (P, Q).It should be pointed out that this process of evaluating RBC path integrals remain valid even for interacting theories, but to gain clarity on the methods, we first consider free theory.The free Hamiltonian H free β is given by and in transformed variables the boundary conditions become Therefore, our original interests in computing path integral with RBC shifts to computing path integral in new variables but with NBC.However, the Hamiltonian in new variables is different which means that one needs to be careful while computing path integral via time-slicing method.We are interested in computing something like eq. ( 14) with the Hamiltonian given in eq. ( 33).Let us consider the following matrix This means that The term in the exponent which can be written as a total derivative (β/2) 1 0 dt dQ 2 /dt.This eventually gives a simplified expression for the Robin-boundary condition transition amplitude.
This can be easily computed after plugging the value of Ḡfree DBC from eq. ( 18) which gives which in limit β → 0 reproduces the expression for the Neumann boundary condition result given. in eq. ( 24).
In the case when interactions are present, the modified Hamiltonian gets amended by potential term V (Q) Then following the same steps as above one arrives at an analogous expression for which is given by Once again the term βQ can be written as a total derivative (β/2) 1 0 dt dQ 2 /dt.Then the expression for the ḠRBC with the potential is given by It should be highlighted that this last line is valid for any arbitrary potential V (Q).For the case of linear potential one can plug the expression for ḠDBC [Q M , t = 1; Q 0 , t = 0] from the eq.( 19) in the above to obtain ḠRBC (q f , t = 1; Notice that in the limit β → 0 this RBC expression goes to the NBC expression given in eq. ( 25).The full Robin path integral for interacting theory is a product of Robin path integral for free theory and a coupling dependent piece.The expression given in eq. ( 41) can be converted into a relation between ḠRBC (Q f , t = 1; P i , t = 0) and ḠNBC (q f , t = 1; p i , t = 0) by exploiting the inverse Fourier transform.After a bit of algebra and computation of Gaussian integrals, it is seen that where P i = p i + β q i .This is an exact relation connecting the two quantities.In the following, we will use this relation to compute the exact expression of gravitational path integral in mini-superspace approximation for Robin BC and investigate it in detail to study the implications of Robin BC on the evolution of the Universe.

III. MINI-SUPERSPACE ACTION
We start by considering the FLRW metric given in eq. ( 3) which is conformally-flat and hence has a vanishing Weyl-tensor (C µνρσ = 0).For a generic d-dimensional FLRW metric the non-zero entries of the Riemann tensor are [46][47][48] Here g ij is the spatial part of the FLRW metric while ( ′ ) denotes the derivative with respect to t p .The non-zero components of Ricci-tensor are while the Ricci-scalar for FLRW is given by Weyl-flat metrics also have a wonderful property that allows one to express the Riemann tensor in terms of Ricci-tensor and Ricci scalar as follows This furthermore helps one to express the square of Riemann-tensor as a linear combination of Ricci-tensor-square and Ricci-scalar-square.This is given by This identity is often used to simplify the expression for the Gauss-Bonnet gravity action for the case of Weyl-flat metrics.
On plugging the FLRW metric of eq. ( 3) in the gravitational action stated in eq. ( 2), we get an action for scale-factor a(t p ) and lapse V D−1 here refers to the volume of the D − 1 dimensional spatial hyper-surface which is given by, In four spacetime dimensions (D = 4) the terms proportional to α (Gauss-Bonnet coupling parameter) becomes a total time-derivative.In D = 4 the mini-superspace gravitational action then becomes the following At this point, one can introduce rescaling of lapse N p and scale factor to bring out a more appealing form of the mini-superspace action by doing the following transformation This co-ordinate transformation re-expresses our original FLRW metric in eq. ( 3) in a different form.
Notice also that in the new co-ordinate frame, the new 'time' is denoted by t.Under this transformation, the original gravitational action in D = 4 given in eq. ( 52) acquires a following simple form Here (˙) represent time t derivative.Note that the GB-term is a total time derivative.Integration by parts of the terms in the action re-expresses the action in a more recognizable form along with some surface terms.
The surface terms consist of two parts: one coming from EH-part of gravitational action while the other is from GB sector.The bulk term is like an action of a one-particle system in a linear potential.From now onwards we will work with the convention that V 3 = 8πG, and in the rest of the paper we will study the path integral of this action.

IV. ACTION VARIATION AND BOUNDARY TERMS
We start by constructing the variational problem by considering the variation of the action in eq. ( 55) with respect to q(t).This will not only yield terms that will eventually lead to the equation of motion but also generate a collection of boundary terms.
In the rest of the paper, we will work with the ADM gauge Ṅ = 0, which implies setting N (t) = N c (constant).To set up the variational problem consistently and avoid missing of any terms we write q(t) = q(t) + ϵδq(t) , where q(t) is some background q(t) and δq(t) is the fluctuation around this.The parameter ϵ is used to keep track of the order of fluctuation terms.On plugging this in the action given in eq. ( 55) and on expanding it to first order in ϵ, we notice that the terms proportional to ϵ are given by following We also notice that in the above expression, there are two total time-derivative pieces which contribute at the boundaries.These boundary pieces need to be canceled appropriately by supplementing the action with suitable surface terms in order to have a consistent boundary value problem for a particular choice of boundary condition.The terms proportional to δq on the other hand will lead to the equation of motion if one demands that the full expression multiplying it vanishes.This gives the dynamical equation obeyed by q(t) This process of construing surface action is based on analyzing the behavior of fluctuations around the classical trajectory (on-shell geometries), and may sometime miss to properly capture the features of topological effects as by definition topological terms don't play a role in the equation of motion.This will become more clear as we proceed further in the paper.
The second-order ODE given in eq(59) is easy to solve and its general solution is given by The constants c 1,2 are determined by requiring the solution to satisfy the boundary conditions.The total-derivative terms in eq. ( 58) lead to a collection of boundary terms where qi = q(t = 0) , qf = q(t = 1) , qi = q(t = 0) , qf = q(t = 1) .
Note that the boundary terms are function of q(i,f) and q(i,f) which are on-shell quantities.These are different from q (i,f ) and q(i,f) which also includes boundary values of the off-shell trajectories.The bulk action given in eq. ( 56) is used to determine the momentum conjugate to the field q(t) In terms of conjugate momentum the boundary terms can be expressed as follows For the consistent variational problem these terms either have to be appropriately canceled by supplementing the original action given in eq. ( 55) with suitable surface terms or they vanish identically for the choice of boundary conditions.
In the following, we will always impose Dirichlet boundary condition at t = 1, while at t = 0 we will either have Neumann or Robin boundary conditions (as imposing Dirichlet at t = 0 leads to unsuppressed perturbations, so it won't be addressed [20,21]).However, as the path integral in two cases (NBC and RBC) are related by eq. ( 43) therefore one can first compute the path integral in NBC case and use the relation in eq. ( 43) to compute the results for the Robin BC.

A. Neumann Boundary condition (NBC) at t = 0
Imposing Neumann boundary condition [23,31] at t = 0 and a Dirichlet boundary condition at t = 1 gives a consistent variational problem.This means that the initial momentum π i and final position q f of all the trajectories are fixed.This will mean Such a boundary condition imposition also has an additional merit that they lead to a wellbehaved path integral where the perturbations are suppressed [2,[28][29][30].For this case, the boundary terms given in eq. ( 64) arising during the variation of action will reduce to the following These residual boundary terms which don't vanish after imposing boundary condition need to be canceled by a suitable addition of surface terms so that the variational problem is consistent.It is seen that if one adds the following terms at the boundary then its variation (as suggested in eq.(57) will precisely cancel the terms given in eq. ( 66).The first term in eq. ( 67) is a Gibbon-Hawking-York (GHY)-term and the second one is a Chern-Simon like term at the final boundary.The constants c 1,2 that appear in the solution to the equation of motion eq. ( 60) can now be determined for the choice of boundary conditions.This will imply where 'bar' over q is added as it is the solution to equation of motion.In this setting at t = 0, we notice that the on-shell value of q i is given by qi For off-shell trajectories, q i can be anything as that is not fixed at t = 0.The surface terms obtained in eq. ( 67) when added to the action in eq. ( 56) leads to the full action of the system.This is given by Furthermore, on substituting the solution to the equation of motion eq. ( 68) and using eq.(69) in the above, we will obtain the on-shell action which is given by This is also the action for the lapse N c .Compare this action with the action computed when Dirichlet boundary condition is imposed at t = 0 [2, 20-22, 28-30] and notice the lack of singularity at N c = 0.The lack of singularity at N c = 0 can be physically understood as in the NBC path integral we sum over all possible transitions from q i (fixed π i ) to fixed q f .This will also include the transition from q f to q f which can occur instantaneously i.e. with N c = 0. Thereby implying that there is nothing singular happening at N c = 0.

B. Robin Boundary condition (RBC) at t = 0
We now consider the situation when we have Robin boundary condition at t = 0 and Dirichlet boundary condition imposed at t = 1.An example of imposing Robin BC is fixing the Hubble parameter at a given time in the cosmological evolution.This boundary value problem poses a consistent variational problem.In this case, we fix a linear combination of the initial scale factor q i and the corresponding initial conjugate momentum π i , while the final scale factor q f is fixed at t = 1.This means For these boundary conditions, it has also been noted that perturbations are suppressed and in path integral their effects are bounded [2,[28][29][30].For the boundary conditions mentioned in eq. ( 72), the boundary terms generated during the variation of action given in eq. ( 64) reduces to the following For the consistent variational problem, these terms need to be removed by supplementing the original gravity action with suitable surface terms.The additional suitable surface terms when varied will precisely cancel the above terms leading to a consistent variational problem.It is seen that if one adds the following terms at the boundary Here the first two terms are the same as for the NBC case, while the rest of the terms arise due to the imposition of RBC at t = 0.In the limit of β → 0 these terms will vanish and we get back the surface terms for the NBC which gives us a check of consistency.The variation of term βq 2 f /2 is although zero (due to imposition of DBC at t = 1), however adding this term helps us later while computing path integral via comparison with results given in section II C. The surface term βq 2 i /2 in eq. ( 74) comes due to the RBC at t = 0 from the Einstein-Hilbert part of the gravitational action and was known from before [31], The term proportional to αβ is new and corresponds to surface terms that need to be added for the Gauss-Bonnet gravity part.
The surface terms obtained in eq. ( 74) when added to the action in eq. ( 56) leads to the full action of the system.This is given by Comparing this action with the action for the NBC case given in eq.(70).In the limit β → 0 as P i → π i .The constants c 1,2 that appear in the solution to the equation of motion given in eq. ( 60) can now be determined for the case of RBC.This will imply where 'bar' over q indicates the solution of the equation of motion.Setting t = 0 in this gives the on-shell value of the q i which is given by qi Off-shell q i could be anything as it is not fixed by the boundary condition imposed in RBC.
The on-shell action can be computed by substituting the solution of the equation of motion and qi into the action for RBC given in eq. ( 75).This is given by This is the action for the lapse N c in the case of RBC at t = 0 and DBC at t = 1.In the limit β → 0 this reduces to the NBC on-shell action given in eq.(71).

V. TRANSITION AMPLITUDE
We now move forward to compute the transition amplitude from one 3-geometry to another (for both NBC and RBC case).The quantity that is of relevance here in the minisuperspace approximation is as follows (see [20,43] for the Euclidean gravitational path integral in mini-superspace approximation) where bd i and bd f are the initial and final boundary configurations respectively, and S tot refers to the mini-superspace action, which for the NBC case is given by eq. ( 70).We will first compute the expression for G NBC [bd f , bd i ] using eq.( 80) and ( 25).We will then use the relation given in eq. ( 43) to compute the expression for G RBC [bd f , bd i ].
A. NBC at t = 0 To compute the expression for the transition amplitude in the case of Neumann boundary condition we make use of the results given in.eq. ( 25).However, to make this connection we first compare the one-particle action written in eq. ( 8) with the mini-superspace action for the NBC case given in eq. ( 70).This shows that the following substitution needs to be made At this point, we are interested in computing This second line can be computed by exploiting the results given in eq. ( 25) and making use of the substitutions mentioned in eq. ( 81).Note that the presence of term q i π i with an integration over the initial q i leads to Fourier transform as expected in the NBC case.This gives the following expression for the Neumann boundary condition where S on−shell tot [q, N c ] is given in eq. ( 71) and ḠNBC [q f , t = 1; π i , t = 0] is.given by Notice that eq. ( 83) this differs from the result obtained in [3] by a numerical prefactor.This difference in normalization is attributed to different styles of doing computation of path integral involving zeta-functions.The process of computing path integral via time-slicing is more reliable and correctly fixes the normalization.
Note S on−shell tot [q, N c ] doesn't become singular at N c = 0 and computation of the path integral over q(t) doesn't give rise to additional functional dependence on N c which can be singular.This allows us to extend the limits of N c -integration all the way upto −∞.This means that the transition amplitude is given by where S on−shell tot [q, N c ] is given in eq.(71).To deal with the lapse integration in the NBC case we first make a change of variables.This is done by shifting the lapse N c by a constant This change of variable correspondingly implies that the action for the lapse S on−shell tot [q, N c ] becomes the following An interesting outcome of this is that after the change of variables the π i (initial momentum) dependence only appears in the constant term.After the change of variables, the transition amplitude is given by where The transition amplitude for the Neumann BC at t = 0 and Dirichlet BC at the t = 1 is a product of two parts: Ψ 1 (π i )Ψ 2 (q f ).Ψ 1 (π i ) is entirely dependent on initial momentum π i and other Ψ 2 (q f ) is function of q f is related to the final size of the Universe.The dependence on two boundaries gets separated, a factorization also noticed in [29] (and also in [3]) where the authors studied the Wheeler-DeWitt (WdW) equation in mini-superspace approximation of Einstein-Hilbert gravity.

B. RBC at t = 0
We now come to the task of computing eq. ( 80) for the case of RBC at t = 0, where S tot [q, N c ] for the RBC case is given in eq. ( 75).This mini-superspace RBC action can be compared with the quantum mechanical RBC problem discussed in subsection II C.This comparison leads to the substitution mentioned in eq.(81).At this point, we are interested in computing where It is noted that the expression for ḠRBC [q f , t = 1; P i , t = 0] is exactly referring to the RBC path integral in terms of DBC path integral with the substitution given in eq. ( 81).Except that N has been added to get rid of e iβq 2 f /2ℏ due to requirement of classicality, where in a WKB sense, the amplitude is expected to become more classical as the Universe expands [20,28].However, using eq.( 43) one can express ḠRBC [q f , t = 1; P i , t = 0] in terms of ḠNBC via the integral transform.This means we have The full path integral in the RBC case also involves integration over lapse N c with limits (0, ∞).This means As the N c integrand is not singular at N c = 0 so one can extend the integration limit all the way up to −∞.The N c -integration thereafter becomes similar to the integral studied in the case of NBC except the lack of Gauss-Bonnet term dependent on α.This means one can write where Ψ 1 and Ψ 2 are given in eq. ( 89) and ( 90) respectively.Here again, as the integrand doesn't have N c = 0 singularity, the integral can be extended all to way to −∞ and including an extra 1/2 factor which is absorbed in the definition of Ψ 1 .The leftover integral is the integral over p.This integral can be cast into a more familiar form by re-definition of p as Such a transformation allows us to rewrite the p integral as an Airy integral.This is given by (97) where Putting all the pieces together give It is noticed that for the RBC case too the transition amplitude gets factorized in two parts: one dependent on the final boundary and one dependent on the initial boundary.In the next sub-section, we will compute these integrals in terms of Airy functions.It is expected that the RBC transition amplitude will be a product of two Airy functions along with exponential prefactors.
In order to see the β → 0 limit clearly it is more cleaner to write the β-dependent terms in eq. ( 93) as follows In this way in the limit β → 0 gives a δ-function δ(P i − p).This limit is more harder to see from the end result of transition amplitude as one has to work with asymptotic forms of the Airy-functions.

C. Airy function
In this sub-section, we will compute the integrals Ψ 2 (q f ) and Φ(P i , β) given in eq. ( 90) and ( 98) respectively.It should be noted that these functions can be identified with the Airy-integrals.These integrals are sensitive to the contour of integration.In the case of Airy-integrals, the regions of convergence are within the following phase angles θ ≡ arg( N ): 0 ≤ θ ≤ π/3 (region 1), 2π/3 ≤ θ ≤ π (region 0), and 4π/3 ≤ θ ≤ 5π/3 (region 2).One can define the following contours: C 0 the contour running from region 0 to region 1, C 1 the contour running from region 1 to region 2, and C 2 the contour running from region 2 to region 0. By making use of the above contours of integration one can define the following integrals Ψ 2 (q f ) can be computed as discussed in detail in the paper [3].It is given by The computation for Φ(P i , β) can be done in an analogous manner.It is given by, Putting these expressions together gives the exact expression for the transition amplitude for the NBC and RBC case.These are given by, G RBC [bd 0 , bd 1 ] = 3πℏ 2iβ exp 4iα Ai 3 This is an exact result for the case when NBC and RBC are imposed at initial time t = 0 respectively.Notice also the correction coming from the Gauss-Bonnet sector of the gravitational action which appears as an exponential prefactor.It is also crucial to notice that this GB correction doesn't appear in the Airy functions which only depend on the boundary conditions.

VI. ℏ → 0 LIMIT
In this section, we will study the ℏ → 0 limit of the exact wave-function computed in the previous section.It is relevant to look at this as in the ℏ → 0 limit the exact computation should reproduce known features and should agree with the results coming from saddle-point approximation which will be discussed in more detail in section VII.The ℏ → 0 limit also highlights the configuration space which gives the dominant contribution in the path integral.This ultimately translates into preferable values for the initial boundary parameters.One can then do the saddle point analysis for these values of initial parameters which gives the dominant contribution in the path integral.
We will start by analyzing the nature of the contribution coming from the Gauss-Bonnet sector and the constraints that come from it in the limit ℏ → 0.  110) vs y.The two saddle points are depicted by green and red squares.For α > 0: green one corresponds to P i = −3i √ k (stable) while the red one corresponds to P i = 3i √ k (unstable).Interestingly, P = −3i is also the value for which geometry become regularised at t = 0 and perturbations are well-behaved.
The first thing we will focus on is the exponential prefactor that includes the effects coming from the Gauss-Bonnet sector.This prefactor is the same in both NBC and RBC cases (P i → π i in limit β → 0).This prefactor is the following exp 4iα where A(P i ) is given by In the limit ℏ → 0 the dominant contribution comes from configuration which extremizes A(P i ), which corresponds to the In the limit β → 0 this corresponds to π i = ±3i √ k!Note that π i = −3i √ k also corresponds to the required value of the initial momentum for regular geometries with well-behaved perturbations in the case of Hartle-Hawking (HH) wave-function for the NBC [3] (3i √ k corresponds to Vilenkin where perturbations are unstable).
It should however be emphasised that ℏ → 0 limit imposes classicalization when the quantum system starts behaving classically.This same end state where the system behaves like a classical system can also be achieved when ℏ is not small but the Gauss-Bonnet coupling is large.In this situation when the Gauss-Bonnet coupling α is large, then the dominant contribution comes from the configuration whose P i lies around P i = ±3i √ k.This situation however is entirely quantum in nature as ℏ is not small while the system acquires the end state of 'classicality' in the limit of large Gauss-Bonnet coupling.This can happen in the very early stages of the Universe.
Let us consider the rotated P i = iy, which is like considering a Euclideanised version action A. Then we have the Euclidean action A E (y) given by On extremization of Euclidean action A E (y) it is seen that for α > 0 the point y = −3 √ k corresponds to a stable saddle point and leads to an exponential with a positive argument (Hartle-Hawking) while y = 3 √ k correspond to unstable saddle point leading to exponential with a negative argument (Vilenkin Tunneling).Coincidently, as was also mentioned in [3,28] √ k also correspond to initial condition for which geometry becomes regularised at t = 0 and perturbations are well-behaved.

B. Choice of β
The next question that one needs to address is that for the given P i = −3i √ k (as required by stability), what are the allowed possibilities for the β in the case of RBC transition amplitude.This issue doesn't arise in the case of NBC as β = 0 for the NBC transition amplitude.To analyze the constraints (allowed values) that get imposed on β from the ℏ → 0 limit, it is sufficient to study the nature of RBC wave-function at t = 0 for the fixed allowed value of P i = −3i √ k.We start by considering the ℏ → 0 limit of Φ(P i , β) which for P i = −3i √ k can be written as For imaginary β the argument of Airy-function is always positive for Λ > 0. This observation allows us to obtain easily the ℏ → 0 limit by exploiting the asymptotic structure of the Airyfunctions with positive arguments where and the ∼ sign implies that we are ignoring the numerical prefactor arising from the asymptotic form of Airy's function, which is not relevant for the following discussion.From the asymptotic structure, one could see that in limit ℏ → 0 the dominant contribution comes from those configurations for which B ′ 1 (β) = 0.This means we have This means that the dominant contribution comes from configuration when β lies around β dom = −iΛ/2 √ k (we don't consider the β → ∞ case).We next consider the exponential prefactor independent of α for P i = −3i √ k.This means that we have to investigate exp where Once again in the ℏ → 0 limit the dominant contribution comes from those configurations for which β lies around the extrema given by B ′ 2 (β) = 0.This given β dom = −iΛ/2 √ k same as before.This eventually give rise to two cases: iβ < iβ dom and iβ > iβ dom (note β and β dom are imaginary).To study these cases carefully let us write This gives Now we consider cases when x ≤ 1 and x > 1 respectively.
In this case, we will have This can be combined with the exponential prefactor coming from the Gauss-Bonnet sector of gravity (α-dependent part) which for P i = −3i √ k is given by exp(8k 3/2 α/ℏ).Combining the two exponential prefactors gives the following This is precisely the Hartle-Hawking state with positive weighting.It is worth emphasizing that the exponential prefactor is independent of β.The positive exponent says that lower Λ is more favorable, i.e., low values of potential are preferred.The analysis also includes the case β = 0 i.e., Neumann BC.For a special value of α = −3/4Λ, this HH-factor becomes unity.

x > 1
In this case, we will have This exponential factor has to be multiplied with the contribution coming from the Gauss-Bonnet sector e 8k 3/2 α/ℏ to get the full overall exponential prefactor.In the large x → ∞ limit we get inverse Hartle-Hawking.However, for other values of x > 1, the nature of exponential depend on the behavior of the functional form of x.Let us call this function We notice that f (1) = 1 and f (∞) = −1.It should be mentioned that f ′ (x) = −6(x − 1) 2 /x 4 < 0 for all values x except x = 1 and x = ∞.This means f (x) is a monotonically decreasing function of x within the range between 1 and −1.At some point the function f (x) changes sign and becomes negative for x > x 0 while it remains positive for 1 < x < x 0 (although it is not equal to 1, thereby implying that it doesn't have the actual Hartle-Hawking factor 6k 3/2 /ℏ Λ as in eq. ( 120)).

VII. SADDLE-POINT APPROXIMATION
The analysis of ℏ → 0 has shown us that certain initial configurations are favorable and give the dominant contribution in the path integral.The ℏ → 0 limit of the Gauss-Bonnet contribution shows that the dominant contribution comes from P i = ±3i √ k.Of which only the P i = −3i √ k correspond to a stable configuration.Coincidentally, this is also the same required value of P i which leads to stable well-behaved fluctuations in the Hartle-Hawking no-boundary proposal of the Universe [2,3,[28][29][30][31].In a sense, the Gauss-Bonnet presence favors the stable Hartle-Hawking no-boundary Universe.
In this section, we will do the saddle point analysis of the gravitational path integral in the mini-superspace approximation.We will do this for the case of Robin boundary condition as the case of Neumann boundary condition has been already investigated in [3].However, the saddle-point study of the path integral will be done in a slightly different manner.The gravitational path integral in the mini-superspace approximation in the RBC case is given in eq. ( 80) where the path integral over the q(t) is given in eq.(91).We use the methods developed in analyzing the quantum-mechanical RBC problem to convert eq.(91) into eq.( 93), which relates the RBC path integral of q(t) to the NBC path integral of the same gravitational theory.Doing lapse N c integration over this gives the full transition amplitude, relating the RBC transition amplitude with the NBC one as given in eq. ( 94).This allows us to write the full RBC transition amplitude as a product of two integrals after we make use of eq.(95).One is integral over p and other integral over N , where N -integral is given by eq.(90).Our strategy in doing the saddle-point analysis is to analyze each of these integrals (integral over N c and integral over p) separately.
We first take a look at the integral over N c given in eq. ( 95).If one shifts N c as in eq.(86) then this integral becomes an integral over N and an exponential prefactor dependent on p.The N -integral depends on q f written as Ψ 2 (q f ).The saddle point analysis of this has already been done in [3] and won't be repeated here.The results from the saddle analysis revealed an interesting feature.For all values of q f there are two saddle points, however, for q f < 3k/Λ only one saddle point which lies on the imaginary axis is relevant, while for q f > 3k/Λ there are two relevant saddle points lying on the real axis.This means that for q f < 3k/Λ Universe is Euclidean as can be seen from the exponential form of the transition amplitude indicating imaginary 'time', while for q f > 3k/Λ Universe is Lorentzian as is noticed from the oscillatory nature of transition amplitude indicating real time [3].
Our task then shifts to analyzing the integral over p given in eq.( 97). where and h(p) is the corresponding Morse-function while the H(p) corresponds to real part of the B-function.By a shift of variable as stated in eq. ( 96), the p-integral can be cast into an Airy integral along with an exponential prefactor as discussed in section V B. Here, we will study integral in eq. ( 124) using Picard-Lefschetz methods (see [2,20,[49][50][51][52] for review on Picard-Lefschtez and analytic continuation).
The saddle-points of p can be obtained by computing the expression dB(p)/dp.The saddle points are computed from the equation This is a quadratic equation in p resulting in two saddle points.The discriminant ∆ of the above quadratic equation is given by For stable configuration referring to Hartle-Hawking no-boundary Universe (P i = −3i √ k) and β given by eq. ( 117), the discriminant ∆ becomes the following For all values of x > 0, the discriminant ∆ < 0 implying that both the saddle points are complex, which for stable Hartle-Hawking no-boundary Universe is seen to be both imaginary.These are given by This shows that while saddle point p1 remains fixed at the same position for all x, the saddle point p2 moves from negative imaginary axis to positive imaginary axis as x increases.It becomes zero for x = 2.The 'on-shell' value of B(p) (B(p) computed at the saddle points) is given by It should be highlighted that for stable Hartle-Hawking no-boundary Universe referring to P i = −3i √ k the 'action' B(p) given in eq. ( 125) is complex.The morse function h(p) at the saddle points for all values of x is given by Higher derivatives of B(p) with respect to p at the saddle point is given by, One can expand the function B(p) around the saddle point pσ where σ = {1, 2}.This gives where δ p = p − pσ .The series stops at cubic order as the highest power of p in B(p) is three.
The second variation at the saddle-point can be written as B ′′ (p σ ) = r σ e iρσ , where r σ and ρ σ depends on boundary conditions.Near the saddle point the change in iH will go like where we write δ p = v σ e iθσ and θ σ is the direction of flow lines at the corresponding saddle point.Given that the imaginary part H remains constant along the flow lines, so this means where k ∈ Z.
For the steepest descent and ascent flow lines, their corresponding θ des/aes σ is such that the phase for δH correspond to e i(π/2+2θσ+ρσ) = ∓1.This implies These angles can be computed numerically in our case.Under the saddle point approximation, the contour integral given in eq. ( 124) can be computed using Picard-Lefschetz methods (see [2] for details of PL methodology).This gives e iθσ e iB(pσ)/ℏ .
where J σ refers to steepest descent line and n σ is the intersection number which will take values (0, ±1) accounting for orientation of contour over each thimble.In the following we will be using this expression to compute the contour integral for various values of x.
A. 0 < x < 1 For 0 < x < 1, both the saddle points p1 and p2 lie on the negative imaginary axis with |p 2 | < |p 1 |, thereby implying that p2 lie below p1 on the negative imaginary axis.2. The gravitational path integral with the Gauss-Bonnet contribution favours P i = −3i √ k as the stable configuration which also happens to correspond to Hartle-Hawking no-boundary Universe.For this value of P i the action for B given in eq. ( 125) is complex.This figure arises in association with the Picard-Lefschetz analysis of the contour integral given in eq. ( 124).The red lines correspond to steepest descent lines (thimbles J σ ), while the thin black lines are steepest ascent lines and denoted by K σ .Here we choose parameter values: k = 1, Λ = 3, and x = 1/2.For this, the saddle point p1 corresponds to blue-square, while saddle point p2 correspond to bluecircle.Only the saddle point p1 (blue square) is relevant.The steepest ascent lines emanating from it intersects the original integration contour (−∞, +∞) which is shown by thick-black line.The Morse-function h is positive for both saddle points: h(p 1,2 ) > 0. The light-green region is allowed region with h < h(p σ ) for all values of σ.The light-pink region (forbidden region) has h > h(p σ ) for all σ.The intermediate region is depicted in yellow.The boundary of these region is depicted in brown lines.Along these line we have h = h(p σ ).

■ ■
For both the saddle points the morse function is positive: h(p 1,2 ) > 0. In figure 2 we plot the various flow-line, saddle points, forbidden/allowed regions.From the graph we notice that only the steepest ascent line from p1 intersects the original integration contour, thereby making it relevant.The thimbles passing through this saddle constitute the deformed contour of integration.It should also be specified that for x < 1, the second-derivative of B(p) computed at saddle point p1 is proportional to +i, which in saddle point approximation gives appropriate Gaussian weight allowing to do Gaussian integration.
The Picard-Lefschetz theory then gives the following in the saddle point approximation as This regime reproduces the same exponential factor as in Hartle-Hawking no-boundary Universe.

B. x > 1
For x > 1 both saddles still lie on the imaginary axis.However, only saddle point p1 remains fixed on the negative imaginary axis for all values of x, while the saddle point p2 moves from negative imaginary axis to positive imaginary axis.It crosses the origin at x = 2 in the case of boundary condition corresponding to Hartle-Hawking no-boundary Universe.

■ ■
3.Here we choose parameter values: k = 1, Λ = 3, and x = 4 and x = 10 respectively.For this, the saddle point p2 corresponds to blue-circle, while saddle point p1 correspond to blue-square.Only the saddle point p2 (blue square) is relevant.The Morse-function h for p1 is always positive, while h(p 2 ) goes from positive to negative as 1 < x < ∞.The crossover happens at x = x 0 given by eq. ( 123).
The Morse function corresponding to saddle point p1 is always positive, while the Morse function corresponding to p2 goes from positive to negative.The crossover take place at x = x 0 where x 0 is given in eq. ( 123).After the cross-over the overall exponential weight becomes negative, thereby implying an inverse Hartle-Hawking regime.For x > 1 it is seen that only the steepest ascent curves emanating from p2 intersect the original integration contour, implying that it is the only relevant saddle point.The thimbles passing through this saddle will constitute the deformed contour of integration.The Picard-Lefschetz analysis then gives the following in the saddle point approximation as This is the degenerate case.The discriminant ∆ = 0 for x = 1 and both saddle points coincide: p1 = p2 = −3i √ k.In this degenerate case both the on-shell action (B(p 1 ) = B(p 2 )) and the Morse-function become equal (h(p 1 ) = h(p 2 ) = 6k 3/2 /Λ).In this situation the saddle-point approximation breakdown as the second derivative B ′′ (p σ ) = 0.One needs to go beyond the second order in the series expansion in eq. ( 133).In this series the third order term is non-zero.This situation is depicted in figure 4. In this case, the contour integration gives the following

VIII. CONCLUSIONS AND OUTLOOK
In this paper we consider the gravitational path integral of Gauss-Bonnet gravity and study it directly in the Lorentzian signature in four spacetime dimensions.Gauss-Bonnet sector of gravity being topological in nature in 4D doesn't contribute to the bulk dynamics of the field but has an active role to play at the boundaries or in situations where boundaries play an important role.One such situation is the path integral which is sensitive to boundary conditions.Past studies have investigated the effects of Gauss-Bonnet sector of gravity on the transition amplitude which is given by the gravitational path integral [2,3].These studies focussed on exploring the consequence of imposing Neumann boundary condition at the initial time.In this paper we study extensively the effects of imposing Robin boundary condition at the initial time and investigate the role played by the Gauss-Bonnet sector of gravity systematically in the Lorentzian gravitational path integral.We setup the platform for raising and addressing these issues in the reduced setup of the mini-superspace approximation.
We start by first considering path integrals of particle in one-dimension in a potential for various boundary condition.The path integral is evaluated for three different choices of boundary condition: Dirichlet boundary condition (DBC) at both initial and final time (fixing position of particle at two end points), Neumann BC (fixing conjugate momenta denoted in paper by p i ) and Robin BC (fixing linear combination of conjugate momenta and position at initial time: p i + βq i = P i , where β is some parameter) with Dirichlet at final time.The transition amplitudes are shown to be inter-related with each other by integral transforms.By exploiting these inter-relations via integral transforms one can compute path integral with NBC or RBC at initial time (and DBC at final time) from the path integral of DBC at initial time.These inter-relations can be further manipulated to express the RBC path integral as an integral transform of the NBC path integral.These results get later utilized in the paper for the computation of gravitational path integral in the mini-superspace.
We take a fresh look at the gravitational path integral in the mini-superspace approximation in the Lorentzian signature with Robin boundary condition at the initial boundary (that is fixing linear combination of conjugate momenta and field at the initial boundary).Past works have shown that NBC and RBC at the initial time leads to stable Universe [2,3,[28][29][30][31], while DBC at initial time leads to unsuppressed perturbations [20][21][22].These motivate us to study RBC gravitational path integral more carefully for the Gauss-Bonnet gravity (the case for Neumann BC was investigated in [3]).The transition amplitude from one 3-geometry to another is given by a path integral over q(t) and a contour integration over lapse N c .
The paper systematically studies the path integral in the mini-superspace approximation for the Gauss-Bonnet gravity.The gravitational action is varied and carefully analysed to setup a consistent variational problem with Robin BC at the initial boundary.This process leads to dynamical equation of motion and a set of surface terms that need to be supplemented to the action to make the gravitational system a consistent variational problem.In this process we construct the surface term needed for the Gauss-bonnet gravity with the Robin boundary condition.This surface action is proportional to β and smoothly reduces to the surface term for the case of Neumann BC problem [3] in the limit of β → 0. To our knowledge this wasn't known in literature earlier.The path integral is then studied for the total action which has been supplemented by these surface terms.
The path integral to be computed is given in eq.80, where S tot for the NBC and RBC case are given in eq. ( 70) and (75) respectively.To compute the path integral over q(t) we make use of the results derived in one-dimensional quantum mechanical problem discussed in section II.In the NBC case: the lapse N c integration and path integral over q(t) can be performed exactly as it was also shown in [3].The numerical prefactor has been correctly computed as we do the computation of path integral via first principles (method of 'timeslicing').In the RBC case: results of sub-section II C and algebraic manipulations allow us to reduce the problem into the product of two Airy-integrals (modulo exponential prefactors).via Picard-Lefschetz methods.
The saddle analysis of the p-integration shows that in the complex p-plane there are always two saddle points, both lying on the imaginary p axis.The saddle point p1 is independent of x (remains fixed), while the saddle point p2 varies with respect to x, and gets pushed to infinity (on the negative imaginary axis) as x → 0 (Neumann limit).For all values of x, only one of the two saddle point is relevant.For 0 < x < 1 the saddle point p1 relevant, while the irrelevant saddle point p2 lies below p1 on the negative imaginary axis in the complex p plane (in the limit x → 0 this saddle point is pushed to infinity and remains irrelevant).In this regime, we get the correct Hartle-Hawking exponential prefactor e 6k 3/2 /ℏ Λ.For x > 1, the saddle p2 becomes relevant.x = 1 is the degenerate situation when p1 = p2 (relevant) giving e 6k 3/2 /ℏ Λ (Hartle-Hawking) beside the numerical pre-factors.For 1 < x < 2 the saddle point p2 still lies on the negative imaginary axis.However, for x > 2 it crosses over and lies on the positive imaginary axis.For 1 < x < x 0 , the argument of the exponential prefactor decreases in magnitude and becomes negative for x > x 0 .The range of x for which only the saddle point p1 is relevant and produces the exact exponential prefactor e 6k 3/2 /ℏ Λ of Hartle-Hawking no-boundary Universe is 0 < x < 1.
Our investigations show the important non-trivial role played by the Gauss-Bonnet sector of the gravitational action in favoring the boundary configurations which lead to the Hartle-Hawking no-boundary Universe.Gauss-Bonnet term although is topological in nature in four spacetime dimensions and is expected to not play any role in the dynamical evolution.But, our analysis clearly shows its contribution in the path integral when not ignored naturally picks and favors boundary configurations which correspond to Hartle-Hawking no-boundary Universe.Although, it still remains to be seen whether the behavior of fluctuations can alter this.

FIG. 1 .
FIG. 1. Plot of Euclidean actionA E given in eq.(110) vs y.The two saddle points are depicted by green and red squares.For α > 0: green one corresponds to P i = −3i √ k (stable) while the red one corresponds to P i = 3i √ k (unstable).Interestingly, P = −3i is also the value for which geometry become regularised at t = 0 and perturbations are well-behaved.