Relativistic stars in bigravity theory

Assuming static and spherically symmetric spacetimes in the ghost-free bigravity theory, we find a relativistic star solution, which is very close to that in general relativity. The coupling constants are classified into two classes: Class [I] and Class [II]. Although the Vainshtein screening mechanism is found in the weak gravitational field for both classes, we find that there is no regular solution beyond the critical value of the compactness in Class [I]. This implies that the maximum mass of a neutron star in Class [I] becomes much smaller than that in GR. On the other hand, for the solution in Class [II], the Vainshtein screening mechanism works well even in a relativistic star and the result in GR is recovered.

in which GR is recovered by the mechanism similar to that in the ghost condensation [25] as well as by the Vainshtein mechanism. However, in these analysis, gravitational fields are assumed to be weak. It has not been cleared whether the Vainshtein screening mechanism holds even in the strong gravitational field (e.g., relativistic star and black hole).
The black hole geometry in bigravity has been concerned, which are classified into non-diagonal ansatz [26], and bi-diagonal ansatz [22,24,[27][28][29]. In the former type ansatz, there are only trivial solutions, which are the same as those in GR 1 . Additionally, the perturbation around the non-diagonal black hole is also identical to GR [32][33][34]. Hence, the massive graviton does not appear in the non-diagonal black hole. To find a non-trivial solution, if it exists, we should assume both metrics can be simultaneously diagonal in same coordinate system.
There exists some special case of the bi-diagonal ansatz such that two metrics are proportional, which we call a homothetic spacetime. The solutions are also given by those in GR. However, in this case, the massive graviton appears in the perturbation around the solutions. As a result, the homothetic Schwarzschild black hole becomes unstable against the radial perturbations if the graviton mass is sufficiently small [35][36][37]. The instability of this black hole implies that there would be a hairy black hole solution as well, and that the homothetic Schwarzschild black hole may transit to the hairy black hole. However, the paper [29] showed numerically that such a hairy black hole does not exist unless the coupling constants satisfy a special condition. One may wonder what we will find in the final stage of gravitational collapse of a compact relativistic star. One may also ask whether there exists a maximum mass of neutron star, beyond which no neutron star cannot exist.
The standard picture in GR is that a star collapses to a black hole when the mass exceeds the maximum value. However, in bigravity, although there exists a Newtonian star solution in the weak gravitational field, no stable black hole solution has been found for generic coupling constants. In order to investigate what happens when a star is compact and relativistic and then the gravitational interaction becomes very strong, we study a relativistic star in the bi-gravity theory. A little attention has so far been paid to a relativistic star in the bigravity. Hence, as a first step, we analyze a star solution with a relativistic effect, and discuss how such a relativistic star behaves in the limit of strong gravity.
In the text, we assume that only g-matter field exists and spacetime is asymptotically flat. We then classify the coupling constants into two classes: Class [I] and Class [II]. For Class [I], we find an example of breaking Vainshtein screening mechanism due to the relativistic effect. The static star solution is found when the pressure of the star is sufficiently small, while the star solution disappears when the pressure is larger than a critical value. Therefore, in Class [I], the maximum mass of the neutron star in bigravity is constrained stronger than one in GR. On the other hand, there is no critical value of the pressure for Class [II]. The result of GR is reproduced even in the strong gravitational field.
The paper is organized as follows. The Hassan-Rosen bigravity model is introduced in Sec. II. In Sec. III, we derive the basic equations in bi-diagonal ansatz of the static and spherically symmetric spacetime. Taking the limit of massless graviton, we discuss behaviours of the solutions deep inside the Vainshtein radius in Sec. IV. We find that the existence of a neutron star solution is restricted depending on the coupling constants. In Sec. V, we numerically solve the basic equations without taking the massless limit, and confirm that the previous solutions with massless limit approximation are valid if the Compton wave length of the graviton mass is sufficiently large compared to the typical radius of the star. We summarize our results and give some remarks in Sec. VI. In appendix A, we summarize the parameter constraint from the existence of a Newtonian star. In Appendix B, introducing a cosmological constant and f -matter field, we discuss solutions with asymptotically non-flat geometry. In Appendix C, we detail the case beyond the critical value of the pressure for Class [I], in which we find a singular behaviour.

II. HASSAN-ROSEN BIGRAVITY MODEL
We focus on the ghost-free bigravity theory proposed by Hassan and Rosen [3], whose action is given by where g µν and f µν are two dynamical metrics, and R(g) and R(f ) are their Ricci scalars. The parameters κ 2 g = 8πG and κ 2 f = 8πG are the corresponding gravitational constants, while κ is defined by κ 2 = κ 2 g +κ 2 f . We assume that the matter action S [m] is divided into two parts: g (g, ψ g ) + S [m] f (f, ψ f ) , (2.2) i.e., matter fields ψ g and ψ f are coupled only to the gmetric and to the f -metric, respectively. We call ψ g and ψ f twin matter fluids [38].
The ghost-free interaction term between the two metrics is given by where {b k } (k = 0 -4) are coupling constants and the 4×4 matrix while U k are the elementary symmetric polynomials of the eigenvalues of the matrix γ, defined explicitly in [11,12]. Taking the variation of the action with respect to g µν and f µν , we find two sets of the Einstein equations: where G µ ν and G µ ν are the Einstein tensors for g µν and f µν , respectively. The γ-"energy-momentum" tensors T [γ]µ ν and T [γ]µ ν are obtained by the variation of the interaction term with respect to g µν and f µν , respectively, taking the form [11,12] The matter energy-momentum tensors T [m]µ ν and T [m]µ ν are given by the variation of matter actions. They are assumed to be conserved individually as where (g) ∇ µ and (f ) ∇ µ are covariant derivatives with respect to g µν and f µν . From the contracted Bianchi identities for (2.5) and (2.6), the conservation of the γ-"energymomenta" is also guaranteed as These equations give non-trivial constraints on solutions, which are absent in GR.

III. STATIC AND SPHERICALLY SYMMETRIC SPACETIMES
To find a non-trivial static and spherically symmetric regular solution, we assume two metrics are bi-diagonal in same coordinate system. Thus, we consider the following metric forms: where the variables {N g , F g , r g , N f , F f , r f } are functions of a radial coordinate r, and a prime denotes the derivative with respect to r. The ansatz has two residual gauge freedoms: One is a rescaling of time coordinate (t →t = ct with c being a constant), and the other is redefinition of the radial coordinate (r →r(r)). The proportional constant factor K is introduced just for convenience. K is one of the real roots of the quartic equation When N g /N f = F g /F f = r g /r f = 1, g-and f -spacetimes are homothetic and the γ energy-momentum tensors turn to be just "effective" cosmological terms. In the text, we focus on asymptotically homothetic solutions, i.e., we assume the boundary condition Solutions with other asymptotic geometrical structure will be discussed in Appendix B.
We introduce new variable µ defined by with µ > −1, which determines the relation between two radial coordinates r g and r f . From the boundary condition, µ should approach zero at infinity.
Introducing new parameters as (3.9) the Einstein equations are reduced to We have two more Einstein equations, which are automatically satisfied since we have two Bianchi identities for g µν and f µν .
In the original Lagrangian, we have six unfixed coupling constants {κ f , b i }, where m is not independent because it is just a normalization factor of b i . In this paper, we use six different combinations of those constants; {m g , m f , Λ g , K, β 2 , β 3 }, in stead of {κ f , b i }, because the behaviours of the solutions within the Vainshtein radius are characterized by β 2 and β 3 as we will see later. The original coupling constants {κ f , b i } are found from {m g , m f , Λ g , K, β 2 , β 3 }.
The energy-momentum conservation laws of twin matters give where we assume that twin matters are perfect fluids. The energy-momentum conservation laws of the interaction terms, which are equivalent to the Bianchi identities, reduce to one constraint equation; Substituting the Einstein equations (3.12) and (3.14) into Eq. (3.17), we obtain one algebraic equation: (3.18) Now we have nine variables N g , N f , F g , F f , µ, ρ g , P g , ρ f and P f , and six ordinary differential equations (3.11)-(3.14), (3.15), (3.16) and one algebraic equation (3.18) with two equations of state P g = P g (ρ g ) and P f = P f (ρ f ). In order to solve those equations numerically, we first take the derivative of (3.18), and then find seven first-order ordinary differential equations: where X = {N g , N f , F g , F f , P g , P f }, and F X and J do not contain any derivatives. Here we have fixed the radial coordinate as r g = r by use of the gauge freedom. We solve these differential equations from the center of a star (r = 0). In order to guarantee that the above set up gives a correct solution of our system, we have to impose the constraint (3.18) on the variables at the center. Note that the proportional factor K is not necessary to be unity. Since K appears only in the form of K 2 ρ f and K 2 P f , however, unless f matter exists, the basic equations are free from the value of K. In what follows, we assume that there is no f -matter just for simplicity. The f -matter effect on the solution will be discussed in Appendix B 3.

IV. REGULAR COMPACT OBJECTS : MASSLESS LIMIT
Before we present our numerical solutions, we shall discuss some analytic features of a compact object. The radius of neutron star is about 10 6 cm, while the Vainshtein radius is given typically by 10 20 cm when the Compton wave length of the graviton mass is the cosmological scale (m −1 eff ∼ 10 28 cm). The magnitude of the interaction term, which is proportional to the graviton mass squared, is much smaller than the density of a neutron star. Hence, the interaction term seems not to affect the structure of a neutron star. If we ignore the interaction terms in the Einstein equations (2.5) and (2.6) (or Eqs. (3.11)-(3.14)), we just find two independent Einstein equations in GR. Then both spacetimes are given approximately by GR solutions, which we can solve easily. In bigravity theory, however, we have one additional non-trivial constraint equation (2.10) (or (3.18) for a static and spherically symmetric case) even in the massless limit. This constraint will restrict the existence of the solutions. In this section, we consider a compact object in this massless limit.
Note that, in this massless limit, the effective action to determine the the Stückelberg variable µ is given by where Λ 2 = m/κ, and g GR and f GR are solutions in GR which act as like external forces to the Stückelberg field 2 . This effective action is indeed the same as the non-compact nonlinear sigma model proposed by [39]. As we will see, the massless limit approximation is valid deep inside the Vainshtein radius. It implies that, inside the Vainshtain radius, the non-compact nonlinear sigma model with a curved metric is obtained as the effective theory for the Stückelberg field. We analyze two models: one is a simple toy model of a relativistic star, i.e., a uniform-density star, and the other is a more realistic polytropic star with an appropriate equation of state for a neutron star.
A. The boundary condition at "infinity" in the massless limit The boundary condition at spatial infinity, which is outside of the Vainshtein radius, is given by Eq. (3.5). Since the radius of a neutron star is much smaller than the Vainshtein radius, there exists the weak gravity region even inside of the Vainshtein radius. We then introduce an intermediate scale where R ⋆ and R V are the radius of a star and the Vainshtein radius, respectively. The space inside the Vainshtein radius can be divided into two regions: the region deep inside the Vainshtein radius (r < R I ) and the weak gravity region (R I < r < R V ), where the gravitational force is described by a linear gravitational potential.
From the analysis for the Vainshtein screening in the weak gravity system [17,23,24], we find that GR (or Newtonian) gravity is recovered in r < R V , while the homothetic solution is obtained outside the Vainshtein radius r ≫ R V . The function µ(r) changes from −1/ √ β 3 at small distance (r ≪ R V ) to 0 at large distance (r ≫ R V ). When gravity is weak, we find µ ≈ −1/ √ β 3 deep inside of the Vainshtein radius. Hence we expect that µ ≈ −1/ √ β 3 at r ≈ R I for a relativistic star. We then obtain the boundary condition for a relativistic star in the massless limit as as r → R I , which we can assume R I ≈ ∞ because R I ≫ R ⋆ . Note that in the massless limit, the Vainshtein radius turns to be infinite.

B. Uniform-density star
First, we consider a uniform-density star. Since the basic equations in the massless limit are just the Einstein equations, we can easily solve them. The g-metric of this g-star is given by the interior and exterior Schwarzschild solutions, while the f -metric is just a Minkowski spacetime: For the interior (r < R ⋆ ), while for the exterior (r > R ⋆ ), where R ⋆ and are the g-star radius and the gravitational mass, respectively. Although we can choose N g (0) (or N f (0)) any value by the rescaling of time coordinate, from the boundary condition N g /N f = 1 at infinity (R I ), we find the ratio as Only one variable µ has not been solved. When we find a regular solution of µ(r) for the constraint (3.18) in the whole coordinate region (0 ≤ r < ∞) with the boundary condition µ → −1/ √ β 3 as r → ∞, we can construct a relativistic star in the bigravity theory.
First we analyze the constraint (3.18) at the center r = 0 (r f = 0), which gives where µ 0 := µ(0). This is the quadratic equation of µ 0 , which does not guarantee the existence of a real root of µ 0 . In order to have a real root µ 0 , we have one additional constraint as We then classify the coupling constants β 2 and β 3 into three cases: In the case (1), the real root µ 0 exists only for the restricted range of P g (0)/ρ g , In fact, there are two critical values; w − and w + (w + > w − ), which are defined by , (4.13) and the real root exists either if On the other hand, for the case (2) and (3), the real root µ 0 always exists for any value of P g (0)/ρ g .
Furthermore, when we take into account the finiteness of the graviton mass, even if it is very small, we find an additional constraint on the coupling constants {β 2 , β 3 } from the existence of non-relativistic star with asymptotically homothetic spacetime [24] (see also Appendix A).
Since the case (2) is completely excluded, we find two classes of the coupling parameters, which provide a relativistic star with asymptotically homothetic spacetime, as follows: 3 where d 1 and d 2 are some complicated functions of β 3 , which are defined by (A10) and (A12) in Appendix A, respectively. Assuming β 3 > 1, which is necessary for the existence of asymptotically homothetic solution, we show the Even if a real µ 0 exists, we may not find a regular solution of µ(r) in the whole coordinate range (0 ≤ r < ∞) because the real root of (3.18) may disappear at some finite radius. In Figs while for Class [II], the parameters are chosen as and Note that there are two real roots for µ 0 . Then we find two branches of µ(r), which we call the branch A and the branch B. The branch A approaches a homothetic solution (µ → −1/ √ β 3 ) as r → ∞ in the massless limit, while the branch B (µ → 1/ √ β 3 ) does not become homothetic at infinity. For the Class [I] example (4.14), µ 0 exists only if figure). We find a regular solution for both branches if P g (0)/ρ g < 1/15. The branch A solutions provide relativistic stars with asymptotically homothetic spacetime, while the branch B solutions are not asymptotically flat.
For 1/15 < P g (0)/ρ g < 1/3, µ 0 does not exist. We find the solution µ(r) only for the region larger than some finite radius, and two branches A and B are connected. The topology of this spacetime is similar to a wormhole, but it has a curvature singularity at the throat (the turning point of µ(r)). For the large value of P g (0)/ρ g , the turning point appears outside of the "star", which means the "wormhole" structure exists even for the vacuum case. (We should analyze the original equations without matter, which will be done in Appendix C). Therefore, the existence of such a wormhole type solution may be caused by the strong gravity effect rather than the effect of the pressure.
The wormhole throat corresponds to the point dµ/dr g = ∞ (i.e., dr f /dr g = ∞). When we have dr f /dr g = ∞, the interaction terms diverges at the point. As a result, the contribution from the interaction term should not be ignored even for the case with a very small graviton mass, and then our assumption is no longer valid at a wormhole throat. Hence, we have to re-investigate whether a relativistic star does not exist for the coupling constants of Class [I]. We shall analyze it in next section.
When P g (0)/ρ g becomes larger, i.e., if P g (0)/ρ g > 1/3, we again find a real µ 0 , but there exists no regular µ(r) for the whole range of r. µ(r) exists in two separated regions; one is smaller than some finite radius (< R ⋆ ) and the other is larger than another finite radius (> R ⋆ ) , In both regions, two branches A and B are connected. We find a kind of closed universe for the smaller-radius inner region, and a kind of wormhole structure for the larger radius outer region. Both spacetime structures contain a curvature singularity at the throats (the turning points of µ(r)).
On the other hand, for the Class [II] example, both branch A and B solutions exist for any value of P g (0) (Fig. 3), and they are not connected each other. Hence we always find a relativistic star with asymptotically homothetic spacetime structure (the branch A solution).
We note that at the boundary of Class [I] and Class [II] (i.e., β 2 = − √ β 3 ). The trivial solution µ = 1/ √ β 3 gives the branch B. While the branch A has a non-trivial solution shown in Fig. 3 (b), which gives a relativistic star for any value of P g (0).
Hence we may conclude that a relativistic star always exists a regular solution for the coupling constants of Class [II]. On the other hand, there does not exist a relativistic star beyond a critical value of the pressure for the coupling constants of Class [I], i.e., for P g (0)/ρ g > w − . Instead, the spacetime may turn to a wormhole geometry with a singularity (or a closed universe with a singularity).
The existence condition of P g (0)/ρ g < w − can be rewritten by the compactness of a star, GM ⋆ /R ⋆ . Using the internal solution (4.3) and (4.5), we find (4.17) Then we obtain the existence condition for Class [I] as This gives the maximum value of the compactness of a relativistic star for given coupling constants β 2 and β 3 .
Since β 2 and β 3 , are restricted as shown in Fig. 1 which is obtained from the existence condition for a regular interior solution in GR because there is no additional constraint in this class. The upper bound of the compactness in Class [I] is almost the same as the observed value (e.g., the compactness is about 0.3 when a radius of a two solar mass neutron star is 10 km, while it is about 0.21 for a two solar mass star with a radius of 14 km [40][41][42].). In order to give a stringent constraint on the theory by observations, we have to analyze more realistic star, which will be discussed in the next subsection. Giving more realistic equation of state, we present a neutron star solution in the bigravity theory. We then discuss its mass and radius in order to give a constraint on the theory or the coupling constants by comparing them with observed values.
We assume a simple polytropic-type equation of state where we set K = 1.5 × 10 5 [cgs]. In the massless limit of the graviton, we have two decoupled Einstein equations. Then the f -metric is given by the Minkowski spacetime because there is no f -matter, For g-spacetime, we have the same neutron star solution as that in GR. We present ρ c -M ⋆ and R ⋆ -M ⋆ relations in Fig. 4, where ρ c = ρ g (0) is the central density. We find that the maximum mass of a neutron star is about 2M ⊙ , where M ⊙ is the solar mass, for the above equation of state. This result is obtained in GR but also it is the case for Class [II] in bigravity because we always find the regular solution for µ(r) in the whole coordinate range (0 ≤ r < ∞). We show some examples for the same coupling constants (4.15) with several values of the central density ρ c in Fig. 5. However, for Class [I], we find the additional constraint to find the regular µ(r) as we expect from the result in the previous subsection. We also present some examples of µ(r) for the same coupling constants (4.14) with several values of ρ c in Fig. 5. This figure shows there is no regular solution of µ(r) in the whole region if the density ρ c is larger than 2.8 × 10 14 g/cm 3 . This upper limit of the density does not reach the central density with the maximum mass of neutron star in GR (see Fig. 4).
Hence this limit of ρ c provides the maximum mass of a neutron star in Class [I], which is much smaller than that in GR (or in Class [II]). In Fig. 1 In this section, we numerically solve the basic equations under the metric ansatz (3.1) and (3.2) with a gmatter field. We find a relativistic star solution and confirm the previous results obtained in the massless limit when the graviton mass is sufficiently small.
We numerically integrate Eqs. (3.15), (3.19) and (3.20) outwards from the center r = 0. The constraint equation (3.18) is used to evaluate the boundary values at the center. Since it must be satisfied in the region of r > 0 too, we use this constraint to check the accuracy of our numerical solutions in r > 0.
Since the equations are seemingly singular at r = 0, We start our calculations from r = 0 + δr with δr ≪ 1. All variables are expanded around r = 0 as where X (n) (0) is the n-th derivative of the variable X at r = 0.
Here, by use of the freedom of time coordinate rescaling, we choose N g (0) = 1 without loss of generality. 4 We determine the values of variables at r = δr by using up to second order of (5.1).
In this section, we focus only on the branch A solution since we are interested in an asymptotically flat spacetime. We will give some remarks for the branch B, which gives an asymptotically AdS spacetime, in Appendix B. 4 Although it gives Ng(∞) = 1, if we wish to find the boundary condition Ng(∞) = 1, we redefine new lapse functions as and new time coordinate as New metrics defined byÑg,Ñ f andt satisfy the boundary con-ditionÑg(∞) = 1 at infinity.

A. A uniform density star
We first discuss a uniform density star, i.e., ρ g = constant. The dimensionless parameters characterizing the star are where we have defined which gives the effective graviton mass on the homothetic spacetime. The first parameter in (5.4) is evaluated as which is much larger than unity because m −1 eff is the Compton wavelength of the graviton and then it must be a cosmological scale.
Once the parameters (5.4) are given, the proper value of µ(0) is determined by a shooting method to adjust the correct boundary condition (3.5) at infinity as well as the asymptotic flatness. Then all coefficients in Eq. (5.1) are fixed by this µ(0) from the expanded basic equations order by order, We use µ 0 as the center value of µ(0) in the case of massless limit. When the value of the graviton mass is sufficiently small, the proper value of µ(0) is close to µ 0 . Hence, we start to search for µ(0) near µ 0 to find a regular solution with the correct boundary condition.
To check the boundary conditions at infinity, we evaluate the eigenvalues of γ µ ν , i.e., If all eigenvalues approach the same constant K as r → ∞, the solution is asymptotically homothetic. Then the γ energy-momentum tensor will become a "cosmological" constant (Λ g ) term at infinity. We find our solution with an asymptotic flatness, if Λ g = 0, which we have assumed for our coupling constants.

Class [I]
As an example in Class [I], we choose the same coupling constants as before, i.e., Λ g = 0 , m g = m f , β 3 = −3 , β 4 = 3 . The branch A solution approaches an asymptotically flat homothetic spacetime. In Fig. 6, we show a numerical solution by setting κ 2 g ρ g /m 2 eff = 2.5 × 10 5 , 5 for which the typical value of the Vainshtein radius is given by GR is recovered within the Vainshtein radius. We note λ 1 is discontinuous at the star surface R ⋆ . It is because the discontinuity of the matter distribution leads the discontinuity of r ′ f as seen in Eq. (3.20). This discontinuity disappears when we discuss a continuous matter distribution such as a polytropic star (4.21) as shown in Fig. 8.
Changing the central value of the pressure P g (0)/ρ g , we find the solution disappears for P g (0)/ρ g > 0.0665. It is consistent with the argument in the massless limit, in where Pg is the numerical solution with a finite mass and P [m=0] g is the solution in massless limit. We set Pg(0)/ρg = 5 × 10 −2 and (5.10) with κ 2 g ρg/m 2 eff = 2.5 × 10 5 (the red solid curve), κ 2 g ρg/m 2 eff = 2.5 × 10 7 (the blue dashed curve) and (5.12) with κ 2 g ρg/m 2 eff = 2.5 × 10 5 (the green dotted curve). We note Pg − P which the critical value is given by P g (0)/ρ g = 1/15 ≈ 0.06667. Hence even in the case with a finite graviton mass, there exists a critical value of the pressure beyond which a regular star solution does not exist.
If we choose the larger value of the parameter as κ 2 g ρ g /m 2 eff = 2.5 × 10 7 , the solution exists for P g (0)/ρ g > 0.0666, which is closer to the value in the massless limit. Hence, we expect that the massless limit approximation is valid for the realistic value κ 2 g ρ g /m 2 eff ∼ 10 43 . If the solution exists, the inner structure of star as well as the gravitational field are restored to the result of GR because of the Vainshtein mechanism. We find differences between our numerical solution and the semianalytic solution in massless limit are very small as shown one example of the pressure P g in Fig. 7. This fact also confirms the validity of the massless limit approximation if the graviton mass is sufficiently small. We conclude that the bigravity for Class [I] cannot reproduce the result in GR beyond the critical value of P g (0)/ρ g .

Class [II]
As an example in Class [II], we choose one of the previous coupling constants, i.e., Λ g = 0 , m g = m f , β 3 = 1 , β 3 = 3 (5.12) and we set κ 2 g ρ g /m 2 eff = 2.5 × 10 5 . In this case, we can find a regular star for any values of P g (0). The solution is almost the same as the massless limit (or GR) as shown in Fig. 7. We conclude that in the bigravity theory in Class [II] the results in GR are recovered and the Vainshtein mechanism holds even in a strong gravity limit. For a neutron star with a realistic equation of state, we can also confirm the above results, i.e. the massless limit is valid. Here we again assume the polytropic equation of state (4.21).
One typical example of the solutions in Class [I] is shown in Fig. 8, where we choose the coupling constants as (5.10) and ρ c = 1.71 × 10 14 g/cm 3 , m −1 eff = 10 4 km , (5.14) We find a neutron star solution with which is the same as those in the massless limit. Our numerical calculation shows that increasing the central density ρ c , the solution exists only for M ⋆ < ∼ 0.882M ⊙ for the coupling constants (5.10). We have obtained M ⋆ 0.886M ⊙ in the massless limit. If we choose the larger value of the Compton wave length of the graviton as m −1 eff = 10 5 km, the mass upper limit increases as M ⋆ 0.884M ⊙ , which is closer to the value in the massless limit.
For Class [II], we always find the same solution as that in GR. As a result, as the case of a uniform-density star, we confirm that the massless limit solution is a good approximation for the sufficiently small graviton mass.

VI. CONCLUDING REMARKS
Assuming static and spherically symmetric spacetimes, We have presented a relativistic star solution in the bigravity theory. For simplicity, we have considered only gmatter fluid and given only asymptotically flat solutions in the text. Some solutions with the other conditions are discussed in Appendix B.
First we obtain the solutions under the massless limit approximation in Sec. IV. Then, by solving the basic equations numerically without the approximation in Sec. V, we confirm such an approximation is valid since the graviton mass, if it exists, must be sufficiently small. We find that the coupling constants are classified into two classes: Class [I] and Class [II]. For both classes, the Vainshtein screening is found in the weak gravitational field. However, when we take into account a relativistic effect, the Vainshtein screening mechanism may not work in some strong gravity regime in Class [I]. In fact, to find a regular function of µ(r) in Class [I], the central pressure is constrained, and as a result, the maximum mass is much smaller than that in GR as shown in Fig. 1. Beyond this maximum mass, the Vainshtein mechanism does not work well since GR solution is not obtained.
On the other hand, there is no additional constraint for Class [II], and the structure of star as well as the gravitational field are restored to those in GR for the expected small graviton mass. The Vainshtein screening mechanism works well in Class [II].
In Table I  The result suggests that Class [II] is favored from the existence condition of a neutron star. As the necessary condition of Class [II], the parameters should satisfy as shown in Fig. 1. However, those parameters should happen to satisfy from the cosmological point of view, which constraint comes from to find a stable solution in the early Universe in bigravity [17]. There is no intersection of the parameters because the boundaries of Class [II] and of the cosmological constraint coincide exactly. If we take the parameters in Class [I] from the cosmological constraint, the equation of state of the star will be strongly constrained to find a two solar mass neutron star. Conversely, if we assume Class [II] from the astrophysical point of view, the problem of ghost or gradient instability may reappear in the early Universe. There is another problem in Class [II] parameters. Since we have started to discuss the bigravity theory in order to explain the present acceleration of the Universe, the parameters (or coupling constants) should predict the existence of a positive effective cosmological constant (Λ g > 0). If we impose the same conditions on the coupling constants as discussed in [11], {b i }'s are given by two coupling constants c 3 and c 4 . The existence condition of de Sitter solution as well as Minkowski solution yields which excludes the possibility of (6.1). Hence, if we assume the Minkowski spacetime is a vacuum solution, Class [II] cannot admit the de Sitter solution as another vacuum solution as well, thus the acceleration of the Uni-verse cannot be explained in the bigravity. In this paper, we have assumed that both static g-and f -spacetimes are static with respect to the same time coordinate t, and the Stückelberg field µ is also static. However there is a possibility such that the existence of the critical value in Class [I] might be caused by the above simple ansatz. The static ansatz of the Stückelberg field may not be necessary to obtain an (approximate) static spacetime. In fact, in the case of cosmology, a homogeneous configuration of the Stückelberg field leads an instability, while the inclusion of an inhomogeneity in the Stückelberg field gives a stable solution, which describes an (approximate) homogeneous spacetime due to the Vainshtein screening [17]. Hence, to draw a final conclude about the existence of a massive neutron star (and also a black hole solution), relaxing the static ansatz of the Stückelberg field, we should extend our analysis to the spacetime with dynamical Stückelberg fields, which we leave for our future work.
with the following conditions: where a prime denotes the derivative with respect to r. From the basic equations, we find a septic equation for µ as where C m 2 , C Λ and C matter are explicitly defined in [17]. These terms have typical magnitudes given by and the last term is given by where M ⋆ is the gravitational masses of the g-matter.
There is a root of Eq. (A5) with µ → 0 as r → ∞, which is the asymptotically homothetic branch. Such a branch should be extended inward without any singularity. As discussed in [23,24], the branch with µ = 0 at r = ∞ reaches to µ → −1/ √ β 3 in the range of r ≪ R V , where we find a successful Vainshtein screening.
Although we cannot find analytic roots µ(r) of the septic equation (A5), we can easily find a inverse function r(µ) because r appears only in C matter as the form (A6). The result indicates that the function r(µ) is a singlevalued function. However, the function µ(r) is not a single-valued function, if there is an extremal value of the function r(µ), i.e., dr/dµ = 0. The point of dr/dµ = 0 corresponds to a curvature singularity. Hence a regular solution must be given by a monotonic function µ(r) in the domain R I < r < ∞, where R I is a typical length, if it exists, below which the weak gravity approximation is not valid.
As discussed in [24], we find the parameter constraint as follows: Since the function µ(r) should be monotonic, the function is approximated by with 1 ≫ δµ > 0 in r ≪ R V . Substituting this expression into (A5), we find Since the right hand side is negative, the necessary condition is given by where However the constraint (A9) is not sufficient, because it does not guarantee that the function µ(r) is a singlevalued function in the domain R I < r < ∞, which is guaranteed by r(µ) has no extremal value in −1/ √ β 3 < µ < 0. We must impose dr(µ)/dµ > 0 for any µ with Three examples of the solution µ(r) are shown in Fig. 9: (a) β 2 = −3, β 3 = 3, (b) β 2 = 1.73, β 3 = 3, and (c) β 2 = 7, β 3 = 3. The case (a) and (b) satisfy while the case (c) satisfies For both (a) and (b), the branch of µ ≃ −1/ √ β 3 in r ≪ R V connects the branch of µ = 0 at r = ∞. However, the case (a) gives the single-valued function µ(r), while the case (b) is not. It indicates that the ratio of two radial coordinates are not single-valued function 6 . For the case (c), there are two curves (c-1) and (c-2) and these are disconnected. Note that, the branch (c-2) can be extended to infinity. This branch is not an asymptotically Minkowski solution, but an asymptotically AdS solution similarly to the branch C which will be discussed in Appendix B 3. As a result, the parameter constraint is approximately given by as shown in Fig. 10. The hatched light-blue region gives a successful Vainshtein screening solution. We can show numerically that there is no regular asymptotically homothetic solution in the narrow region along β 2 = √ β 3 (the red region), in which µ(r) is not a single-valued function such as (b) in Fig. 9, and should then be excluded.
ity. For the ansatz (3.1) and (3.2), we find the eigenvalues {λ 0 , λ 1 , λ 2 , λ 3 } of γ µ ν as Then, the γ energy-momentum tensor is given by We then find in following three cases that the γ energymomentum tensor turns to be a cosmological constant: Case (i) Case (ii) Case (iii) We note that the equation b 1 + 2b 2 λ 2 + b 3 λ 2 2 = 0 is equivalent to where we use λ 2 = K(1 + µ). Case (i) gives an asymptotic homothetic spacetime, i.e., an asymptotic de Sitter or anti-de Sitter spacetimes as well as an asymptotic Minkowski spacetime. In addition, as we will show in the next subsection, we also find a solution with a cosmological constant given by Case (ii).

Relativistic star with g-matter
Just for simplicity, we discuss a uniform-density star only with g-matter fluid. We use the parameters (4.14) as an example for Class [I], and parameters (4.15) for Class [II]. We then choose In the text, we consider the branch A without a cosmological constant, in which case, the branch A solution approaches the Minkowski homothetic spacetime. Here, we discuss asymptotic structures of branch A when we introduce a non-zero cosmological constant.
For the branch A, the results are the same both in Class [I] and in Class [II]. When we introduce a negative cosmological constant, the solution approaches the homothetic anti-de Sitter spacetime at infinity as shown in Fig. 11 (N g /N f , F g /F f , r g /r f → 1). For a positive cosmological constant, when 2Λ g 3m 2 eff (the Higuchi bound) is satisfied, the solution seems to approach a homothetic de Sitter spacetime. Since we cannot solve the basic equations beyond the cosmological horizon, we cannot conclude definitely that the solution is asymptotically homothetic, but as shown in Fig. 12, the solution seems to approach a homothetic spacetime because the eigenvalues coincide around r ≈ m −1 eff before the horizon. However, if 2Λ g 3m 2 eff , a regular solution disappears as discussed in the appendix of [17]. As a result, the branch A always approaches a homothetic spacetime if the cosmological constant satisfies 2Λ g 3m 2 eff .  Fig. 13, where ℓ AdS := −3/Λ g is the AdS curvature radius. Note that this solution is not asymptotically homothetic. The eigenvalues λ 0 and λ 2 approach the same value with satisfying 1 + 2β 2 µ + β 3 µ 2 = 0, for which the interaction term becomes just a cosmological constant as discussed in Appendix B 1. Although the g-and f -spacetimes are not homothetic at infinity, both spacetimes approach asymptotically to some AdS spacetimes.

Relativistic star with f -matter
Here, we discuss the effect of the f -matter field. For simplicity, we assume ρ f ≫ ρ g , for which we regard that the g-spacetime is almost vacuum.
The action of the bigravity is symmetric for g-and f -spacetimes under the transformation Then the case only with f -matter is equivalent to the case only with g-matter for corresponding coupling constants under the transformation (B12), i.e.,   [17]. Although the result presented in [17] is only the case of Class [I], we find the same behaviour even for Class [II]. One exceptional case is a homothetic solution.
there exists a homothetic solution, i.e., N g = N f , F g = F f and µ = 0, for which the solution is identical to that in GR in the whole space region.

a. Massless limit approximation
In the massless limit, the interior solution is given by FIG. 14: The same figure as Fig. 2 in the case of f -star.
where we assume a uniform density for f matter fluid.
The g-spacetime is just a Minkowski solution. The exterior solution is given by where we define the gravitational mass by and R ⋆ is the radius of the f -star measured in fspacetime. Similarly to the argument in Sec. IV, the ratio must be The center value of µ is given by a root of thus there are two branches (the branch C and D) similar to the case of g-star. The branch C approaches a homothetic spacetime as we will see later.
We chose the coupling constants as (5.10) in Class [I]. The solution in the massless limit is shown in Fig. 14. For the case of the f -star, the wormhole geometry is not found. Now we solve the basic equations for each branch without the massless limit approximation.
b. Branch C We set and ρ f = constant. For the above parameter setting with K = 1, an asymptotically AdS solution is found for −0.76m 2 eff Λ g (1) 0.05m 2 eff . This solution in the branch C is asymptotically homothetic because the eigenvalues λ 0 , λ 1 , λ 2 converges to the same constant although its value is not unity as shown in Fig. 15.
The reason is as follows: When we fix parameters {m g , m f , Λ g , K, β 2 , β 3 }, the original coupling constants {κ f , b i } are determined. Once the original coupling constants are given, all homothetic solutions given by are characterized by the proportional factorK which is one of the roots of the quartic equation In the range of −0.76m 2 eff Λ g (1) 0.05m 2 eff , there are four real roots forK. For instance, when we set Λ g (1) = 0, we find K = −0.604, 1, 1.44, and 3.83 , and find four homothetic solutions (one Minkowski, one de Sitter, and two AdS spacetimes). It turns out that the solution we solved approachesK = 1.44 homothetic spacetime. Since Λ g (1.44) < 0, it is the asymptotically AdS spacetime. Note that when we assume Λ g (1) −0.76m 2 eff , there are only two real roots ofK, e.g, K = −0.586, and 1 , for Λ g (1) = −m 2 eff . In this case, we cannot find a regular solution for the branch C in Λ g (1) −0.76m 2 eff . In the case of 3m 2 eff /2 ≫ Λ g (1) 0.05m 2 eff , there are four homothetic solutions, e.g., K = −0.621, 1, 1.11, and 4.85 . (B32) for Λ g = 0.1m 2 eff . The solution may approach thẽ K = 1.11 homothetic solution with Λ g (1.11) > 0. However, because of a numerical instability, we cannot confirm that there is a regular solution approaching de Sitter spacetime for 3m 2 eff /2 Λ g (1) 0.05m 2 eff . Finally, we give a comment for the case of Λ g (1) 3m 2 eff /2. In this case, the Jacobian J = dr f /dr g diverges before reaching the cosmological horizon. Therefore, this solution has the curvature singularity as discussed in Appendix C.
c. Branch D For the branch D, we cannot construct any regular solution with or without a cosmological constant by our numerical approach. Although the solution is regular below the Vainshtein radius, there is a singularity at a radius near the Compton wavelength of the massive graviton. Thus we will not discuss the branch D furthermore.

Appendix C: Wormhole-type solution
In Class [I], as shown in Fig. 2 (b), we cannot find a regular solution beyond the critical value of the pressure. The solution turns to a closed spacetime or a wormholetype spacetime beyond the critical value.
In this appendix, we shall discuss what kind of wormhole type structure is obtained in the bigravity theory. To find a solution with a wormhole-type structure, we should integrate the basic equations from the wormhole throat. As mentioned in the subsection IV B, a wormhole throat corresponds to the point of J = ∞, where the function J is the Jacobian for the radial coordinate transformation from r g to r f . When we find J = ∞ at some radius, such a coordinate transformation is singular. That is, we cannot define the transformation r g → r f at the point. Similarly, we cannot define the transformation r f → r g at the point of J = 0.
When the coordinate transformation r f = r f (r g ) is not well-defined (i.e., J = ∞) at some point, we cannot integrate beyond such a singular point as a function of r g . However, the inverse function r g = r g (r f ) is well-defined at J = ∞. As a result, we can solve the equations and find the solution as a function of r f by using the radial coordinate r f , i.e., the basic equations to be solved are Although the point of J = ∞ is a curvature singularity as we will see, we can continue to solve the equations and find the solution beyond such a singularity. For simplicity, we assume vacuum spacetimes, i.e., there is neither g-matter nor f -matter. A wormhole throat of g-spacetime is given by J −1 = dr g /dr f = 0 , at which we assume the variables N g , F g , N f , F f , µ are finite. Setting the radial coordinate as r = r f , we find the derivatives of g-variables are finite at J −1 = 0 because Eqs. (3.11) and (3.12) yield for J −1 → 0. Furthermore, Eqs. (3.13) and (3.14) indicate that the derivatives of f -variables are also finite at J −1 = 0, and Eq. (C2) indicates dµ/dr f is finite. Hence, the first derivatives of all variables are finite even at J −1 = 0. Since the differential equations are first order, we can solve the equations numerically beyond J −1 = 0 by use of the r f coordinate. Since two metric are symmetric in the bigravity theory, the above argument is also applied to the point of J = 0, which is a wormhole throat in f -spacetime, At J = 0, the coordinate transformation r g = r g (r f ) is not welldefined, but the solution is obtained as a function of r g beyond this singularity.
In the case of Λ g = 0, the branch B solution contains a singularity at some radial point. To find a regular wormhole-type solution, we should introduce a negative cosmological constant. Here we set the parameters as and Λ g = −75m 2 eff (ℓ AdS = 0.2m −1 eff ) .
We first use the g-radial coordinate r g . Suppose that a wormhole throat exists in the f -spacetime (which we call the f -throat), so J = 0 at a radius r g = a f The value of N g on the throat is arbitrary by the rescaling freedom of the time coordinate, and the value of F g gives the gravitational field strength at the throat, which characterize the property of the wormhole. Since we have two algebraic equations at the f -throat as where C is the constraint equation defined by Eq. (3.18), when we give the values of F g and N g at r g = a f , the values N f (a f ), F f (a f ) are determined by Eqs. (C7) as functions of µ(a f ). We first solve variables outward on the r g coordinate system, and find an asymptotically homothetic AdS spacetime by tuning the value of µ(a f ). Next, we solve variables inward with respect to the r g coordinate. When we find the point of J −1 = 0 at a radius r g = a g , which is the wormhole throat in g-spacetime (the g-throat) 7 . we cannot continue to integrate the basic equations numerically on the r g coordinate. Then we switch the radial coordinate from r g to r f , and solve variables with respect to the r f coordinate beyond the point of J −1 = 0. Finally we find a global wormhole-type solution, which example is given in Figs. 16  µ(a f ) is tuned as µ(a f ) = 0.03847, which gives the asymptotically AdS spacetime. Here we have introduced a typical length scale of the wormhole r S by where we define a mass function M g (r) by Fig. 16 shows the relation between two radial coordinates. The top panel gives r f /r g in terms of r g coordinate. It shows that has r f /r g takes two different values at the same radius r g . One branch (r f /r g → 1) approaches the homothetic AdS spacetimes, while another branch (r f /r g → 1.183) approaches the non-homothetic AdS spacetime. Two different asymptotic structures are connected by the wormhole. Fig. 16 shows that the g-throat and the f -throat are located at the different points.
We depict the Ricci curvature scalar of the f -metric as well as one of the g-metric in Fig. 17, where we have used the variable r f /r g to parametrize the radial coordinate, instead of either r g or r f , because either coordinate r g or r f is not a single-valued function near the throats. The g-throat (J −1 = 0) is located at r f /r g = 1.2563 and the f -throat (J = 0) is founded at r f /r g = 1.03847. The Ricci curvature scalar of the g-metric diverges at the g-throat. It is caused by the divergence of the γ energymomentum tensor at the wormhole throat. As shown in Fig. 17, Ricci scalar goes to +∞ as r f /r g → 1.2563 − ǫ, while it goes to −∞ as r f /r g → 1.2563 + ǫ with 0 < ǫ ≪ 1. Note that f -spacetime curvature is finite even at the g-throat of J −1 = 0. Only the g-spacetime Ricci curvature diverges. Inversely, only the f -spacetime Ricci scalar diverges at the f -throat. This behaviour is quite similar to the case of the cosmology [11].
Finally, we discuss the Vainshtein screening. Since the γ energy-momentum tensor cannot be ignored at the throat point, the Vainshtein screening mechanism is no longer guaranteed. We may find a deviation from the GR result. In fact, the geometry of the vacuum spacetime turns to a wormhole geometry, which does never appear in GR. To see the differences of the metric functions from GR, we show the variation rates of the mass function M g and the ratio N g /F g in Fig. 18. In GR, two functions are exactly constant. In the bigravity, although two functions are not exactly constant, these are almost constant. Hence, the metric functions are well-approximated by the Schwarzschild-AdS metric (up to their first derivatives) although the topology of the solution is different from the Schwarzschild-AdS spacetime.