Vacuum decay and bubble nucleation in $f(R)$ gravity

In this work we study vacuum decay and bubble nucleation in models of $f(R)$ higher curvature gravity. Building upon the analysis of Coleman-De Luccia (CDL), we present the formalism to calculate the Euclidean action and the bounce solution for a general $f(R)$ gravity in the thin wall approximation. We calculate the size of the nucleated bubble and the decay exponent for the Starobinsky model and its higher power extensions. We have shown that in the Starobinsky model with a typical potential the nucleated bubble has a larger size in comparison to the CDL bubble and with a lower tunneling rate. However, for higher power extension of the Starobinsky model the size of the bubble and the tunneling exponent can be larger or smaller than the CDL bubble depending on the model parameters. As a counterintuitive example, we have shown that a bubble with a larger size than the CDL bubble but with a higher nucleation rate can be formed in $f(R)$ gravity.


Introduction
Tunneling is a quantum mechanical phenomenon which does not occur in classical physics. A minimum of the potential energy of the classical system is stable. However, in a quantum system where is not zero, uncertainty in the momentum causes the state of the system to become unstable. This is true for all minima of the potential energy, except the lowest one which remains stable. While the latter is called the "true" vacuum the others are called "false" vacua.
Vacuum tunneling plays important roles in various physical theories. For example, it is important for the Standard Model of particle physics since in the absence of new physics in higher energy scales, the vacuum of the Higgs boson can become unstable [1]. A similar but different issue is the degenerate vacua of non-abelian gauge theories [2]. In addition, it is believed that the landscape of string theory after compactification to four dimension posses many vacua. It is an ongoing debate how our primordial universe can emerge from tunneling in this vast landscape of vacua. Our main interests here are focused on tunneling in early universe. Since the energy scale of early universe after big bang is high new physics with multiple vacua can play important roles. Therefore, it is a well motivated question to study the effects of vacuum decay in early universe. In particular, tunneling and vacuum bubble nucleation played crucial roles in the development of old inflation scenarios [3,4].
The primary question regarding vacuum decay is the probability of its occurrence which is normally addressed using semi-classical approximations. The decay rate for a scalar field theory in flat spacetime was first studied by Coleman et. al. [5,6]. It was then generalized to include the effects of gravity by Coleman and De Luccia (CDL) [7], namely the "CDL" instanton (see also in [8]). It was shown that gravity can have significant effects on vacuum decay and bubble nucleation, specially for large bubbles, i.e. bubbles with sizes comparable to the horizon of the de Sitter background.
One natural question then is to study vacuum decay and bubble nucleation in modified theories of gravity. In particular, the Starobinsky's model of inflation based on higher order curvature gravity [9] is one of the most successful scenario of inflation in early universe which is well consistent with the cosmological observations [10,11]. In addition, models of modified gravity are studied as candidates to explain the mysterious natures of dark matter or dark energy, for a review see [12,13]. Among different classes of modified gravity, here we consider f (R) theories in which the action is given by a general well-posed function of the Ricci scalar where M P = 1/ √ 8πG is the reduced Planck mass. f (R) theories are simple but at the same time general enough to capture features of higher order curvature terms. To give some examples, Starobinsky's model of inflation is given by f (R) = R+αR 2 with α > 0 [9]. Theories with with f (R) = R − α/R n can cause late-time acceleration, though with drawbacks [13].
Related to our work, the effects of Gauss-Bonnet term on tunneling is considered in [14,15]. Also [16] studied the effect of non-minimal coupling to gravity. More recently, tunneling in Jordan-Brans-Dickie theory was studied in [17]. The paper is organized as follows. In section 2 we briefly review the formalism of vacuum decay in flat spacetime. In section 3 we present our analysis of vacuum decay in f (R) theory which is used in section 4 for the thin wall limit. Section 5 is devoted to f (R) = R + αR n as an example of our setup. In section 6 we reformulate the equivalent analysis in the Einstein frame and finally we conclude in section 7.

Vacuum decay in flat spacetime
Here we briefly review the vacuum decay in flat spacetime. For an extensive review see [18] and [19]. Consider the canonical scalar field where the potential U (φ) has two unequal minima as shown in Fig. 1. We denote the field values at the false and true vacuum respectively by φ 1 and φ 2 and the corresponding potential values by U 1 and U 2 . Consider the situation where the field value is initially at φ 1 in all space. As mentioned above, there is a nonzero probability at every point that the field jumps to the turning point quantum mechanically. This is allowed as tunneling respects energy conservation. The conventional approach to compute the probability is to find the wave function of the system by solving the time-independent Schrodinger equation. The probability amplitude of tunneling is proportional to the ratio of the wave functions at the turning point and the false vacuum. For a field theory, the configuration space is infinite dimensional and we must solve for a wave functional. The wave functional can be found via WKB approximation under the potential barrier. However, the WKB equation is not easily solved for a multidimensional system even at zeroth order. The idea due to [20] is that the wave functional is maximized on a path in the configuration space, namely the most probable escape path (MPEP), and exponentially suppressed in its neighbourhood. As a result, the multidimensional tunneling problem is reduced to a one dimensional case which is easily solved.
It can be shown that there is a parameterization τ of the MPEP curve in the configuration space that is a solution of the Euclidean equations of motion with appropriate boundary conditions [21,22]. This corresponds to replacing t → −iτ in the equation of motion, obtaining The boundary condition is that at τ → −∞ and τ = 0, φ is at the initial and final configurations respectively. The tunneling probability amplitude can be shown to be related to the Euclidean action S E of this solution. Since we have spacetime translation symmetry, the probability per unit volume of spacetime is with the decay exponent B = S E (bounce) − S E (false vacuum). The Euclidean action is computed on a solution that passes over MPEP twice, the second time in reverse from the turning point to the false vacuum, hence the "bounce" solution. It must then be subtracted from the false vacuum solution, i.e. in the absence of the bounce, where φ = φ 1 everywhere in spacetime. The prefactor can be calculated using the path integral methods [6] but its effect is subdominant compared to the exponential factor. Coleman [5] assumed and then confirmed in [23] that the MPEP must be an O(4) symmetric solution of (3). As a result, we must have with ξ = √ τ 2 + x 2 being the Euclidean distance from the origin. If we intuitively interpret τ as the time, the solution to this equation can be described as a growing bubble in which its interior is in true vacuum while its exterior is filled with the false vacuum. The bubble attains its maximum radius at τ = 0 and then shrinks and disappears in τ → ∞.
If the difference between the false and true vacuum energy ≡ U 1 −U 2 is small compared to the barrier height, it can be shown that the thickness of the wall can be neglected with respect to its radius. This is the thin wall approximation which can simplify the analysis considerably. In the limit of thin wall approximation the energy of the wall can be parameterized by the surface tension σ which can be calculated analytically. In this limit, the size of the bubble at the time of nucleation is given by ρ 0 = 3σ/ . The tunneling exponent is then obtained to be After formation, the bubble wall expands, approaching the speed of light asymptotically.
This formalism was generalized by Coleman and De Luccia (CDL) [7] to include the effects of Einstein gravity. It was shown that gravity can play important roles in vacuum tunneling and bubble nucleation, specially for large bubbles with sizes comparable to the Hubble radius of the background de Sitter spacetime.

Euclidean action for f (R) gravity
In this section we present our analysis of vacuum decay in f (R) gravity in which the CDL solution can be obtained in the special case where f (R) = R.
The starting action is in which the matter Lagrangian is the same as in Eq.
(2) containing a scalar field with the false and true vacua respectively at φ 1 and φ 2 . The boundary term S bnd. is added to take care of the total derivative terms. In Einstein theory of general relativity (GR) this is done by adding a total divergence term, the Gibbons-Hawking-York term. There is not a unique similar way in f (R) theories [12]. However, in correspondence with scalar-tensor theories one can obtain [25] where h is the determinant of the induced metric on the boundary, K is trace of the extrinsic curvature and F (R) ≡ df dR . The equations of motion are obtained by varying the action with respect to metric, yielding [13] Σ µν = 8πG T µν , in which is the generalization of the Einstein tensor and is the energy momentum tensor. For the scalar field we have Another useful equation is the trace of (9) which reads where T is the trace of the energy momentum tensor.
The Klein-Gordon equation also reads As discussed above, we have to write the Euclidean action with positive definite signature metric. This is given by where If we define the Euclidean energy momentum tensor then Eqs. (9) and (13) where dΩ 2 3 and dΩ 2 are the line elements of S 3 and S 2 respectively. The curvature of the ξ = const. three sphere is given by 6 ρ(ξ) 2 so ρ(ξ) represents the radius of the sphere in curved spacetime. In comparison, in flat spacetime ρ(ξ) = ξ.
The field φ is only a function of ξ and the equation becomes where a prime means differentiation with respect to ξ. The angular part of the integration in the Euclidean action (15) is trivial and we obtain where the Ricci scalar is given by This can be simplified by eliminating the kinetic term with the help of equations of motion On the other hand, from (10) we have Using these equations, the action is simplified to If we do the integration by part on the last term in the braces we obtain The surface term can be shown to be canceled by the boundary term since the trace of the extrinsic curvature of the hypersurface of ξ = const. is K = 3ρ ρ . In addition, the second term in the integral can be rewritten using the expression of the Ricci scalar Eq. (21) and we finally obtain For the Einstein GR the above action reduces to which is the expression used originally in [7]. We need to compute the Euclidean action for the false vacuum and the bounce solution. In the remaining of this section we perform the calculations for the false vacuum solution and devote the next section to the analysis of bounce solution.
We are looking for the solution of the equations of motion when the field is sitting at its false vacuum where φ(ξ) = φ 1 . Since by assumption the symmetry group of the solution is O(4), we have 6(= 4×3 2 ) Killing vectors for different rotations at each point. Additionally, for the case when the field is constant over the spacetime we also have 4 more Killing vectors for spacetime translation invariance. So, we have 10 Killing vectors overall which is the maximum number of symmetry transformations for a four-dimensional spacetime, i.e. the spacetime is maximally symmetric.
It can be shown that for any maximally symmetric spacetime the Ricci scalar is a constant. If we look at the trace of the generalized Einstein field equation, Eq. (13), we get where T = −4U 1 .
Let us introduce Λ 2 1 ≡ 3M 2 P U 1 as the length scale associated with the potential where we have assumed that U 1 is positive. Equation (28) is an algebraic equation for the Ricci scalar. The solution to this equation, whatever it is, we call it R 1 ≡ 12 to account for the dimension of the curvature. Thus L 1 is the length scale that characterizes the curvature of the spacetime. It is important to note that for a general f (R) theory L 1 = Λ 1 . However, one can check that for the Einstein GR, they are equal.
Eq. (28) may not have any positive solution for some theories or even it may contain more than one solution. However, here we assume that we have one and only one positive solution for both false and true vacuum.
The function ρ(ξ) can be deduced using the ξξ component of the generalized Einstein equation, Eq. (9) 3 where ). Combining the above equation with Eq. (28) we obtain The most general solution to this equation is in which a constant of integration is eliminated by an appropriate shift of coordinate such that ρ(ξ = 0) = 0. This is the well known expression for de Sitter solution in GR. Note that in our f (R) setup L 1 , and not Λ 1 , appears in ρ(ξ). This is understandable since L 1 determines the curvature scale of the spacetime. Similar expressions are also valid if U 1 is zero or negative. Let us assume that Eq. (28) has Minkowski solution with R 1 = U 1 = 0 or anti de Sitter (AdS) solution with R 1 = − 12 . The corresponding solution then becomes, Note that the spacetime is closed with no boundary for dS and the integration in Eq. (26) is over [0, πL 1 ] which is bounded by the zeros of ρ(ξ). However, for the other cases, the spacetime is open with boundary at infinity and the integration runs over [0, ∞). Finally, we can obtain Minkowski and AdS solutions from dS by L 1 → ∞ and L 1 → iL 1 . Thus, in what follows, we assume the potential is positive everywhere and the tunneling is from dS to dS spacetimes.

Decay exponent in the thin wall limit
To compute the decay rate, we need the Euclidean action for the bounce solution. This is the nontrivial solution of the equations of motion Eqs. (13) and (14) with the boundary conditions which for the closed spacetime are where ξ m is the second zero of the function ρ(ξ). It is not easy to find such a solution for a general potential, nor is it to prove its existence. In analogy to the CDL analysis, we assume that an instanton solution exists in which the variation of φ occurs in a small region near the bubble wall. We also assume that everywhere inside the bubble is filled with true vacuum while outside the bubble the state is in false vacuum. This is illustrated in Fig. 2.
Quantitatively, in the thin wall limit we assumeρ ρ ∆ξ, where ∆ξ is the wall thickness andρ = ρ(ξ) is the bubble radius withξ being the coordinate position of the wall. As a result, we can find φ(ξ) inside the wall approximately. More specifically, Eq. (14) can be written as However, in the thin wall limit, the integral over the friction term above can be discarded and we approximately have For the later use, this yields dξ = dφ The function ρ(ξ) is approximately known everywhere where L 2 is defined (similar to L 1 ) as the solution of Eq. (28) for the true vacuum with T = −4U 2 . We have assumed ρ wall (ξ) ≈ρ since we neglect the variation of ρ inside the wall. The form of the solution is such that the conditions ρ(0) = ρ(ξ m ) = 0 are trivially satisfied. By the continuity of the metric we must havē thus from the three constantsρ,ξ and ξ m only one is independent (neglecting the wall thickness). Here, we assume for the arguments The meaning of these conditions become more evident when we consider embedding of a dS spacetime as a sphere inside a five dimensional Euclidean spacetime. As a result, the bounce solution is a composite of two spheres with different radii (see Fig. 2). Conditions (39) mean that the position of the wall is further from the hemisphere in the false vacuum region. This is reasonable since we expect that the size of the bubble at the moment of formation to be much smaller than the size of the dS horizon. In the classification of [18,24] this is the type A bounce solution. Note that only the type A bounce solution exists in the flat spacetime solution in the limit of weak gravity. Therefore, in this work, we consider only the type A bounce solution.
In GR, solving φ(ξ) fixes the Ricci scalar algebraically through the trace of the Einstein equations. However, here we need to solve a differential equation where we have used Eq. (35) to eliminate the the kinetic term. It may not be easy to solve this equation in general.
The decay exponent is the Euclidean action for the bounce solution minus the Euclidean action for the false vacuum solution. The Euclidean action is given in Eq. (26). We can split the integral in Eq. (26) into three regions, inside the wall, on the wall and outside the wall, both for the bounce and the false vacuum solution, where the vertical lines show the integration limits. Below we calculate each contributions separately.
Since outside the wall the bounce solution is at false vacuum, by a change of the variable of the integral, we can write Note that if the spacetime is open as in Minkowski or AdS cases this is zero since the integral is up to infinity. For dS case, from Eqs. (38) and (39), we can write so the integration domain is of O( ) hence we may neglect B out . For B in , from (26), we get where we have used Eq. (29) in the second line. Finally, the contribution from the wall is where the integration is over the wall i.e.ξ < ξ <ξ + ∆ξ. The first term under the first integral, using Eq. (36), can be written as which is the surface tension of the bubble in the absence of gravity. Note that in GR, this is also the only surviving term in (47). However, in a general f (R) theory gravity has nonzero contributions to B w even in the thin wall limit. To put Eq. (47) into a more useful form, let us use Eq. (13) to write Since everything is only a function of ξ, we have The integration of the first term can be done and gives F 1 2 which is zero since the Ricci scalar and any function of it are constant in either sides of the wall. The integral of the second term can also be neglected, by the same argument that we have neglected the friction term in (34). Finally, using Eq. (35), the trace of the energy-momentum tensor can also be written in terms of the potential, so at the end we obtain where we have introduced S 1 and S 2 to account for the new terms. For GR, the integral in the first bracket is 3σ and the second integral vanishes i.e. S 1 = 3σ − 2σ = σ and S 2 = 0 and everything is consistent. Combining the contribution from B w and B in from Eqs. (46) and (51), the final expression for the decay exponent can be written as In order to findρ, the size of the bubble at the time of nucleation, we need to extremize the decay exponent given in Eq. (52), where we have introduced lowercase s i defined via Upon twice squaring, and defining x ≡ρ 2 , the equation can be written as where In general, the solution of Eq. (55) is given by where I and II are some complicated functions of the coefficients. The decay exponent is then computed if we substitute this solution into Eq. (52). In general this procedure is too complicated to be insightful. In the next section we employ the above formalism to the important case of f (R) = R + αR n in which the above equations can be solved analytically. Before that, it is useful to check our results with those of CDL in GR.

CDL solution
For GR, C 3 = C 4 = C 5 = 0, so the solution becomes which, after a straightforward calculation, results in the expression for the radius of the bubble in the CDL tunneling from dS to dS where we have defined R s ≡ 2M 2 P σ . If the true vacuum is a Minkowski space i.e. U 2 = 0 and Λ 2 → ∞, then where ρ 0 ≡ 2Λ 2 1 /R s = 3σ/ is the bubble radius in the absence of gravity and γ ≡ 1 + ρ 0 2Λ 1 2 is the prefactor due to GR. Note that in the the current case where U 2 = 0, Λ 2 Having obtained the extremum value of the bubble radius withρ = ρ cdl , the decay exponent B is obtained to be where B 0 is the value of the tunneling exponent in flat spacetime given by Eq. (6). It should be noted that the type A bounce in GR is valid if the following condition is met which is translated to the condition that γ < 2 in the above formulas.
5 Example: f (R) = R + αR n In this section we employ the formalism developed above to theories with f (R) = R + αR n in which α is a constant of dimension M 2(1−n) . Models with α > 0, n > 1 can be used to obtain inflation. For example, the Starobinsky model corresponds to n = 2 with α ∼ (3×10 −5 M P ) −2 in order to be consistent with the observations. In the regime α < 0, 0 < n < 1, it can be used for dark energy models. In this section, however, we restrict ourselves to models with n > 1 as in inflationary models.
Our assumption is that the correction to GR from f (R) is subleading so the theory to leading order is given by the conventional GR. This is motivated from the fact that the corrections to GR becomes important in high energy limit so at lower energies one should recover the standard GR. In our model with f (R) = R + αR n this approximation corresponds to αR n−1 1 .
The first step is to solve the de Sitter solution of the theory, namely solving Eq. (28) which results in The exact solution to this equations is not easy. However, employing the approximation Eq.
(62) we can solve it perturbatively in α. The solution to leading order in α is Curiously we note that for n = 2 the solution is the same as GR, i.e. L = Λ. Inside the wall we have to solve Eq. (40) which is Here we have defined which is the value of R(ξ) in the GR limit when α = 0 so the terms in the bracket in Eq. (65) represent the correction in R(ξ) to leading order in α. In addition, in obtaining the above result we have approximated the d'Alembertian with the second derivative inside the wall in the limit of thin wall approximation. The next step is to compute S 1 and S 2 . For S 1 from Eq. (51) we have where the subscript 1, as always, means the value of the corresponding variable at the false vacuum value. Here we have used the fact that the Ricci scalar at either sides of the wall is constant so its derivative is zero. We can write this in more compact form defining the parameter z via so that Similarly, one can show that with the new parameter y defined as The parameters z and y depend on n and the potential U (φ). In particular, for the Starobinsky model with n = 2, we have y = 6 regardless of the potential while z must be computed or each potential. In order to get some insights on their values and their dependence on the model parameters, we calculate them for a specific potential later on. We can computeρ to first order in α from Eq. (53). Since the expression is too complicated in its general form, here we present our analytical results for the case of tunneling from a de Sitter vacuum to a Minkowski vacuum. In this case, we obtain where ρ cdl is the size of the bubble in CDL analysis given in Eq. (59). Correspondingly, the tunneling exponent is obtained to be in which B cdl is the decay exponent in CDL analysis given in Eq. (60). In particular, for Starobinsky theory with n = 2, one obtain simpler expressions and where we have used the fact that y = 6 in this case. We see clearly that the size of the bubble at the time of nucleation and the tunneling exponent B depend non-trivially on α (calculated to first order) and n. We may deduce some general conclusions from Eq. (72). Note that z and y are always positive. Also we have 1 ≤ γ ≤ 2 so the sign of the expression in the brackets may be deduced from the sign of nγ 2 − (6n − 4)γ + 6n − 4. For 1 < n ≤ 2, which includes Starobinsky's theory, it is always nonnegative so we haveρ > ρ cdl e.g. in these theories the radius of the bubble at the moment of formation is bigger than the CDL bubble. It is worth noting that the radius of the bubble in GR is smaller than in the flat spacetime. As a result, higher order curvature terms compensate the effects of conventional gravity.
For n > 2, similar deduction holds in the interval However, for greater γ it should be checked numerically whether or not the bubble radius is bigger than the CDL bubble. Unfortunately, similar general conclusions can not be stated for the tunneling exponent Eq. (73) since the relative sign of the last two terms is negative and it is not trivial to compare their magnitudes. However, for the case n = 2 from Eq. (75) one can show that in the interval 1 ≤ γ ≤ 2( √ 3 − 1) ≈ 1.4 , B > B cdl while for greater γ, nothing general can be deduced. It is interesting to note that one can tune the potential so that the bubbles have bigger sizes compared to CDL bubbles while at the same time the tunneling rate is higher. This is because the domain of the parameter space where B > B cdl does not match in general with the domain whereρ > ρ cdl . This is rather counterintuitive since for both cases of flat spacetime and GR, bigger bubbles have lower rate (greater decay exponent). However, we see that this is not true in general in f (R) theories.
So far we did not specify the form of the potential. However, in order to get some insights about the numerical values of the various physical parameters, similar to [31,32], we consider the following potential  occurs at φ 2 = v with the condition λv 4 > 3 . Parameters and U 1 have their understood meaning and λ is a dimensionless scaling. The advantage of this form compared to [32] is that interesting quantities like can be varied without affecting others. Fig. 3 shows the potential for different values of λ. It is evident that by increasing λ, the top of the barrier is increased while other important properties do not change.
For this potential, the parameters z and y can be computed for different values of and λ which are shown in Fig. 4 and 5. Other parameters are specified in the captions of the figures. In addition, we have set h ≡ v M P = 0.1. It can be seen that both parameters have negligible dependence on while they normally increase by increasing λ.
The dependence of the bubble radius and the decay exponent to are shown in Fig. 6    where the quantities (ρ/ρ cdl − 1) and B/B cdl − 1 are plotted with λ = 100 and h = 0.1. In order to have a better visualization, the vertical axes are in log scales. The apparent singularities for the decay exponent (right panels) are because the corresponding quantity vanishes at ≈ 1.26v 4 . The range of parameters are such that γ covers the whole interval 1 ≤ γ ≤ 2. We have set α = 1 in appropriate units of M P since α scales everything equally to first order.
The same quantities are shown in Fig. 7, now varying λ while keeping fixed at = 0.5v 4 . Now, we have the apparent singularities at λ ≈ 40 where B/B cdl − 1 crosses zero.
It is evident that the radius of the bubble at the moment of materialization is enhanced due to the modification of gravity in our specific model. Specially, increasing while holding λ fixed result in a larger bubble radius. However, the dependence of the bubble radius on λ is different for different values of n. For Starobinsky theory, n = 2, it is monotonically decreasing, while for n = 3 and 4 it has more non-trivial dependence. The behavior of the decay exponent is more complicated. For n = 2, it is bigger than the case of GR and monotonically increases by increasing . For n = 3 and 4, it becomes smaller than GR for 1.26v 4 where (B/B cdl − 1) becomes negative in Fig. 6. This is a manifestation of the fact that in these theories bubbles with bigger sizes and yet with a higher rate than in GR can be formed. Similarly, decay exponent as a function of λ, drops below that of GR for λ 40.

Einstein frame formulation
The formalism developed above were based on Jordan frame. However, it is well-known that the f (R) higher curvature gravity is equivalent to a theory containing the standard Einstein gravity coupled minimally to a scalar field. Intuitively speaking, the new scalar field degree of freedom plays the role of the non-trivial function F (R) = 1. With this discussion in mind, one may ask how our previous result in Jordan frame can be translated in Einstein frame. One expects that the physical results to remain invariant under the change of frames, it is only the interpretation which changes. In this section we redo our analysis of tunneling in Einstein frame and verify that indeed the physical results are the same as in Jordan frame.
The action in the Einstein frame with the metricĝ µν is given by where F (A(ψ)) = exp 2 3 and F is the derivative of f , the original f (R) theory. The factor 2 3 is just a consequence of the spacetime dimension. This is the conformal factor between the two frames i.e.ĝ µν = F g µν . This equation must be solved to find A(ψ) and then the potential is given by For GR, V (ψ) = 0 and ψ becomes just a constant and the action reduces to the Einstein-Hilbert action coupled minimally to matter. For f (R) = R + αR n it is where F (ψ) is given by (79). We are interested in the Euclidean action with scalar field as matter where we have dropped the hat over the metric for convenience. Note that there is a coupling between ψ and φ via the exponential function (79) which is universal for any f (R) theory.
Assuming O(4) symmetry and conformal transformation, the line element is now given by which is conformally related to (18). Correspondingly, equations of motion for the scalar fields are φ + 3ρ and The first one is the familiar equation of motion for the scalar field with no change as expected. The second one is nothing but (13) if we change the variable from F to ψ, validating the wellknown result that the new scalar field degree of freedom is related to the derivative of f (R). As before, we can eliminate the kinetic terms of φ and ψ from (82), and after doing the integration by part to eliminate the second derivative of ρ and F , we obtain Note that although in the Einstein frame we have two scalar fields in presence of Einstein gravity, the Euclidean action is different since the coefficient of ρ term is 3F rather than only 3. The reason is that the field ψ appears directly in the metric and this results in nontrivial effects. Using Eq. (80), we see that the action (86) is exactly the same as the action (26), only expressed in terms of a different variable. The rest of the analysis goes as in Jordan frame.

Summary and discussions
In this paper, we studied the vacuum decay for a scalar field in f (R) gravity. The higher curvature terms are expected to appear in the effective theories of gravity in high energy physics. In addition, many theories of high energy physics predict various effective potential with multiple vacua. Therefore, the question of tunneling and vacuum bubble nucleation in theories of higher curvature gravity is well motivated. Among the theories of higher curvature gravity, f (R) theories are the best candidates which are free of the pathologies such as the ghost instabilities associated with higher curvature theories. In addition, they can be a viable model of inflation such as the Starobinsky model.
Our formalism was general, valid for any theory of f (R) gravity. However, to present compact analytical formula, we have considered the special class of f (R) = R + αR n with α as a perturbative parameter. Our motivation is that the higher curvature terms should appear as the high energy corrections to the Einstein gravity so we expect αR n−1 1. We have calculated analytically the size of the nucleated bubble and the decay exponent B in the case of tunneling from a dS space to a Minkowski background. In particular, for the Starobinsky model with n = 2, we have seen that the nucleated bubble has a bigger size compared to bubble created in GR. In addition, numerical calculations show that the rate of its materialization is lower. This is obtained for a typical potential for a range of parameters. Although we have not confirmed this analytically, it is expected to be true for potentials similar to the example studied here. For other values of n, depending on the model parameters, the bubble is bigger or smaller than the CDL bubble and also the decay rate is different. An interesting phenomenon which we have observed is that a bubble with a size bigger than the CDL bubble but with a higher nucleation rate can be formed. This is a counterintuitive phenomenon which happens as a consequence of the higher curvature corrections. Intuitively one may associate the effects of higher curvature terms as changing the effective Newton constant. However, our results show that the effects of higher curvature terms are more non-trivial than simply rescaling the effective Newton constant.
Before we close this paper, there are a number of comments worth discussing.
• It is believed that a dS spacetime has a temperature related to its curvature scale L via .
Since in f (R) theory L is in general different than Λ, then the associated temperature is modified in comparison to GR. For instance, for f (R) = R + αR n it is given by to first order in α with T 0 = 1/2πΛ. In [33] the tunneling rate using thermal fluctuations in the static patch in the fixed background approximation has been calculated. Furthermore, a plausible interpretation of CDL tunneling specially for big bubbles has been presented. We believe that their thermal interpretation can be applied to our case with the modified temperature Eq. (88).
• So far we have concentrated only on the CDL instanton. One may wonder the implications of our analysis to the Hawking-Moss instanton [34]. This is important since when the potential is too flat then the CDL bounce does not exist. In this case the decay exponent can be computed easily, yielding For f (R) = R + αR n theory it yields in which L top and Λ top represent the scales (as defined in section 3) calculated at the top of the potential.
Curiously we note that the correction term vanishes for Starobinsky model with n = 2. The reason is that the Hawking-Moss solution may be viewed to have thermal origin. For n = 2, we have L = Λ (as we have observed before) so the associated de Sitter temperature does not change and the nucleation rate remains unchanged.
• A bounce solution is a viable solution if it has one and only one negative mode i.e. the second variation of the Euclidean action has only one negative eigenvalue [6]. This has been investigated for GR in [35,36,37,38,39] and also for the case of non minimal coupling to gravity and vacuum decay seeded by black holes (inhomogeneities) [40,41,42] in [31]. Similar calculations must be done for f (R) gravity to check the viability of our results.
• In this work we have studied the tunneling procedure working in the Euclidean spacetime and left out the problem of evolution of the bubble after nucleation in the Lorentzian spacetime. It is normally assumed that the Lorentzian evolution is given by the analytic continuation of the Euclidean solution. Although this procedure may not be straightforward in curved spacetime (see fro example [26]), in analogy to the flat spacetime, the line element after nucleation is taken to be ds 2 = −ρ 2 dτ 2 + dξ 2 + ρ 2 cosh τ 2 dΩ 2 , where we have replaced χ → iτ . Another approach to this problem is to solve directly the equations of motion for the nucleated bubble in the presence of gravity. In GR this is done via Israel junction conditions [27,28]. The junction conditions for f (R) theories of gravity are developed in [29,30]. One important difference is that the Ricci scalar and the trace of the extrinsic curvature must be continuous over the junction surface. This is in contradiction with the fact that inside and outside of the bubble are in different vacua. The evolution of the bubble after nucleation in f (R) theories is an open question which we leave for future studies.
• There is another formalism to derive equations of motion from the action of a gravitational theory, namely the Palatini formalism. In this approach one treats the metric and the connection coefficients as independent fields. For the Einstein-Hilbert action this makes no difference compared to the conventional approach and it turns out that the connection coefficients must be the Levi-Civita connection. However, for f (R) theories the equations of motion are modified. Specifically, it can be shown that in the Einstein frame the Euclidean action is which is the same as (82) replacing F with ψ everywhere. The important difference is that there is no kinetic term for the ψ field i.e. f (R) theories in Palatini formalism are Brans-Dicke theories with ω BD = − 3 2 [12] so that the field ψ is non-dynamical. The equation of motion for ψ, assuming SO(4) symmetry, is This will change the solutions inside the wall of the bubble in the thin wall limit and will affect the final result.