Fermion Proca Stars: Vector Dark Matter Admixed Neutron Stars

Dark matter could accumulate around neutron stars in sufficient amounts to affect their global properties. In this work, we study the effect of a specific model for dark matter -- a massive and self-interacting vector (spin-1) field -- on neutron stars. We describe the combined systems of neutron stars and vector dark matter using Einstein-Proca theory coupled to a nuclear-matter term, and find scaling relations between the field and metric components in the equations of motion. We construct equilibrium solutions of the combined systems, compute their masses and radii and also analyse their stability and higher modes. The combined systems admit dark matter (DM) core and cloud solutions. Core solutions compactify the neutron star component and tend to decrease the total mass of the combined system. Cloud solutions have the inverse effect. Electromagnetic observations of certain cloud-like configurations would appear to violate the Buchdahl limit. This could make Buchdahl-limit violating objects smoking gun signals for dark matter in neutron stars. The self-interaction strength is found to significantly affect both mass and radius. We also compare fermion Proca stars to objects where the dark matter is modelled using a complex scalar field. We find that fermion Proca stars tend to be more massive and geometrically larger than their scalar field counterparts for equal boson masses and self-interaction strengths. Both systems can produce degenerate masses and radii for different amounts of DM and DM particle masses.


I. INTRODUCTION
The nature of dark matter (DM) is one of the large remaining open questions in physics.Even though it constitutes roughly 26.8% of the total energy density of the universe [1] and has a long observational history [2], its properties remain largely unknown.We currently know that DM likely is a particle that is only interacting gravitationally and weakly with standard model particles, and that is invisible through electromagnetic radiation.Large-scale structure formation in the universe further suggests that DM is mostly cold, i.e., slowly moving [2][3][4][5].This makes it an integral part of the standard model of cosmology.
Neutron stars (NSs) are used to probe a large range of physical phenomena.They are dense and compact remnants of heavy stars.
Their high densities make them excellent laboratories for probing gravitation and nuclear physics under extreme conditions.They are characterized using the nuclear matter equation of state (EOS).The EOS describes the relation between pressure and energy density of the matter found inside NSs.It is needed to close the system of differential equationsthe Tolman-Oppenheimer-Volkoff (TOV) equations [6,7] -that describe the density distribution of a spherically symmetric static NS and the spacetime curvature.A significant constraint on the EOS is the ability to produce NSs with masses larger than two solar masses, 2 M ⊙ .The most massive NS known to date is PSR J0952−0607 with a mass of M = 2.35 +0.17  −0.17 M ⊙ [8].The lighter companion of the binary system observed in the GW190814 gravitational wave event [9] was also proposed to be the heaviest NS, with a mass of around 2.6 M ⊙ .But there is evidence that it might be the lightest known black hole instead [10].High maximum NS masses require stiff EOS, where the nuclear matter is difficult to compress and the energy density rises sharply with increasing pressure.Other constraints include the measurements of the pulsars PSR J0030+0451 [11] and J0740+6620 [12] by the NICER telescope.They also favor a stiff EOS.In contrast, the gravitational wave event GW170817 [13,14] favors soft EOS which produce smaller NSs that are more compact and more difficult to tidally disrupt.
Additionally, it has been proposed to probe the DM properties using NSs.For example, DM can form a cloud or accumulate inside NSs as a core.In sufficient amounts, it can modify the NS properties such as mass, radius and tidal deformability.These properties have been measured using telescopes such as NICER and the gravitational wave detectors LIGO, Virgo and KAGRA.This allows us to probe the properties of DM such as its particle mass and self-interaction strength (see, e.g., [15][16][17][18][19][20][21]).There exist numerous candidates for DM particles.A possible DM candidate is an additional bosonic field (scalar field or vector field), as was studied in [22][23][24][25][26].
The idea that an astrophysical object consists of a mixture of fermionic and bosonic matter goes back to [27,28].A multitude of different models of these fermion boson stars (FBSs) have since been investigated (see, e.g., [18,[29][30][31][32] for reviews).In the simplest case, the fermionic and bosonic components interact only gravitationally through the effects of their matter-energy content and without an explicit coupling (i.e., they are minimally coupled).This makes FBSs interesting objects in the context of DM research (see, e.g., [15,33,34]).They have been studied in connection to NSs, where the NS provides the fermionic component and a bosonic field provides the bosonic component of the FBS [15,33].The bosonic component can be modelled via, e.g., scalar and vector fields.Another possibility is to study the bosonic field as a fluid or gas of particles.Which treatment is used depends mainly on the mass range of the supposed DM particle.With DM masses above the eV scale, it is treated as a collection of particles or as a fluid.For DM masses ≪ eV the correlation length becomes larger than the average particle separation.Dark matter is then best described as a macroscopic wave.We follow the second approach in this work.We refer to [35] for a review of observational prospects of DM at different mass scales.
FBSs have been studied with regard to their stability [28].Their dynamical properties were explored in [36][37][38][39][40][41].Numerical simulations aiming to understand the gravitational wave signals were performed by [41].In all these cases, the NS component was modelled using a perfect fluid and a classical complex scalar field was used for the bosonic component.However, understanding vector fields is equally relevant for a number of reasons.If DM is a spin-1 particle, it would be described using a vector field.Some theories of modified gravity also feature vector fields with similar behavior [42][43][44][45][46].In these cases, the vector field is usually directly coupled to the curvature.But if the coupling is weak enough, these fields could behave very similar to minimally coupled fields.In this work, we therefore explore the effect of minimally coupled vector fields on NSs.
Fermion boson stars can form in a variety of ways.But in essence, the problem reduces to how one can accumulate a large amount of scalar or vector fields in and around a NS.One common motivation for these fields is bosonic dark matter.It could arrange itself around NSs as a cloud or inside NSs as a core.
Points 3) and 4) in particular are highly model dependant.Accretion and decay rates depend on the DM particle interaction cross sections and available decay channels.Previous works, for example [57], have shown, for the case of non-interacting scalar DM, that old isolated neutron stars set strong bounds on the allowed scattering cross section between light quarks and DM.This could in practice strongly disfavour the accumulation of significant amounts of DM in NS through accretion.NSs with clouds could form in a similar way, given that either the DM is the dominant contribution to the FBS or that the DM properties only allow low-compactness configurations (e.g., when the particle mass is small [15]).The fermionic and bosonic components could conceivably be separated from one another, e.g., during a supernova NS-kick [58][59][60][61].There, the stellar remnant gets ejected and rapidly moves away from the remaining stellar envelope.This process could allow for NSs with a large range of possible DM-fractions.The DM particles most interesting for FBSs are generally (self-interacting) ultralight DM particles, weakly interacting massive particles, dark photons [24,25] (as a candidate for vector DM) and axions [15,35,[62][63][64][65][66][67][68][69][70].
Another formation channel is motivated through theories of modified gravity.One way of producing large amounts of scalar (or vector) fields is superradiance [42,71].Spontaneous scalarization [30,72] also provides a way of producing significant scalar [30,73] and vector1 [43,44] field amplitudes.It has also been studied explicitly in NSs [45,46,72] and could be a way of forming systems with scalar and vector fields.Scalarization might also take place dynamically in the late stages of the evolution of binary NS systems [74], forming either a black hole or a FBS after merger (depending, e.g., on the initial masses of the binary objects).
Self gravitating vector fields have already been investigated.These objects are called Proca stars.They are modelled by a complex vector field and were first proposed by [75].They can be thought of as macroscopic condensates of spin-1 particles [30].Proca stars have been studied by a number of groups analytically [76][77][78] and numerically [79,80], such as in merger simulations [81,82].Different types of Proca stars with charge [83], rotation [75] and with a quartic self-interaction potential [84] were also considered.Other works [85][86][87] studied shadow images of Proca stars in different scenarios.
In this work, we study the combined system of a vector field and NS matter, which we call fermion Proca stars (FPSs).Starting with an action for complex vector fields coupled minimally to gravity and nuclear matter, we derive a system of differential equations and solve them numerically (section II A).We also pedagogically motivate the boundary conditions (section II B), find an analytical bound for the vector field amplitude and derive scaling relations in the equations of motion (section II C).The equations are solved using a shooting method and the integrator implemented in our code (for the code, see [88]).The numerical methods are also explained in section II E. We show radial profiles of FPSs (section III A) and then compute global quantities such as mass and radius and compare them to astrophysical observations (section III B).In section III C, we compare FPSs to their counterpart with a scalar field.In the following, we refer to the scalar case as "fermion boson stars" (FBSs).Finally, we compute higher modes of FPSs and compute configurations with different EOS (section III D).We find that the vector field significantly affects the NS properties and thus produces detectable signatures.FPSs admit DM core and cloud solutions.Small DM masses lead to DM clouds, and large masses form DM cores.Core solutions compactify the NS component.Cloud solutions lead to less compact configurations.Some solutions appear to violate the Buchdahl limit when only observing the NS component.We then compare FPSs (with a vector field) to FBSs (with a scalar field).FPSs tend to be more massive and geometrically larger than FBSs for equal boson masses and self-interaction strengths.For a given measurement, this would favor larger vector DM masses (compared to scalar DM), because larger DM masses produce smaller and less massive objects.We find a significant amount of degenerate solutions between different choices of FBSs, FPSs, the DM properties and the EOS.For different boson masses and DM-fractions, FPSs and FBSs can both be degenerate with each other and also be degenerate with pure NSs with a different EOS.Using scaling relations for pure boson stars and Proca stars, we show that FBSs and FPSs are virtually indistinguishable if the boson masses differ by a factor of 1.671 and the DM has no self-interactions.We confirm the existence of FPSs in higher modes which are stable under linear radial perturbations.
Throughout this work, we use units where G = c = M ⊙ = 1 (also see Appendix A).The Einstein summation convention for tensors is implied.This paper is based on the Master thesis of Cédric Jockel [89].

A. Equilibrium Solutions
Fermion Proca stars (FPSs) are combined systems of fermions and vector bosons, which interact only gravitationally.They can be seen as a macroscopic Bose-Einstein condensate which coexists with a NS at the same point in space.We model FPSs using a relativistic fluid for the NS component and a complex vector field for the bosonic component.FPSs are described by the Einstein-Proca system minimally coupled to a matter term L m , where R is the Ricci curvature scalar, g is the determinant of the spacetime metric g µν and κ = 8πG/c 4 is a constant.The bar denotes complex conjugation.F µν = ∇ µ A ν − ∇ ν A µ is the antisymmetric field strength tensor and V (A ρ Āρ ) is the vector field potential.The latter depends solely on the magnitude of the vector field A ρ Āρ .By taking the variation of Eq. ( 1) with respect to the inverse spacetime metric δg µν , one obtains the Einstein equations where µν are the energy-momentum tensors describing the NS matter and the vector field matter, respectively.The energy-momentum tensor of the NS matter is taken to be that of a perfect fluid: P and e are the pressure and the energy density of the fluid, respectively.The energy density e is related to the rest mass density ρ through e = ρ(1 + ϵ), where ϵ is the internal energy.u µ is the four-velocity of the fluid.The energy-momentum tensor Eq. ( 3) and the fluid flow J µ := ρu µ are conserved (implying conservation of energy-momentum and of the rest mass, respectively).This leads to the conservation equations The conservation of the fluid flow J µ allows us to define the conserved total rest mass of neutron matter, which we call the fermion number N f .We obtain the fermion number by integrating the right part of Eq. ( 4) over space, The energy-momentum tensor of the vector field is given by where the derivative of the potential V is The equations of motion (Proca equations) of the vector field and the complex conjugate are computed from the action Eq. ( 1) using the Euler-Lagrange equations for a complex vector field.One obtains The covariant derivative of Eq. ( 8) is zero, i.e., ∇ µ ∇ ν F µν = 0.This leads to a dynamical constraint on the field derivative, resembling the Lorentz condition used in the Maxwell and Proca equations (also see [30,75]): This constraint could be useful in numerical simulations to track the numerical error and assess constraint violations of a given numerical scheme.The global U (1)-symmetry in the Lagrangian Eq. ( 1) under the transformation of the vector field A µ (and Āµ ) gives rise to a conserved Noether current The conserved quantity (i.e., the Noether charge) associated to Eq. ( 10) is obtained by integrating the conservation equation ∇ µ j µ = 0 over space, N b is called the boson number and is related to the total number of bosons present in the system.It can equivalently also be interpreted as the total rest mass energy of the bosonic component of the FPS.
We proceed by solving the Einstein equations Eq. ( 2) and the Proca equations Eq. ( 8) for spherically symmetric and static configurations in equilibrium.For that, we consider the spherically symmetric ansatz for the spacetime metric We further assume the perfect fluid to be static, such that the four-velocity can be written as For the vector field, we employ the harmonic phase ansatz and a purely radial vector field (see [75,78,83,84,90]).The vector field is then given by where ω is the vector field frequency and E(r), B(r) are purely radial real functions.
Using the spherical symmetric metric ansatz Eq. ( 12) together with the harmonic phase ansatz Eq. ( 14) for the vector field, we solve the Einstein equations and obtain the equations of motion.One obtains an expression for the radial derivative of a(r) by rearranging the tt-component of Eq. ( 2).We then divide the tt-and rr-components of Eq. ( 2) by α 2 and a 2 , respectively.We add both terms and find a direct relation between the first radial derivatives of a(r) and α(r).We use this to solve for the derivative of α(r).
The evolution equations for the vector field components can be computed from the Proca equations Eq. ( 8).It does not matter which equation of Eq. ( 8) is used since the complex phase will cancel out and will leave only the radial functions in both cases.The ν = r component yields the equation of motion for E(r).The ν = t component of Eq. (8) gives us the equation of motion for B(r).Finally, the r-component of the conservation equation for the energy-momentum tensor (left side of Eq. ( 4)) provides a differential equation for the pressure P (r).For a more detailed derivation, we refer to [89].The full equations of motion for the Einstein-Proca system coupled to matter are thus: This system of equations is closed by providing an equation of state P (e) (or P (ρ, ϵ)) for the nuclear matter part.Note that all equations are first-order differential equations.This is different to scalar FBSs where an additional variable has to be introduced to make the system first-order (see, e.g., [15,33]).Another difference is that no derivative of the potential enters the equations of motion for the metric components in the scalar field case, but it does for the vector field case.For the considered system and ansatz for the metric Eq. ( 12) and vector field Eq. ( 14), the expressions for the fermion number Eq. ( 5) and boson number Eq. ( 11) simplify to R f denotes the fermionic radius (i.e., the radius of the NS component).It is defined by the radial position at which the pressure P of the NS component reaches zero.It is also possible to define the bosonic radius R b .It is defined as the radius where 99 % of the bosonic restmass N b (see Eq. (16b)) is contained.Using these definitions we gain the ability to discriminate between DM core and cloud solutions.Core solutions have R f > R b and cloud solutions have R f < R b .The total gravitational mass is defined in the limit of large radii, imposing that the solution asymptotically converges to the Schwarzschild solution

B. Initial Conditions
We derive the boundary conditions of equations Eq. (15a)-Eq.(15e) at r = 0 and at r = ∞.The values at the origin will later serve as initial conditions for the numerical integration.We first consider the equations of motion in the limit r → 0 while imposing regularity at the origin (i.e., the solution must not diverge).We first analyze equation Eq. (15a).The term proportional to 1/r dominates at small radii and will diverge if r → 0. Thus, the only way to maintain regularity is to set a(r = 0) = 1.It directly follows that a ′ (r = 0) = 0. Similarly, equation Eq. (15b) leads to α ′ (r = 0) = 0.The exact value of α(r = 0) = α 0 is a priori undetermined and can be chosen in a way thought suitable.We will elaborate on this in section II B. The initial conditions for the vector field components E(r) and B(r) can be obtained in a similar manner.We first consider Eq. (15d).In the limit r → 0, the term proportional to 1/r dominates and regularity then demands that E ′ = ωB.It follows that B ′ (r = 0) = 0.This result can be inserted into Eq.(15c), which leads to the relation Since at r = 0, α(r = 0) ̸ = 0 and V ′ ̸ = 0 in general, this relation can only be fulfilled if we demand that B(r = 0) = 0. Plugging this relation into Eq.(15c) yields E ′ (r = 0) = 0.The central value of the field E ′ (r = 0) = E 0 is therefore undetermined by the equations of motion, and thus is a free parameter of the theory.
A similar analysis at large distances reveals the boundary conditions at r → ∞ for all variables.We impose an asymptotically flat spacetime.This requires that a(r → ∞) = α(r → ∞) = 1.All terms proportional to r in Eq. (15a) and Eq.(15b) must vanish at infinity to fulfill the flat-spacetime limit.Therefore, the vector field components must vanish at infinity, E(r → ∞) = 0 and B(r → ∞) = 0. Pressure P (r), energy density e(r) and rest mass density ρ must be zero outside the NS component of the FPS.This will happen at the fermionic radius R f .We summarize all boundary conditions in the following: The initial condition for the metric component α(0) = α 0 is fixed by its behavior at infinity.We also note a popular and widely employed alternative to the self-consistent wave-treatment of DM performed in this work -the two-fluid formalism (see e.g.[21,91,92]).In there, the nuclear matter and the DM are both modeled as perfect fluids, which only interact gravitationally.One can then solve the Einstein equations and obtain a set of modified TOV equations that describe the density distribution of both fluids.Although simplistic, this formalism has the advantage that it is easy to implement, numerically cheap and is applicable to a wide range of fermionic and bosonic DM models.It is also possible to use arbitrary effective EOS for the dark matter.A disadvantage of this model is that it ignores possible wave properties of ultralight DM, that we study in this work.It is also only possible to study nonexcited ground states of the wave-like DM.Further, some emerging properties like the maximal bound on the vector field amplitude (see Eq. ( 24)) would not be captured by the two-fluid formalism.

C. Analytical Results
For a scalar (fermion) boson star, one can scale the field frequency ω to absorb the initial value of α 0 so that it may be set to one (see, e.g., [15]).We investigate whether a similar scaling relation also exists for FPSs.We find that the equations of motion Eq. (15a)-Eq.(15e) are invariant when simultaneously scaling the following variables as The potential V (A ρ Āρ ) is always invariant with respect to this scaling because The invariance of Eq. (15a)-Eq.(15e) under the scaling relation Eq. ( 20) thus allows us to choose σ in such a way that the initial condition for α(0) = α 0 may be set to α 0 = 12 .We will make use of this relation in the numerical analysis.All pre-scaling physical values can be recovered from the asymptotic behavior of α(r → ∞) by performing the inverse transformation to Eq. ( 20).Note that the expression for total gravitational mass Eq. ( 17) is not affected by this scaling.
In contrast to the scaling relation of boson stars with a scalar field, where only the frequency ω and the metric component α are re-scaled, the vector field component E is also affected in the case of Proca stars.To our knowledge, this is the first time the scaling relation Eq. ( 20) has been mentioned explicitly (apart from the Master thesis [89] which precedes this work).[84] briefly mentioned scaling the frequency but not the vector field component.
We also report an analytical bound on the central vector field amplitude E(0) = E 0 .Equations Eq. (15c) and Eq.(15d) govern the dynamics of the vector field.Note that the term in the denominator of the equation of motion for B(r) Eq. (15d) could in some cases lead to singularities.We analyze the behavior of the denominator by setting it equal to zero.This leads to a remarkable behavior when considering a quartic self-interaction potential V of the form where m is the mass of the vector boson and λ is the self-interaction parameter.We insert the potential Eq. ( 22) into the singular term in Eq. (15d) and obtain This expression holds for all radii.We analyze its behavior in the limit r → 0 by applying the initial conditions given in Eq. ( 19).One obtains a critical value for the central field amplitude E 0 : We here also defined the dimensionless interaction parameter Λ int = λ/8πm 2 .This expression constitutes an analytical upper bound for the central amplitude of the vector field.This means that any FPS with initial conditions for the field larger than E 0,crit will be physically forbidden, since Eq.(15d) will become singular and diverge.This result matches the analytical bound found by [84].
The relation implies that for strong self-interaction strengths Λ int , the allowed range for Proca stars becomes increasingly small and vanishes in the limit of very strong self-interactions.This fact could conceivably be used to constrain the vector field parameters m and λ.For example, a maximal vector field amplitude implies a maximal amount of accretion of vector bosons until the system becomes unstable.The field would then either dissipate to infinity, shed the excess vector field component, or collapse into a black hole.We leave a thorough investigation for future work.

D. Stability Criterion
Every FPS solution is characterized by the initial values for the central density ρ c and the central value of the vector field E 0 .When studying them in astrophysical contexts, the question of stability of FPSs naturally arises.The stability of pure Proca stars and NSs to radial perturbations is well known (see [75] for Proca stars).The stable and unstable solutions are separated by the point at which the total gravitational mass reaches its maximum with regard to the central density ρ c (for NS) and the central field E 0 (for Proca stars).Since FPSs are two-parameter solutions, the stability criterion needs to be modified.It was first presented for scalar FBSs by [93] (also see [30] for a review).But the criterion is more general and can also be applied to systems of two gravitationally interacting fluids.This is why we apply it here for FPSs.The idea behind the generalized stability criterion is to find extrema in the total number of particles (fermion number N f or boson number N b ) for a fixed total gravitational mass.The transition between stable and unstable configurations is given by the point at which where d/dσ denotes the derivative in the direction of constant total gravitational mass (see [93]).Up to a normalization factor, Eq. ( 25) can be written as If one is only interested in the precise points where FPSs become unstable, the unspecified normalization factor in Eq. ( 26) becomes irrelevant, since the whole relation is set to zero.In summary, the stability criterion Eq. ( 25) can be used to discriminate between astrophysically stable and unstable FPS solutions.When perturbed, unstable solutions will either collapse to a black hole, dissipate to infinity or migrate to a stable solution through internal re-configuration (see [30]).

E. Numerical Methods
In this work, we solve the equations Eq. (15a)-Eq.(15e) numerically to obtain self-consistent FPS solutions.We have implemented the algorithm in the code [88] which was developed by the authors of [15].The equations have one parameter undetermined by the boundary conditions Eq. ( 19), namely the vector field frequency ω.We use a shooting-algorithm to find ω numerically.
For given ρ c and E 0 , there exist only discrete values of ω, such that the boundary conditions at infinity Eq. ( 19) are fulfilled.These discrete values are called eigenvalues or modes.There are infinitely many of these modes.They are characterized by the number of roots (i.e., zero-crossings) the field E(r) has.Usually we are only interested in the lowest mode, since only it is believed to be dynamically stable [30].The lowest mode of the vector field always has one root in E(r).The following algorithm can however be used to find any desired mode.We integrate the system of ordinary differential equations Eq. (15a)-Eq.(15e) using a fifth order accurate Runge-Kutta-Fehlberg solver for some fixed value of ω.The vector field will then diverge towards positive or negative infinity at some finite radius.The system only converges at infinity if any mode is hit directly.But this is impossible to achieve numerically with finite precision.We thus make use of this diverging property to find the wanted frequency mode.When the frequency ω is close to the wanted mode, the divergence will happen at increasingly large radii, the closer the chosen value for ω is to the mode.A higher accuracy in finding ω will therefore push the divergence to larger radii.When ω is not exactly tuned to the mode, the vector field profile E(r) will diverge towards +∞ or −∞ and change its direction of divergence when ω passes a mode.The direction of divergence depends on which mode is solved for.For modes with an even number of roots, the field will diverge to +∞ if the frequency ω is below the mode, and it will diverge to −∞ if ω is above the mode.This will be reversed for all modes with an odd number of roots.By making use of the direction of divergence, we gain a binary criterion to find the correct mode.The value of ω can then be adapted -increased or decreased -based on the direction of divergence and the wanted mode.This procedure requires to integrate the system of equations multiple times with different values for ω, until the correct value is found.We implement this method in our code [88] using a bisection algorithm, which converges exponentially fast.We start with upper and lower values of ω, which are guaranteed to be smaller/larger than the wanted value of ω at the mode.In practice, lower and upper bounds of ω bound = [1,10] have proven to be numerically robust.We then perform the bisection search by taking the middle value of ω in this range and counting the number of roots in E(r) at each step.This also allows us to discriminate between different modes and to target specific modes by demanding a certain number of roots in the field E(r).The bisection is complete when the current value of ω found through bisection is close enough to the value of the mode.In our experience, the absolute accuracy needed to obtain robust solutions is on the order of ∆ω = |ω mode − ω bisection | ≈ 10 −15 .
Once a sufficiently accurate frequency ω is found, we modify the integration, such that E(r) and B(r) are set to zero at a finite radius r * B .This radius r * B is defined at the point where the field E(r) and its derivative E ′ (r) are small.This roughly corresponds to the last minimum of E(r) before it diverges.Also note that this is different to the bosonic radius R b defined previously.The condition can be summarized as the point where E(r * B )/E 0 < 10 −4 and E ′ (r * B ) ≪ 1.This is necessary because the interplay of the vector field and the NS matter can complicate the numerical solution.In some parts of the parameter space, especially for small initial densities ρ c , the vector field could diverge while still inside the NS component, i.e., before the pressure P (r) reaches zero (within numerical precision, we consider the pressure to be zero when P < 10 −15 ).This divergence would make finding physical values such as the fermionic radius R f impossible.Therefore, we artificially set E = B = 0 for r > r * B .This allows us to circumvent the divergence and accurately resolve the rest of the NS component.Note that the divergence of the vector field only happens because it is impossible to perfectly tune ω to the exact value, within numerical precision.If ω could be found exactly, the divergence would not happen.Setting the vector field to zero at some radius r * B is thus simply a way to maintain numerical stability of our algorithm.The condition was chosen so that the remaining contribution of the vector field to the other quantities (i.e., the metric components) is minimized.We have tested this method for different thresholds and confirmed that all extracted results are the same.
After integrating the solution to radii outside the matter sources (i.e.where E = B = P = 0), we can extract global observables such as the total gravitational mass and radius.The outside of the source is located at radii r larger than both the fermionic radius R f and r * B .In this regime, neither the NS matter nor the vector field contribute significantly.There, we can extract the total gravitational mass M tot Eq. ( 17) and then compute the integrals Eq. (16a) and Eq.(16b) to obtain the fermion/boson numbers N f , N b .The vector field convergence condition E(r * B )/E 0 < 10 −4 cannot be fulfilled for some configurations due to numerical precision limits.This generally happens for small initial field values E 0 ≲ 10 −4 , where the vector field extends far outside the NS component.In these cases, we extract the total gravitational mass M tot = 1 2 r ext (1 − a −2 (r ext )) at the point where its derivative has a global minimum.When the vector field diverges, also the metric components do, and with it also M tot .By taking the point where the derivative of the mass has a global minimum, which roughly corresponds to where the vector field and its derivative is closest to zero, we get the best possible estimate of the mass of the system before the divergence.
During our numerical analysis, we encountered the phenomenon that the bisection algorithm to find the frequency ω could fail for some specific initial conditions for E 0 and ρ c .We found this to be the case due to the bisection algorithm jumping over multiple modes in one iteration step.The wanted mode was then skipped and ended up outside the bisection bounds.The bisection then converged on an unwanted ω-value, or ended up failing entirely.We solved this problem by employing a backup algorithm that activates if the bisection fails.It restarts the bisection for ω but with different lower and upper bounds of ω bound .We tested the backup algorithm for 4800 FPS configurations with different vector field masses m and self-interaction strengths Λ int = λ/8πm 2 with equally distributed initial conditions for E 0 and ρ c .We found that 330 (≈ 6.8 %) of all configurations needed one restart of the bisection, and only 3 (≈ 0.06 %) of all configurations needed two restarts.In none of the tested cases, the bisection had to be restarted three times or more.

III. RESULTS
We consider FPSs with a quartic self-interaction potential of the same form as in Eq. (22).We further define the effective self-interaction parameter Λ int = λ/8πm 2 .The parameter Λ int is a useful measure for the self-interaction strength and parametrizes scaling relations for the total gravitational mass M max ≈ 1.058M 2 p /m [75] (for small Λ int ) and [84] (for large Λ int ).These scaling relations and numerical pre-factors can be derived numerically by fitting the configurations of maximum mass for a given m and Λ int to postulated scaling behaviors.Note that the parameter Λ int was originally introduced in the context of pure Proca stars and thus the scaling relations will not be generally valid for the mixed system.They can however be useful to understand the limiting cases where the FPS is dominated by the bosonic component.Nonetheless, we regard Λ int to be a useful measure to compare different choices of the mass and self-interaction strength.The self-interaction parameter Λ int in our work differs from the one used in [84] by a factor of two, even though they are defined in the same way.This is because a different normalization was used for the vector field.
We hereafter investigate models with parameters in the order of m ≈ 1.34•10 −10 eV and Λ int ≈ 0−100.This mass range is chosen because in this work we want to study DM which behaves as a macroscopic wave on typical neutron star length scales.Therefore, the correlation length of the ultralight DM particle is on the scale of km.Due to our units of c = G = M ⊙ = 1, lengths are measured in units of half the Schwarzschild radius of the Sun (≈ 1.48 km).Then the reduced Compton wavelength of the bosonic field is also measured in these units.m = 1 in our code units thus corresponds to 1.336 • 10 −10 eV (see a detailed explanation in Appendix A).The range for the self-interaction parameter was chosen so that it fulfills the observational constraints for the DM cross-section of 1 cm 2 /g obtained from the Bullet Cluster [94,95].The choice of Λ int is thus consistent with observations as long as it fulfills: .
Note that, in our conventions, the units of Λ int are [solar gravitational radius] divided by [1.336 • 10 −10 eV ].We used that the total cross-section of a self-interacting (scalar) particle is σ = λ 2 /16πm 2 .This should give a sufficient order of magnitude estimate for the cross-section for the vector particle, too.
For most calculations, we use the DD2 equation of state (with electrons) [96], taken from the CompOSE database [97], to describe the NS component.It was chosen because it is widely used by a number of groups and thus is well known in the literature.The DD2 EOS is based on a relativistic mean-field model with density-dependent coupling constants, which has been fitted to the properties of nuclei and results from Brueckner-Hartree-Fock calculations for dense nuclear matter.Therefore, the DD2 EOS describes also the EOS of pure neutron matter from chiral effective field theory (see [98]).For the purpose of our investigations, the particular choice of the nuclear equation of state is not of importance and has no effect on our general conclusions.

A. Radial Profiles
We compute the radial profiles of FPSs.In particular, we consider the radial dependence of the pressure P (r) and the vector field components E(r), B(r).Even though the radial distribution of physical quantities can not yet be observed directly (although one could infer the DM distribution using the geodesic motion of light [99]), a good understanding of the internal structure of FPSs can be used to deduce their global quantities and vice versa.Knowledge about the internal distribution is also relevant for numerical applications.Another reason we include the radial profiles here is to facilitate reproducibility of this work and for the sake of code-validation for future works.Radial profiles of pure Proca stars have already been discussed by [75] and for the case of a quartic self-interaction potential like Eq. ( 22) by [84].We used the results of [84] in particular to verify that our code [88] reproduces the results correctly and consistently.
In Figure 1, we show radial profiles of the pressure P (r) (orange) and the vector field components E(r) (black), B(r) (blue) of the zeroth mode of different FPSs with potential Eq. ( 22).In the left panel, we take a boson mass of m = 1.34 • 10 −10 eV and an interaction strength of Λ int = 0.The FPSs have varying central vector field amplitudes E 0 and central densities of ρ c = 5ρ sat .Here we take ρ sat = m n n nuc ≈ 2.5 • 10 14 g/cm 3 to be the nuclear saturation density, with m n being the neutron mass and n nuc = 0.15 fm −3 the average nuclear number density.The radial profile of a pure NS is A kink can be seen in the profile for B(r) at roughly 11.5 km.This corresponds to the point where the fermionic radius is located.This illustrates the gravitational back-reaction between the vector field and NS matter.24), the maximal amplitude is roughly E0,crit ≈ 0.0282.The limited field amplitude strongly limits the effect on the fermionic component.
shown with the orange continuous line and has no corresponding vector field (because it would be zero everywhere).The presence of the DM can be seen to compactify the NS component with increasing central field amplitude E 0 .The field forms a DM core configuration.
In the right panel of Figure 1, all parameters are left equal except for the vector boson mass, which is set to m = 1.34 • 10 −11 eV .Due to the low DM mass, the correlation length increases, which increases the size of the vector field component and forms a DM cloud configuration.Since the amount of energy density of the vector field is distributed inside and outside the NS component, the effect on the radius is small.At around r = 11.5 km, a kink can be seen in the radial profile of the field component B(r).This point coincides with the point where the fermionic radius of the FBS is located.This illustrates the gravitational back-reaction between the vector field and the NS component of the FBS.
In Figure 2, we show radial profiles of the pressure P (r) (orange) and the vector field components E(r) (black), B(r) (blue) of an FPS.In the left panel, we show an FPS in the first mode, which can be identified by the fact that the E(r) component crosses the x-axis twice and B(r) crosses it once.The boson mass is m = 1.005 • 10 −10 eV and Λ int = 0.This time, the central density is taken to be ρ c = 4ρ sat and the central vector field amplitudes vary.The right panel of Figure 2 shows an FPS in the zeroth mode with a vector boson mass of m = 1.34 • 10 −10 eV and a self-interaction strength of Λ int = 50.The maximal amplitude is roughly E 0,crit ≈ 0.0282 due to the analytical bound on E 0 , see Eq. ( 24).The limited field amplitude strongly limits the possible effect on the fermionic component and thus on the fermionic radius, especially in the limit of large Λ int .It may therefore be difficult to detect strongly self-interacting vector DM within a NS if one only considers measurements of the fermionic radius.It is also conceivable that the maximum amplitude E 0,crit implies a maximum amount of possible accretion of vector DM, which could be used to set bounds on the DM self-interaction strength.We leave the analysis of this aspect for a future work.

B. Stable Solutions
We compute a grid of FPSs with different central densities ρ c and central vector field amplitudes E 0 .Using the array of solutions, we compute the stability curve using the stability criterion Eq. ( 26).The stable solutions can then be filtered and analyzed further.This can be seen in the left panel of Figure 3, where we compute FPSs with a quartic self-interaction potential Eq. ( 22) with m = 1.34 • 10 −10 eV and Λ int = 0. We additionally compute the stability curve using the stability criterion Eq. ( 26).The stability curve defines the boundary between stable and unstable configurations under linear radial perturbations.The shape of the stability curve for FPSs is qualitatively very similar to the case of scalar FBSs (compare to [15]).For pure neutron stars and Proca stars, respectively, the curve converges on the ρ c -and E 0 -axis at the point, where the non-mixed configurations have their maximum gravitational masses.We take only the FPSs inside the stability region, enclosed by the stability curve, and plot them in a mass-radius (MR) diagram.This leads to the graph in the right panel of Figure 3.We see that stable FPS configurations form an MR region instead of an MR curve (which would be the case for single-fluid systems).The stable configurations form core or cloud solutions, depending on their DM-fraction N b /(N b + N f ).A DM core is present if the bosonic component is geometrically smaller than the fermionic component (i.e. when R b < R f ).The opposite is true for cloud configurations.The FPSs with high DM-fractions have masses of roughly 1 M ⊙ .This is higher than for scalar FBSs with equal boson mass m (compare to [15]).This can be explained through the different scaling relations for pure Proca stars and boson stars.
Another point where FPSs differ from FBSs is the existence of a maximal amplitude E 0,crit Eq. ( 24) for the vector field.When increasing the self-interaction strength Λ int , the maximal possible vector field amplitude shrinks.This affects the shape of the stability curve.In Figure 4 (left panel), we show such a case where the self-interaction strength is Λ int = 5.The stability curve does not reach the E 0 -axis anymore, but instead rises vertically from the pure NS configurations until it reaches the FPSs with maximal central vector field amplitude E 0,crit ≈ 0.089.We have manually extended the stability curve so that it proceeds horizontally until it reaches the E 0 -axis.It is noteworthy that this behavior starts at surprisingly small self-interaction strengths and persists up to higher Λ int .In principle, also a third behavior of the stability curve of FPSs is conceivable.For some specific Λ int , it should be possible that the stability curve does not admit one continuous shape like in Figure 3 or Figure 4, but that the stability curve is cut into two parts.Namely, one part which starts at the E 0 -axis and then rises to reach the edge where E 0,crit is located, and another part which starts at the ρ c -axis and then rises roughly vertically until it too reaches the analytical bound for the vector field amplitude E 0,crit (think of a horizontal line cutting through the stability curve in Figure 3 at, e.g., E 0 = 0.06).During our testing, we did not find any case where the stability curve follows this behavior.However, there is also no reason that we are aware of why such a behavior of the stability curve should be forbidden.This is why we presume that such a case might exist.We compute various FPSs with different values of the vector boson masses m = {1, 10, 0.1} × 1.34 • 10 −10 eV and self-interaction strengths Λ int = {0, 10, 100}.We chose the same parameter values as in [15] to allow for easy comparability.In Figure 5, we show the mass and fermionic radii of all stable FPS configurations in an MR diagram.In Figure 6, we show the mass plotted against the effective gravitational radius R G .It is defined as the radius where The gray region marks the Buchdahl limit, where no stable NS can exist.Observing only R f of these systems would appear to violate the Buchdahl limit, even though the FPS as a whole does not.RG is the radius where 99% of the total rest mass is contained.The rows correspond to bosonic masses of m = {1, 10, 0.1} × 1.34 • 10 −10 eV , columns correspond to self-interactions of Λint = {0, 10, 100} respectively.We use the DD2 EOS for the fermionic part.Notice the different scales of the bottom plots.For pure NSs, because the crust has comparatively low density, RG is significantly smaller than R f (compare to Figure 5).RG tends to be higher as compared to scalar FBSs for equal boson masses and self-interaction strength (compare to Figure 3 in [15]).
99 % of the total rest mass N f +N b is contained.The stable solutions have been obtained using the stability criterion Eq. ( 26).Note the different axis scaling in the figures.It was chosen such that the relevant trends and features of the solutions can be seen well.We hereafter discuss some general trends and compare the results to the one obtained for scalar FBSs.
The following analysis should thus be explicitly compared to figures 2 and 3 in [15].We find that many of the general conclusions regarding FBSs can also be applied to FPSs.FPSs with small DM-fractions are dominated by the fermionic component, leading to only small changes in the fermionic radius.In the case of DM dominated FPSs, the solutions behave similar to pure Proca stars.This leads to higher masses as compared to FBSs, where the total gravitational mass of pure boson stars will be roughly half that of a Proca star with the same boson mass, as can be seen well for the cases where m = {1, 0.1} × 1.34 • 10 −10 eV .FPSs can thus reach higher total gravitational masses as compared to FBSs with the same DM mass and selfinteraction strength.For m = 1.34 • 10 −9 eV , the bosonic component is concentrated inside the fermionic one and forms a DM core.Even small amounts of DM can have a significant impact on the fermionic radius, since the whole vector field is concentrated entirely inside the NS component.More massive DM particles can thus have larger effects on the fermionic radius compared to low-mass DM at similar DM-fractions.This is due to the cloud-like structure of low-mass DM.For small DM masses, the majority of the DM will be concentrated outside the NS part -due to its larger correlation length -and will thus have smaller effects on the fermionic radius.The smaller the mass and the larger the self-interaction strength, the more likely the formation of a DM cloud is.The opposite is true for DM core solutions.FPSs tend to produce configurations with larger total masses compared to scalar FBSs.Their halos also extend to larger radii, as can be seen from the gravitational radius in Figure 6.
In general, the gravitational radius of FPSs is larger in size as compared to scalar FBSs (compare to Figure 3 in [15]).The larger gravitational radius suggests that FPS have larger tidal deformabilities, compared to their scalar field counterparts (FBS) with equal m and Λ int .This is because objects with larger radii are generally favored to tidally disrupt.This could favor higher vector boson masses compared to the corresponding scalar boson mass in the case of FBSs.A future quantitative analysis of the tidal deformability of FPSs is needed to definitively verify this hypothesis.When considering the gravitational radius of FPSs with small boson masses of m = 1.34•10 −11 eV (bottom row of Figure 6), the transition between DMdominated and NS-dominated configurations appears more abrupt than in the FBS case (compare to Figure 3 in [15]).For example, when starting with a system with a DM-fraction of roughly 0% or 80%, increasing the DM-fraction by small amounts can massively impact the total mass and gravitational radius of the combined system.Finally, note the outlier points in Figure 6 for m = 1.34 • 10 −11 eV and Λ int = 100 at roughly R G = 350 km.These are likely to be numerical artifacts and should thus not be regarded as physical.This is to be expected since for small DM masses and large self-interactions, the numerical solution gets increasingly difficult.This problem could be avoided by using smaller step-sizes and higher numerical precision.But this would also lead to longer run-times of the code.

C. Comparison with Scalar FBS
We show MR relations of FPSs and scalar FBSs with fixed DM-fractions In the left panel of Figure 7, we show different FPSs with constant DM-fractions.The DD2 EOS [96] was used for the NS component.For the vector boson, we chose masses of m = {1, 0.1} × 1.34 • 10 −10 eV and no self-interactions.This figure should be explicitly compared to Figure 5 (left panel) in [15] as the same masses and DM-fractions were chosen.The MR curve of a pure NS with the DD2 EOS (black line) is shown as a reference.Depending on the boson mass, FPSs can have increased or decreased maximum total gravitational mass when there is vector DM present.FPSs tend to produce configurations with larger gravitational masses compared to FBSs with equal parameters (mass, self-interaction and DM-fraction).This is not surprising when considering the scaling relations of pure boson stars and Proca stars, respectively.The gravitational mass scales like M max ≈ 0.633M 2 p /m for pure boson stars and like M max ≈ 1.058M 2 p /m for pure Proca stars, where m is the mass of the scalar/vector boson, respectively.The presence of light bosonic DM can help to increase the total gravitational mass of a NS.This can make EOS which do not fulfill the observational constraints for the maximum NS mass viable again.Vector DM has a larger effect on the gravitational mass than scalar DM and thus smaller amounts of vector DM are needed to produce an equal increase in the total gravitational mass.In the right panel of Figure 7, we show different FPSs (orange and green lines) and FBSs (blue lines)  ).The solid black-white line shows the mass-radius curve for pure fermionic matter.The vector boson mass was chosen so that in the limit of pure boson stars/Proca stars, the same total gravitational mass is produced.Both diagrams show only marginal differences.
for different boson masses, no self-interactions and constant DM-fractions N b /(N b + N f ).We used the DD2 EOS [96] for the NS component.The parameters were chosen in a way to illustrate the degeneracies that can arise from different DM models or EOS for the NS component.For example, FPSs and FBSs with boson masses of m = 1.34 • 10 −11 eV (dashed lines) produce virtually indistinguishable mass-radius relations, when the FPSs and the FBSs have a DM-fraction of 60% and 75% respectively.
A similar behavior can be seen for the cases where the boson mass is m = 1.34 • 10 −10 eV (dot-dashed lines).Here, FBSs with 15% DM-fraction produce similar MR curves to FPSs with 20% DM-fraction.In addition, the resulting MR curves are comparable to the curve corresponding to a pure NS with the KDE0v1 EOS [100].They also match the curve corresponding to an FPS with 10% DM-fraction and a vector boson mass of m = 2.24 • 10 −10 eV (green line).In conclusion, FPSs can produce degenerate results in the MR plane with both FBSs and pure NS, given that different DM-fractions and EOS are allowed.Additional observables, such as the tidal deformability, are needed to break the degeneracy.However, it seems difficult to prevent degenerate solutions from existing in general, since FPSs themselves can be degenerate with other FPS-solutions that have different boson masses and DM-fractions.
We further explore the degeneracy between FPS and FBS solutions.In Figure 8, we show the stable FBS and FPS solutions in an MR diagram.We used the scaling relations of the maximum mass for pure boson stars (M max ≈ 0.633M 2 p /m) and pure Proca stars (M max ≈ 1.058M 2 p /m) to match the boson masses in a way that both FPSs and FBSs will have the same gravitational mass in the pure boson star/Proca star limit.To guarantee matching solutions in this limit, we chose a scalar boson mass of m = 1.34 • 10 −10 eV and we chose a mass of 1.058 ÷ 0.633 ≈ 1.671 times the mass of the scalar boson -i.e.m = 2.24 • 10 −10 eV -for the vector boson.We find a high degree of similarity between the MR region of FBSs and FPSs with the scaled masses.This makes both solutions almost indistinguishable.The small differences present between the left and right panel of Figure 8 can be attributed to a slightly different grid-spacing used for the initial conditions ρ c , ϕ c (and ρ c , E 0 ).This can be seen in the MR regions at small total gravitational masses M tot < 0.5 M ⊙ and also at radii R f > 15 km.The color shading further reveals a different distribution of DM-fractions for a given M tot and R f , even though the difference is small.We expect a similar behavior to hold when considering different scalar and vector boson masses (with zero-self-interaction), given that they differ by the same factor of ≈ 1.671.This adds further confidence to the observation that FBSs and FPSs might be difficult to distinguish since a given solution might be another system but with different boson mass (or DM-fraction).Similar scaling relations also exist for boson stars and Proca stars in the limit of large self-interactions Λ int .A similar procedure might therefore be possi-ble when also matching the self-interaction strength appropriately.An independent measurement of the DM particle mass would break this degeneracy to a certain degree.But it would also be necessary to constrain the self-interaction strength and the DMfraction through other means.For example using correlations of the DM abundance in the galactic disk (see [101,102]) or using the bound on the maximal vector field amplitude.The scaling behavior between (fermion) boson stars and (fermion) Proca stars also suggests another application.If it persists for large non-zero selfinteractions, it might be possible to use the effective bosonic EOS derived by Colpi et al. [103] also to model (fermion) Proca stars.Since the EOS by Colpi et al. was originally derived for a scalar field, one would then have to scale the boson mass by a factor of 1.671 and the self-interaction by an appropriate amount.The necessary scaling for the self-interaction will be dictated by the scaling relations for pure boson stars (M max ≈ 0.22 ) at large self-interaction strengths.We however note that great care is needed since Proca stars technically do not exist in the limit of large self-interactions (see the analytical bound on the vector field amplitude Eq. ( 24)).We plan to study this aspect in the future.

D. Higher Modes and Different EOS
We broaden our analysis to FPSs with different EOS for the fermionic component and to FPSs where the bosonic component exists in a higher mode.Higher modes are usually assumed to be unstable, but as numerical simulations of scalar boson stars have shown [37,104], higher modes might be dynamically stable when gravitationally interacting in a multi-component system.We therefore start by considering FPSs in the first and second mode in Figure 9.In the left panel of Figure 9, we show the total gravitational mass and the fermionic radius of stable FPS configurations, where the bosonic component is in the first mode (as opposed to the ground mode, which is the zeroth mode).The vector boson mass is m = 1.34 • 10 −10 eV , and the self-interaction is set to zero.We first note the fact, that stable solutions under linear radial perturbations, according to the stability criterion Eq. ( 26), exist at all.This is a non-trivial statement as higher modes of Proca stars (and also of scalar boson stars) are usually believed to be unstable.Note however that the stability criterion Eq. ( 26) is merely a necessary condition for stability and that there could be additional conditions that must be fulfilled for a solution to be stable in higher modes.Also, our stability analysis does not consider the dynamical stability of the higher modes.They might thus be unstable in non-static scenarios.It is however possible that the higher modes of the bosonic part might be stabilized through the gravitational interaction with the fermionic part of the FPS.In general, solutions in higher modes need energy to be supported or excited.Otherwise they would decay to the ground mode.The authors of [104] explicitly studied the stability of higher modes for scalar boson stars with multiple scalar fields.They investigated cases where one field is in the ground mode and another is in the excited state.They found that the excited mode is stable if the charge of the conserved Noether current (which can also be interpreted as the particle number or the total restmass, see Eq. ( 11)) of the ground mode is larger than the Noether charge of the excited mode: N ground > N excited .If this logic also holds for mixed systems of fermions and bosons, this would imply an additional stability condition that the fermion number N f must be larger than the boson number N b of the FPS in the excited mode.We further refer to chapter 3.7 in [30] for a more detailed review.The FPSs in the first mode exhibit higher gravitational masses in the configurations dominated by the bosonic component, compared to FPSs in the zeroth mode (compare to Figure 3).The numerical value of the frequency ω in the higher mode is also larger than the frequency in lower modes.This behavior is con-sistent with earlier works, which studied pure Proca stars analytically [84] and numerically [79].They also observed that higher frequencies lead to larger total gravitational mass.The left panel of Figure 9 shows a number of outlier points at around 11 km and 2.3 M ⊙ .These are likely numerical artifacts due to the increased difficulty of finding accurate numerical solutions for higher modes.The right panel of Figure 9 shows stable FPS configurations in the second mode.The vector boson mass is m = 1.34 • 10 −10 eV and the self-interaction is set to zero.Here also, the existence of stable solutions is to be acknowledged.In the limit of high DM-fractions, the FPSs converge to the solution of pure Proca stars and reach total gravitational masses of roughly 2.5 times that of Proca stars in the zeroth mode (compare to Figure 3).In comparison to the case in the first mode (left panel of Figure 9), the quality of the overall solution can be seen to deteriorate further.We believe the outlier points at roughly < 13 km and 1 M ⊙ to be non-physical numerical artifacts.The outlier points coincide with the solutions in the zeroth mode.This suggests that our solver did not find the second mode in these cases and converged on the zeroth mode instead.Solutions of FPSs in even higher modes should therefore be considered with great care.The difficulty of obtaining accurate numerical solutions is likely to increase further for higher modes.The quality of the solution is however sufficient to gain a qualitative understanding of FPSs in higher modes.In conclusion, higher modes are stable under linear radial perturbations and increase the total gravitational mass of FPSs by substantial amounts.
We investigate the effect that different EOS have on FPSs.In Figure 10, we use the APR EOS [105] for the fermionic part.We chose a vector boson mass of m = 1.34•10 −10 eV with no self-interaction for the bosonic part.In the left panel, we notice that the shape of the stability curve (black curve) is affected by the choice of the EOS.On the ρ c -axis, it converges to a value of around 7.5ρ sat .This is higher than the corresponding value of ρ c when the DD2 EOS is used (compare to Figure 3) because the APR EOS is softer than the DD2 EOS.This means that the nuclear matter is easier to compress and higher central densities can be supported by the EOS.The easier compressibility also shows itself through smaller NS radii (see the right panel).In the limit of pure Proca stars, the stability curve converges to the same value as it does when the DD2 EOS is used (compare to Figure 3).The MR region shows a similar qualitative behavior as in the DD2 case.The high DM-fraction limit in particular shows a convergence to the solution to pure Proca stars.The APR EOS also allows higher central amplitudes of the vector field E 0 , compared to the DD2 EOS with equal boson mass and self-interaction strength.Figure 11 shows different FPS configurations where the FSG EOS [96] was used for the fermionic part.For the bosonic part, we used a boson mass of m = 3.01•10 −11 eV and no self-interaction.The FSG EOS is a soft EOS and thus reaches higher central densities ρ c for pure NSs.It is excluded by current observational constraints (see Figure 7), as it cannot produce pure NSs with masses of M = 2.35 +0.17 −0.17 M ⊙ [8].However, adding DM to the pure NSs can significantly increase the maximum gravitational mass of the combined system.The FSG EOS is then able to reach the observational bound on the maximum NS mass in the presence of DM.In fact, the MR curve of the pure DD2 EOS is entirely contained within the stability region of the FPSs with the FSG EOS.This again raises the point that some FPS solutions are degenerate with some NS solutions (see Figure 8), when allowing for different DM-fraction and DM masses.Another factor complicating the identification of the EOS vs. the DM effects is that the presence of dark matter might change the fermionic radius produced by a given EOS.For example, [92] found that the presence of self-interacting and repulsive fermionic dark matter can lead to nearly indistinguishable fermionic radii for different EOS.To figure out whether and which types of mixed DM-NS systems might exist, it will be crucial to perform sophisticated parameter searches of the system and obtain more measurements to constrain the DM and NS properties in future studies.

IV. CONCLUSIONS
In this work, we studied the impact that bosonic dark matter (DM) has on the mass and radius of neutron stars (NSs).DM was modeled as a massive, self-interacting complex vector field.DM was further assumed to only interact gravitationally with the fermionic neutron star matter.We derived the equations of motion describing static spherically symmetric fermion Proca stars (FPSs) and computed their properties numerically.We also found a scaling relation between the frequency, vector field and metric components, and we derived an analytical upper bound on the vector field amplitude.
We showed that the presence of the vector field can lead to core-like and to cloud-like solutions.Core-like solutions can increase the compactness of the NS component.For some configurations, observing only the fermionic radius and the total gravitational mass would appear to violate the Buchdahl limit.We found core-like solutions for vector boson masses of m ≳ 1.34 • 10 −10 eV and small self-interactions Λ int = λ/8πm 2 .Cloud-like solutions appeared when m ≲ 1.34 • 10 −11 eV and Λ int is large.For some small boson masses m ≲ 1.34 • 10 −11 eV , the presence of DM can significantly increase the total gravitational mass while leaving the fermionic radius approximately constant.We computed radial profiles of FPSs and found that the existence of a maximum possible vector field amplitude limits the effect of DM on the NS when the self-interaction Λ int is large.The maximum amplitude implies a maximum possible amount of vector boson DM accretion and could thus be used to set bounds on the DM properties.
We also compared FPSs to FBSs with a scalar field.We used the same parameters as in [15] to simplify the comparison.For stable FPS configurations, we found that many of the general qualitative trends that apply to FBSs also apply to FPSs.But vector DM leads to higher FPS masses and larger gravitational radii for equal m and Λ int .This could also imply a larger tidal deformability of FPSs compared to FBSs.Also, a measurement of the gravitational radius would favor larger vector boson masses compared to scalar boson masses.For FPS configurations of constant DM-fraction, we found that the effect of vector DM on the NS properties (total gravitational mass and fermionic Additionally we show contours of constant gravitational mass.The black line corresponds to the stability curve, which separates stable solutions (in the lower left region) from unstable solutions (everywhere else).The qualitative behavior of the stability curve of is similar to the case with the DD2 EOS (see Figure 3) Right panel: Mass-radius diagram displaying the fermionic radius vs. the total gravitational mass for FPS configurations that are within the stability region displayed in the left panel.Each point corresponds to a single configuration and is color-coded according to the DM-fraction N b /(N b + N f ).The solid black-white line shows the mass-radius curve for pure fermionic matter.In both cases, a vector field with a mass of m = 1.34 • 10 −10 eV and no self-interactions was considered in addition to the APR EOS [105] for the fermionic part.).The solid black-white line shows the mass-radius curve for pure fermionic matter.In both cases, a vector field with a mass of m = 3.01 • 10 −11 eV and no self-interactions was considered in addition to the FSG EOS [96] for the fermionic part.
radius) is larger compared to FBSs with equal DM-fraction, mass m and self-interaction strength Λ int .One therefore needs a larger amount of scalar DM to cause the same effect as vector DM.For different boson masses and DM-fractions, we found that FPSs and FBSs can both be degenerate with each other and also be degenerate with pure NS with a different EOS.We found an especially high degree of similarity between FBS solutions with no self-interaction and a boson mass of m = 1.34 • 10 −11 eV with FPS solutions where the vector boson mass is larger by a factor of 1.671.We expect the similarity in the behavior to hold also for different boson masses (and also for non-zero self-interactions), as long as the vector boson mass is scaled accordingly by the right factor.These similarities also hint towards a possibility to use the effective EOS by Colpi et al. [103] also for (fermion) Proca stars.We however note that great care is needed since Proca stars do not exist in the limit of large self-interactions (see the analytical bound on the vector field amplitude Eq. ( 24)).The similarities between FBSs and FPSs might also be useful for numerical applications.Scalar (fermion) boson stars are easier to implement and numerically cheaper to solve than FPSs.One could then simply solve the equations for scalar (fermion) boson stars with a re-scaled mass (and self-interaction parameter Λ int ) to compute the properties (M tot , R f ) of (fermion) Proca stars.The prevalence of degenerate solutions highlights the importance of measuring additional observables, such as the tidal deformability, to break the degeneracies.
We confirmed the existence of higher modes that are stable under first-order radial perturbations.We found that higher modes lead to higher total gravitational masses of the mixed FPS systems.Using FPSs with different EOS for the fermionic part, we explicitly confirmed that for certain DM masses, previously excluded EOS are able to fulfill observational bounds if DM is present.Mixed systems of bosonic DM and NS matter can therefore be consistent with all current observational constraints if suitable boson masses and self-interaction strengths are chosen.

1 FIG. 1 .
FIG.1.Left panel: Radial profiles of the pressure P (r) (orange) and the vector field components E(r) (black), B(r) (blue) of the zeroth mode of different FPSs with potential Eq.(22).The boson mass is m = 1.34 • 10 −10 eV and Λint = 0.The FPSs have a central density of ρc = 5ρsat and varying central vector field amplitudes E0.The pressure has been re-scaled by a factor of 3 for convenience.The DM forms a core and compactifies the fermionic component.Right panel: Same as in the left panel, but this time the vector boson mass is set to m = 1.34 • 10 −11 eV .The DM forms a cloud around the fermionic component.The radius of the fermionic component is barely affected by the field.A kink can be seen in the profile for B(r) at roughly 11.5 km.This corresponds to the point where the fermionic radius is located.This illustrates the gravitational back-reaction between the vector field and NS matter.

FIG. 2 .
FIG. 2. Left panel:Radial profiles of the pressure P (r) (orange) and the vector field components E(r) (black), B(r) (blue) of the first mode of different FPSs with potential Eq.(22).The boson mass is m = 1.005 • 10 −10 eV and Λint = 0.The FPSs have a central density of ρc = 4ρsat and varying central vector field amplitudes E0.The pressure has been re-scaled by a factor of 3 for convenience.Right panel: Radial profiles of the pressure P (r) (orange) and the vector field components E(r) (black), B(r) (blue) of FPSs in the zeroth mode with potential Eq. (22).The boson mass is m = 1.34 • 10 −10 eV and the self-interaction strength is Λint = 50.The FPSs have a central density of ρc = 5ρsat and varying central vector field amplitudes E0.The pressure has been re-scaled by a factor of 3 for convenience.Due to the analytical bound on E0 Eq. (24), the maximal amplitude is roughly E0,crit ≈ 0.0282.The limited field amplitude strongly limits the effect on the fermionic component.

FIG. 3 .
FIG.3.Left panel: Total gravitational mass of different FPSs as a function of the restmass density ρc and central vector field amplitude E0.Additionally we show contours of constant gravitational mass.The black line corresponds to the stability curve, which separates stable solutions (in the lower left region) from unstable solutions (everywhere else).Right panel: Mass-radius diagram displaying the fermionic radius vs the total gravitational mass for FPS configurations that are within the stability region displayed in the left panel.Each point corresponds to a single configuration and is colour-coded according to the DM-fraction N b /(N b + N f ).The solid black-white line shows the mass-radius curve for pure fermionic matter.In both cases, a vector field with mass of m = 1.34 • 10 −10 eV and no self-interactions was considered in addition to the DD2 EOS for the fermionic part.

FIG. 4 .
FIG. 4. Left panel: Total gravitational mass of different FPSs as a function of the restmass density ρc and central vector field amplitude E0, with m = 1.34 • 10 −10 eV and Λint = 5.Additionally we show contours of constant gravitational mass.The black line corresponds to the stability curve, which separates stable solutions (in the lower left region) from unstable solutions (everywhere else).The stability curve reaches configurations with the maximum possible vector field amplitude E0,crit ≈ 0.089.This is a feature unique to FPSs.Right panel: Mass-radius diagram displaying the fermionic radius vs the total gravitational mass for FPS configurations that are within the stability region displayed in the left panel.Each point corresponds to a single configuration and is colour-coded according to the DM-fraction N b /(N b + N f ).The solid black-white line shows the mass-radius curve for pure fermionic matter.A vector field with mass of m = 1.34 • 10 −10 eV and Λint = 5 was considered in addition to the DD2 EOS for the fermionic part.

FIG. 5 .
FIG. 5. Relation between total gravitational mass Mtot and fermionic radius R f for different FPSs.The rows correspond to bosonic masses of m = {1, 10, 0.1} × 1.34 • 10 −10 eV , columns correspond to self-interactions of Λint = {0, 10, 100} respectively.We use the DD2 EOS for the fermionic part.Notice the different scale of the bottom plots.The gray region marks the Buchdahl limit, where no stable NS can exist.Observing only R f of these systems would appear to violate the Buchdahl limit, even though the FPS as a whole does not.

mFIG. 6 .
FIG.6.Relation between total gravitational mass Mtot and effective gravitational radius RG for different FPSs.RG is the radius where 99% of the total rest mass is contained.The rows correspond to bosonic masses of m = {1, 10, 0.1} × 1.34 • 10 −10 eV , columns correspond to self-interactions of Λint = {0, 10, 100} respectively.We use the DD2 EOS for the fermionic part.Notice the different scales of the bottom plots.For pure NSs, because the crust has comparatively low density, RG is significantly smaller than R f (compare to Figure5).RG tends to be higher as compared to scalar FBSs for equal boson masses and self-interaction strength (compare to Figure3in[15]).
FIG. 7. Left panel: Mass-radius relations of FPSs with the DD2 EOS [96] for vector boson masses m = {1, 0.1} × 1.34•10 −10 eV , no self-interactions and constant DM-fractions N b /(N b +N f ).This figure should be compared to Figure 5 (left panel) in[15] as the same masses and DM-fractions were chosen.The orange band marks the observational constraint of J0952-0607[8] and the percentage numbers denote the respective DM-fractions.Right panel: Massradius relations of FPSs (orange and green lines) and FBSs (blue lines) with the DD2 EOS for different boson masses, no self-interactions and different DM-fractions.The black lines correspond to the pure NSs with the DD2 EOS and KDE0v1 EOS[100] respectively.FPS and FBS solutions with different masses and DM-fractions can both be degenerate with each other, or also degenerate with pure NSs with a different EOS.

FIG. 8 .
FIG. 8. Left panel: Mass-radius diagram displaying the fermionic radius vs. the total gravitational mass for stable FBS configurations with scalar boson mass of m = 1.34 • 10 −10 eV and no self-interaction.Each point corresponds to a single configuration and is color-coded according to the DM-fraction N b /(N b + N f ).The solid black-white line shows the mass-radius curve for pure fermionic matter, modeled by the DD2 EOS.Right panel: Mass-radius diagram displaying the fermionic radius vs. the total gravitational mass for stable FPS configurations with vector boson mass of m = 2.24 • 10 −10 eV and no self-interaction.Each point corresponds to a single configuration and is color-coded according to the DM-fraction N b /(N b + N f).The solid black-white line shows the mass-radius curve for pure fermionic matter.The vector boson mass was chosen so that in the limit of pure boson stars/Proca stars, the same total gravitational mass is produced.Both diagrams show only marginal differences.

FIG. 9 .
FIG. 9. Left panel: Mass-radius diagram displaying the fermionic radius vs. the total gravitational mass for stable FPS configurations in the first mode with vector boson mass of m = 1.34 • 10 −10 eV and no self-interaction.Each point corresponds to a single configuration and is color-coded according to the DM-fraction N b /(N b + N f ).The solid black-white line shows the mass-radius curve for pure fermionic matter.Right panel: Mass-radius diagram displaying the fermionic radius vs. the total gravitational mass for stable FPS configurations in the second mode with vector boson mass of m = 1.34 • 10 −10 eV and no self-interaction.Each point corresponds to a single configuration and is color-coded according to the DM-fraction N b /(N b + N f ).The solid black-white line shows the mass-radius curve for pure fermionic matter.

FIG. 10 .
FIG.10.Left panel: Total gravitational mass of different FPSs as a function of the rest mass density ρc and central vector field amplitude E0.Additionally we show contours of constant gravitational mass.The black line corresponds to the stability curve, which separates stable solutions (in the lower left region) from unstable solutions (everywhere else).The qualitative behavior of the stability curve of is similar to the case with the DD2 EOS (see Figure3) Right panel: Mass-radius diagram displaying the fermionic radius vs. the total gravitational mass for FPS configurations that are within the stability region displayed in the left panel.Each point corresponds to a single configuration and is color-coded according to the DM-fraction N b /(N b + N f ).The solid black-white line shows the mass-radius curve for pure fermionic matter.In both cases, a vector field with a mass of m = 1.34 • 10 −10 eV and no self-interactions was considered in addition to the APR EOS[105] for the fermionic part.

FIG. 11 .
FIG.11.Left panel: Total gravitational mass of different FPSs as a function of the rest mass density ρc and central vector field amplitude E0.Additionally we show contours of constant gravitational mass.The black line corresponds to the stability curve, which separates stable solutions (in the lower left region) from unstable solutions (everywhere else).Right panel: Mass-radius diagram displaying the fermionic radius vs. the total gravitational mass for FPS configurations that are within the stability region displayed in the left panel.Each point corresponds to a single configuration and is color-coded according to the DM-fraction N b /(N b + N f ).The solid black-white line shows the mass-radius curve for pure fermionic matter.In both cases, a vector field with a mass of m = 3.01 • 10 −11 eV and no self-interactions was considered in addition to the FSG EOS[96] for the fermionic part.