Fifth-Force Screening around Extremely Compact Sources

Many non-linear scalar field theories possess a screening mechanism that can suppress any associated fifth force in dense environments. As a result, these theories can evade local experimental tests of new forces. Chameleon-like screening, which occurs because of non-linearities in the scalar potential or the coupling to matter, is well understood around extended objects. However, many experimental tests of these theories involve objects with spatial extent much smaller than the scalar field's Compton wavelength, and which could therefore be considered point-like. In this work, we determine how the fifth forces are screened in the limit that the source objects become extremely compact.


I. INTRODUCTION
Non-linear scalar field theories are popular candidates to help explain many cosmological mysteries, including dark matter [1,2], dark energy [3] and inflation [4]. They also appear in many modified theories of gravity, which often introduce new scalar degrees of freedom [3] -so called scalar-tensor theories. Of course, the Standard Model Higgs sector [5][6][7] is also described by a non-linear scalar field theory, although the energy scales involved are significantly higher than those of interest to late-universe cosmology.
In the absence of a symmetry principle [8][9][10][11], it is difficult to forbid a new light scalar from coupling to matter, and, as a result, we expect it to mediate a long-range fifth force. A canonical scalar field with a linear equation of motion has its mass and matter coupling fixed, and thus the existence of such forces is difficult to reconcile with terrestrial and Solar System tests of gravity [12]. However, the phenomenology of non-linear scalar field theories is much broader and allows for more interesting behaviour. If the scalar theory is non-linear then the behaviour of fluctuations of the field, which mediate the fifth force, will depend on the properties of the background scalar evolution. In this way, the properties of the force, such as its coupling strength and range, can vary with the environment. Of particular interest are theories in which the fifth force is dynamically suppressed in and around dense and compact sources. This phenomenon is known as screening, and such screening mechanisms can explain why the effects of new forces might evade laboratory experiments or in precision Solar System tests of gravity.
In this work, we focus on theories that screen through non-linearities in the field equations for the scalar φ. The non-linearities may be terms in the scalar potential, in the way the field couples to matter or both. They should be contrasted with the non-canonical kinetic terms that effect Vainshtein screening [13]. The two key benchmark models relevant to this paper are the chameleon model [14][15][16], which relies on non-linearities in the potential to make fluctuations more massive in dense environments, and the symmetron model [17,18] (for related earlier work, see Refs. [19][20][21][22][23]), which instead relies on a non-linear potential and a quadratic coupling to matter to vary both the mass of the field and the strength of the coupling to matter depending on the local energy density.
For screening models of this type, we know that, in any given environment, compact objects will either be screened, meaning the fifth force that they source or experience will be weaker than gravity, or unscreened. This can be understood in terms of the scalar "charges" of these objects. If the object is screened then the charge sourcing the scalar field around it will be smaller than the object's gravitational mass due to the so-called "thinshell effect" [15]. If the object is unscreened then the scalar charge will correspond to the gravitational mass.
The condition that determines whether or not an object is screened is typically phrased in terms of the object's mass (or density) and size. However, in many situations, there is a large hierarchy between the (small) spatial extent of the source object and the (large) Compton wavelength of the scalar field. We would like to be able to treat these sources as effectively point-like, but, to understand their screening properties, it is currently necessary to know the details of their structure. A deeper understanding of screening for compact and point-like sources is particularly relevant to laboratory experiments that use small test particles, such as atoms [24][25][26], neutrons [27,28] or micro-spheres [29], where the scalar Compton wavelength may be of order the size of the vacuum chamber; and to astrophysical environments, where we need to know the screening properties of stars and planets, but the scalar Compton wavelength may be of order galactic scales [30]. Herein, we have differentiated between point-like and compact sources. While the former refers only to the spatial extent of the source, we consider a compact source to be one that is both small, compared to the Compton wavelength of the field, and high density.
A comprehensive treatment of the screening of point-like objects also represents an important first step towards a robust understanding of when distributions of discrete sources may be "averaged over" and treated as a single extended source with some continuous density.
Such an understanding is also important for knowing how to include the effects of screening when computing interparticle forces in N-body simulations, e.g., of structure formation in the universe [31][32][33][34][35][36]. In this work, we consider the screening of fifth forces around compact, point-like objects in flat space, paving the way to further study in curved spacetimes, such as around black holes [37][38][39][40].
Exact analytic solutions can often be obtained for systems that vary only in one spatial dimension, having applications to experiments that take place between parallel plates, such as torsion-pendulum, Casimir, and neutron experiments [42,[45][46][47][48], as well as the nested cylinders of the MICROSCOPE mission [49]. The solution for a point source in one spatial dimension is also known analytically for particular choices of potential and coupling [50].
Beyond these examples, however, our knowledge is more limited. Therefore, in this work, we perform a detailed analysis of the impact of the dimensionality of the system on the problem of taking the point-like limit of an object in models with screening.
In this article, we aim to solve for the static configurations of the scalar field around spherical and cylindrical matter sources in the limit where the radius of the source tends to zero. We will see that it is important to differentiate between the point-like limit where the spatial extent of a source is sent to zero, while holding its mass fixed, and the compact limit where the size of the source is shrunk and the mass of the source is increased. Non-trivial static vacuum solutions of non-linear classical field theories are well known, and they are often referred to as (topological) defects. These solutions may be stable if they're supported by some non-trivial topology of the vacuum manifold in field space [51], and their stability properties can be characterised using Derrick's theorem [52]. The solutions that we consider in this work differ in that they are supported by matter distributions. We present a new extension of Derrick's theorem, which, for the first time, characterises the stability properties of field profiles supported by sources. Moreover, in three spatial dimensions, it allows us to infer the behaviour of the scalar field profiles in the point-like limit. While we might expect the point-like limit to be divergent, carefully studying the limiting procedure allows us to understand how these divergences are regulated 1 and to determine how the fifth forces due to compact, point-like sources are screened.
The key results of this article are: • An extension of Derrick's theorem to scalar profiles supported by external sources.
• Piecewise analytic approximations for finding the field profiles around point-like sources.
• A numerical treatment of the scalar field profiles around point-like and compact sources. 1 A similar approach is taken in theories of infinitesimally thin branes in higher-dimensional spacetimes [53,54].
• A comparison of the screening factors between the full numerically calculated field profiles and the approximate solutions that are currently used in the literature to quantify the amount by which the fifth force is screened.
A summary of the behaviour of screening in this limit for the two prototypical theories that we study, and for spherically and cylindrically symmetric sources, can be found in Table II. This paper is organised as follows. Section II provides an overview of the particular models that we will examine in later sections, as well as an introduction to the screening factors that have previously appeared in the literature. Section III describes how Derrick's theorem is modified in the presence of classical sources and discusses its implications for the behaviour of the scalar field profiles as the source becomes increasingly compact. In Sec. IV, we turn our attention to the computation of the field profiles. Section IV A focuses on improved piecewise approximations to the potentials and field profiles. Section IV B subsequently describes the numerical procedure that we employ to solve for the field profile without making any approximations to the form of its potential. In Sec. V, we describe the implications of our numerical results for screening and estimates of screening around extremely compact sources. In Sec. V C, we comment on the reliability of our classical calculations, establish the regimes in which we trust the classical predictions and identify when quantum corrections may become large. Finally, we provide our conclusions in Sec. VI.
Extra details relating to the numerical procedures and an explicit example of a construction of a piecewise field profile are included in the appendices.
Throughout this work, we use the mostly plus metric signature convention (−, +, +, +). A summary of our notation is provided in Table I.

II. SCREENING MODELS
In this section, we introduce the scalar-tensor theories that will be the focus of this article. The fifth forces that the additional scalar fields mediate are subject to screening mechanisms, which arise due to non-linear terms in their "effective potentials". 2 The models that we consider have an action of the form 2 The term "effective potential" used here should not be confused with the effective potential that includes loop corrections to the tree-level potential.
M mass scale determining the strength of the matter coupling µ mass parameter in the potential of the screening scalar α quartic self-coupling of the screening scalar ρ energy density of the matter source λ s screening factor φ(ρ) position of the density-dependent minimum of the effective potential φ ∞ ≡φ(0) position of the minimum of the effective potential for ρ = 0 r s source radius m s source mass ρ in constant density of the source ρ out constant ambient density external to the sourcē φ in(out) ≡φ(ρ in(out) ) position of the minimum of the effective potential for ρ = ρ in(out) m in(out) effective mass of the scalar fluctuations inside(outside) the source x ≡ µr dimensionless radius x s ≡ µr s dimensionless source radius r shell thin-shell radius d number of spatial dimensions ρ ≡ ρ/ρ c normalised density: ρ c = µ 2 M 2 for the quadratically coupled model and ρ c = µ 3 M/ √ α for the linearly coupled model  which describes a canonically normalised scalar field φ that is subject to a potential V (φ) and conformally (Weyl) coupled to matter, which we denote by the set of fields {ψ}. Herein, R is the Ricci scalar and G is Newton's gravitational constant. Matter fields move on geodesics of the Jordan-frame metricg coupling to the field φ through the coupling function A(φ). Keeping only these terms, the coupling functions of the models that we study are and leading to linear and quadratic couplings to matter, respectively. A linear coupling to matter arises in the chameleon model [14], whereas the quadratic coupling is used in the symmetron model [17]. The omission of a linear term in φ/M in Eq. (4) follows directly from requiring a Z 2 symmetry of the coupling function, i.e., a symmetry under φ → −φ.
Because of the key role that non-linear effects play in their behaviour, we consider theories that exhibit spontaneous symmetry breaking (SSB), akin to the symmetron model [17]. We take the potential for the scalar field to be of the form: where the dimensionless self-coupling α > 0 and squared mass parameter µ 2 > 0 are real.
The constant µ 4 /(4α) ensures that the field has zero energy density in the Minkowski vacuum. While this constant is irrelevant from the point of view of the scalar-field equations of motion, it is relevant for the applicability of Derrick's thoerem, as described in Sec. III.
Spontaneous symmetry breaking can also be realised in the context of symmetron screening [55] via the Coleman-Weinberg mechanism [56], and we expect the behaviour of this variant to be qualitatively similar.
It is possible to find analytic solutions for models with the potential in Eq. (5) in one spatial dimension only. Even a semi-analytic analysis of the behaviour of the solutions around maximally symmetric sources in higher dimensions is extremely challenging. Given the importance of non-linear effects, these models therefore provide an illustrative proving ground within which to develop and apply the techniques discussed in this work.
The coupling of the scalar field to matter means that its dynamics can be understood as being due to an "effective potential" where T µ µ is the trace of the covariantly conserved energy-momentum tensor of matter (see, e.g., Ref. [3]). In non-relativistic environments, matter can be treated as a perfect fluid, and we can write T µ µ = −ρ, where ρ is the energy density of matter. The coupling of the scalar field to matter fields leads to two effects. First, matter fields act, through the trace of their energy-momentum tensor, as a source for the scalar field.
Second, gradients in the scalar field give rise to a "fifth force", such that matter does not move on geodesics of the Einstein-frame metric. The scalar fifth force per unit mass has the form manifesting as a deviation from general relativity. As discussed in the introduction, the non-linearities of the scalar field can make the strength of this fifth force vary depending on the environment, and the effects of the fifth force can be suppressed or screened.
The scalar field theories that we consider in this work have a non-vanishing mass m out exterior to the source. As a result, the fifth force mediated by the scalar fluctuations is subject to a Yukawa suppression at radii larger than 1/m out . However, the focus of the present work is the additional suppression of the fifth force due to non-linear effects in addition to the Yukawa suppression.
In order to quantify screening, it is common to introduce a "screening factor" λ s , and a number of different definitions of this factor have appeared in the literature; see, e.g., Refs. [15,16]. In phenomenological applications, this screening factor is often defined in terms of the ratio of the magnitudes of the scalar fifth force F φ to the Newtonian gravitational force F N . However, we instead wish to use a screening factor in which the impact of nonlinearities is not confounded with the Yukawa suppression described above. To this end, we where F Yukawa is the magnitude of the Yukawa force due to scalar fluctuations of the same (effective) mass but with a quadratic potential. Well inside the Compton wavelength of the scalar field, the screening factor in Eq. (8) is directly proportional to previous definitions in We will now review the standard argument for how this screening takes places around static, compact and spherically symmetric sources, focusing on the quadratically and linearly coupled models described above. Here, "linearly coupled" and "quadratically coupled" refer to the coupling as it appears in the effective potential, not as it appears in the equation of motion.
A. Quadratically coupled model Taking the matter degrees of freedom to be non-relativistic, the effective potential of the quadratically coupled model, with the coupling function given in Eq. (4), takes the form When the local density ρ ≥ µ 2 M 2 , the effective potential only has one minimumφ(ρ) at φ(ρ) = 0. In regions of lower density with ρ < µ 2 M 2 , the effective potential has two minima atφ(ρ) = ±φ ∞ 1 − ρ/µ 2 M 2 , whereφ ∞ = µ/ √ α is the field value at the minimum of the effective potential when ρ = 0. Since the classical fifth force given by Eq. (7) is proportional to the ambient field value, we haveφ(ρ) = 0 in sufficiently dense regions and the fifth force vanishes.
In order to treat the screening behaviour analytically, we consider a compact, maximally symmetric source with finite radius r s , uniform density ρ in and total mass m s , embedded in a diffuse background with uniform density ρ out that is less than the critical value µ 2 M 2 .
In the literature, e.g., Ref. [18], the screening of the fifth force sourced by this object is determined by dividing space into two regions: one inside and one outside the source object.
In addition, we make the following assumptions: (i) the field inside the source remains close to zero and (ii) the field outside the source remains close to the value that minimises the effective potential outside the source, which we denote bȳ where we have arbitrarily (but without loss of generality) selected the positive vacuum solution at spatial infinity. We then approximate the effective potential as quadratic about its minimum in each of the two regions. The effective mass of the field inside and outside the source is given by whereφ in =φ(ρ in ). We impose the boundary conditions that the field be regular at the origin and tend to a constant at infinity, and that φ and dφ/dr be continuous at the surface of the source.
After making the assumptions described above and further approximating the fifth force as F = −φ ∇φ/M 2 ≈ −φ out ∇φ/M 2 , 3 the external field profile around our object, which we centre at the origin, can be written as [57] where the screening factor λ s is given by The cylindrically symmetric (i.e., two-dimensional) case can be treated in a similar manner, with the exterior field falling off as the modified Bessel function of the second kind K 0 (m out r): with where we have again assumed that the fifth force can be approximated as −φ out ∇φ/M 2 .
Lastly, we review the one-dimensional case that was analysed in detail in Ref. [50]. For the quadratically coupled model, the equation of motion takes the form where ϕ = φ/φ ∞ , x = µr, x s = µr s and primes indicate derivatives with respect to x. The exterior solution is and the interior solution can be written in terms of the Jacobi elliptic function nc as where The values of the field ϕ at the origin and at the surface of the source, ϕ(0) and ϕ(x s ) respectively, are obtained by matching the solutions and their first derivatives at x s , which can be done only semi-analytically for a source of finite radius.
In the point-source limit, x s → 0, we have ρ(x) = m s δ(x), and the solution reduces to with The screening factor for a point source is thus with the limiting behaviour and we see that the fifth forces sourced by more massive objects are screened.

B. Linearly coupled model
For the case where the scalar coupling to matter is linear [cf. Eq. (3)], we again take the matter degrees of freedom to be non-relativistic so that the effective potential is Unlike the quadratically coupled model, there is no critical density or phase transition between screened and unscreened regimes. In this case, the mass of fluctuations in the scalar field increases with the ambient matter density. This makes the perturbations, and hence the force, sourced by regions deep inside a large dense object short-ranged [16] .
Before continuing, we remark further on our truncation of the expansion of A(φ) to leading order in φ/M 1. At leading order, we have dropped the term ρφ 2 /(2M 2 ), since it is subleading compared to ρφ/M . However, the power counting in φ/M 1 does not ensure that ρφ 2 /(2M 2 ) is subleading compared to the mass term −µ 2 φ 2 /2. In fact, this will not be the case for much of the parameter space that we study later in this work. We therefore need to justify keeping the mass term in the potential, while dropping the quadratic term from the expansion of the coupling function. One possibility would be to consider a coupling function with two scales: M 1 , which controls the odd power series in φ, and M 2 , which controls the even power series in φ. 4 We can then arrange for the quadratic term in the expansion of A(φ) to be subleading compared to the mass term in V (φ) without affecting the magnitude of the linear coupling. Similar fine-tuning arguments apply in respect of the quartic operator in V (φ), and the cubic and quartic operators arising from the expansion of A(φ).
Considering again a uniform-density, maximally symmetric source of mass m s and radius r s , embedded in a low-density background, the standard approach to computing the screening factor is to divide space into three regions: one outside and two inside the source [15]. In the innermost region, the field is so massive that the field value is approximately constant.
The surface of this region defines a "thin shell" radius r shell with 0 ≤ r shell ≤ r s . In the region r shell ≤ r ≤ r s , the effective potential is dominated by the matter coupling.
To leading order in φ/M , the field value that minimises the effective potential inside the source must satisfy and the mass of the scalar fluctuations aboutφ in is The exterior field profile around a homogeneous spherically-symmetric source of radius r s The resulting screening factor is given by and, taking m out r 1, the thin-shell radius is Only the matter within the shell near the surface of the object sources perturbations in the field. This form of the solution holds for all models with a linear coupling to matter, as different self-interaction potentials V (φ) merely lead to different expressions forφ out .
When r shell = 0 and λ s = 1, which occurs for 2M (φ out −φ in ) ≥ ρ in r 2 s , the field cannot reach its minimum anywhere and the object is entirely unscreened. When r shell → r s and λ s → 0, occurring when ρ in r 2 s 2M (φ out −φ in ), the field quickly drops to its minimum near the surface of the source and the object is maximally screened.
The cylindrically symmetric or two-dimensional case is similar, with [58] and the thin-shell radius r shell is given by where and W k denotes the k th real branch of the Lambert W function. The argument of the W function becomes non-negative when X < −Y , and Eq. (31) no longer gives an r shell between 0 and r s . As the argument approaches zero, the W function tends to −∞, r shell → 0, and the source is unscreened.
In the one-dimensional case, the equation of motion takes the form The interior solution is not of a simple form, nor is it particularly illuminating, and we therefore omit it here. On the other hand, the exterior solution is the same as in Eq. (17), so, in the point-source limit, we have [50] We encounter a branch point in the square root at m s = √ 2µMφ ∞ , at which the solution becomes imaginary, signifying a breakdown phenomenon [50]. As we will see in Sec. IV, this behaviour is not observed in the linearly coupled model in spatial dimensions d > 1 .
The screening factor of a point source is thus with the limiting behaviour We are interested in stationary localised solutions, which must correspond to a minimum of the energy. For the quadratically and linearly coupled models that we study here, the energy is given in d spatial dimensions by where n = 1, 2. This energy is finite as long as we shift the potential so that the field has zero energy density in vacuum, see Eq. (5). We can write the contributions from the kinetic, potential and matter-coupling parts respectively as For a spherically symmetric, uniform source with ρ(x) ≡ ρ(r) = ρ in θ(r s − r), where ρ in is a constant and θ is the unit step function, we can write where dΩ d is the integral over the solid angle subtended by the d-dimensional sphere.
In order to determine whether any given solution φ is stable (with respect to growth or collapse), we are interested in the first and second variations of the energy with re- Writing We now suppose that we want the field profile to be generated by a point source of finite mass, with r s → 0. As long as φ and its spatial derivatives remain finite, this removes the > 0 for theories with a canonical kinetic term and spatially varying profiles, this is satisfied in d = 1 and 2, but not in d > 2. Therefore, for any model with a canonical kinetic term, the derivative of any stable and stationary field has to diverge at r s in d > 2 as r s → 0.
Whether the field itself also diverges cannot be derived from the arguments above, but we might expect it to depend on the form of the potential and the matter coupling, such that the field will tend towards the minimum of its effective potential. This may or may not diverge as ρ → ∞, depending on the choice of V (φ) and the value of n. With the choice of bare potential in Eq. (5), an even value of n brings the minimum of V eff (φ) to zero as ρ → ∞, while an odd value causes it to diverge, as we will see later in Sec. IV B (see Fig. 2).
In the remainder of this work, we take d to be the effective dimensionality of the source.
The case d = 3 corresponds to a spherically symmetric source in three spatial dimensions; the case d = 2 instead corresponds to a cylindrically symmetric source in three spatial dimensions, i.e., a cylinder with infinite extent along its symmetry axis. For the latter case, we assume that the solution φ(r, z) in three spatial dimensions has a separable solution, φ(r, z) = φ(r)Z(z) (with some abuse of notation), such that we can integrate out Z(z) to obtain an effective two-dimensional action, absorbing the mass scales and separation constant into the model parameters, such that φ(r) is dimensionless (having therefore canonical dimensions within a two-dimensional action).
Returning to our generalisation of Derrick's theorem, it is easily verified that it is satisfied for the one-dimensional cases described in Secs. II A and II B. For the quadratically coupled model, we have from Eq. (42a) that where are the dimensionless contributions from the kinetic and potential terms, where the constant ensures that the field has zero energy density in vacuum. Their values are and we therefore require which is indeed the case. We also see explicitly that For the linearly coupled model, since the solution is of the same form as the quadratically coupled case, it is easily shown that Derrick's theorem is again satisfied, as long as At the branch point in Eq. (34), the energy becomes complex, signalling an instability, and Derrick's theorem breaks down. However, this feature is unique to one spatial dimension.
In this work, Eqs. (42a) and (42b) will be used to perform quantitative checks of the numerically obtained field profiles (see Sec. IV B). However, as we have seen, they provide us with a qualitative prediction for the behaviour of those numerical solutions in the pointlike limit in d = 3.

SOURCES
In this section, we will describe how to compute the static scalar field profiles around sources as they become more compact and point-like. We begin by describing how the piecewise analytic approach to solving the field profiles can be extended to more compact sources. We then proceed to describe our numerical approach to computing these profiles, and the properties of the solutions that we find. We will consider solutions around an infinite cylinder (d = 2) and around a sphere (d = 3). Both sources will have a radius r s and uniform density ρ, such that ρ in = ρ inside the source (r ≤ r s ) and ρ out = 0 outside the source (r > r s ). We note that ρ is a mass density with units mass/length 3 for a sphere and units mass/length 2 for an infinite cylinder.
Since we consider extended sources of finite density, we demand that the field profile is regular at r = 0 and that it approaches one of the field values that minimise the effective potential as r → ∞. In vacuum, the effective potential for both models is invariant under φ → −φ, and the field is equally likely to sit in the minimum at positive or negative φ. If the field couples to matter quadratically then the symmetry under φ → −φ is preserved in the matter coupling, and, up to an overall sign, the form of the field profile is insensitive to the choice of vacuum at r → ∞.
The linearly coupled model requires more careful consideration, as the symmetry is broken explicitly inside the source by the matter coupling. For certain combinations of boundary conditions and sources, we will be unable to find piecewise analytic solutions. For example, if the boundary condition is φ → −φ ∞ as r → ∞ then the field may enter a regime close to the surface of the source where the φ 4 term of the full effective potential dominates. It will not be possible to find a piecewise analytic solution in this case.
In what follows, and for the linearly coupled model, we will consider configurations where the field is in theφ ∞ vacuum at infinity but where the field becomes large and negative inside the source. This raises the question of whether such a configuration is truly stable or whether it can tunnel to nucleate a bubble of the −φ ∞ vacuum around the source. We leave further investigation of this interesting possibility for future work.

A. Piecewise approximations
In Sec. II, we outlined the approach that is common in the literature to finding piecewise analytic solutions for the scalar field profiles around compact objects. For the quadratically coupled model, we divided space into two parts: one inside and one outside the source; for the linearly coupled model, we further divided the space inside the source into two regions.
We will refer to these as two-and three-part solutions, respectively. When we outlined these approximations to the field profiles and the corresponding screening factors, however, we failed to check whether the values of the field in each spatial region were consistent with the approximations made to the effective potential. It is easy to see that for more point-like sources, these approximations do fail, and the solutions in Sec. II are inconsistent. For a wide range of sources, however, we can continue to make successive linear and quadratic approximations to the scalar potential, before matching the solutions to these approximate equations together.
The effective potentials in Eqs. (9) and (24) depend on three fixed scales, which we take to be the Compton wavelength 1/µ (in vacuo), the vacuum expectation value (vev) φ ∞ = µ/ √ α and the matter coupling ρ/M n (n = 1, 2, depending on whether the coupling to matter is linear or quadratic in the scalar field). This gives us the freedom to make two rescalings of the equations of motion and express them only in terms of dimensionless quantities. The solutions will depend on the dimensionless distance x = µr (so that the source radius is x s = µr s ) and the dimensionless source densityρ := ρφ n−2 ∞ /µ 2 M n . For the quadratic coupling, this dimensionless density isρ = ρ/µ 2 M 2 ; it is the ratio of the source density to the critical density. For the linear coupling, it isρ = √ αρ/µ 3 M .
When taking the point-source limit x s → 0, we keep the source mass fixed. We therefore parametrise the source mass asm s :

The effective potential
We now briefly outline when linear and quadratic approximations can be made to the effective potential for our models. We begin by considering the effective potential outside the source, which is the same for both the quadratically and linearly coupled models.
Requiring continuity of the effective potential and its first derivative outside the source, we find that the effective potential is well approximated by As a result, we can consider up to three spatial regions outside the source.
Inside the source, ρ > 0, and we must differentiate between the behaviour of the linearly and quadratically coupled models.
2. Linear coupling. For the linearly coupled model, there are two cases to consider:ρ < √ 12/9, when the effective potential inside the source has two non-degenerate minima; andρ > √ 12/9, when the effective potential has only one minimum.
When the effective potential has two minima, the first minimum φ + satisfiesφ ∞ / √ 3 < φ + <φ ∞ (moving from the upper to the lower limit of this range as the density increases). The second minimum φ − lies in the range φ − < −φ ∞ (decreasing monotonically asρ increases). The only situation that we can approximate analytically is when φ remains close to φ + inside the source, in which case we can approximate the effective potential as a quadratic around φ + , i.e., When only one minimum is present, the field may pass through a region where |φ| ρ 1/3φ ∞ before it is able to reach the minimum of the potential. This region defines a "shell" near the surface of the source. In this case, we find where whereφ(ρ) :=φ(ρ)/φ ∞ . Although not obvious from the expression, ζ < 0 by construction, sinceφ(ρ) < ζφ ∞ < 0. Asρ → ∞, ζ → −ρ 1/3 /2 =φ(ρ)/(2φ ∞ ).
The construction of Eq. (50) is similar to the much simpler construction of the linearly coupled solutions in Sec. II B, where instead only theρ term is kept in the shell and the interior region has a constant φ(r) =φ(ρ). By Eq. (29), the equivalent of ζ is then 1 − ρr 2 s /(2Mφ ∞ ). Note that, in the case of quadratic coupling, we only have one approximation to the potential and one region inside the source; for the linearly coupled model, we have at most two. This is the same as in the calculation of the screening factors outlined in Sec. II. However, in what follows, we will be careful to check the validity of these approximations. Figure 1 provides examples of piecewise continuous field profiles around spherically symmetric sources (in three spatial dimensions) and cylindrically symmetric sources (in two spatial dimensions) for both the quadratically (Fig. 1a) and linearly coupled models (Fig. 1b).

Piecewise scalar field profiles
The constants of integration in the piecewise solutions to the scalar equations of motion are determined by ensuring continuity of the field and its first derivative at the matching surfaces. An example of this construction is given in Appendix A. For cylindrically symmetric sources, analytic solutions to the equations determining the constants of integration are not available, and they must instead be solved numerically. The source parameters for the various examples shown in Fig. 1 were chosen to give sources of sufficient compactness to require a two-part (solid blue), three-part (solid orange) or four-part (solid green) piecewise construction. The dotted lines are the "standard" two-(for quadratic coupling) and three-part (for linear coupling) approximations to the field profiles around a source of the same mass and radius, as outlined in Secs. II A and II B, and references therein. For the linearly coupled case, the sources haveρ > √ 12/9 so that the effective potential inside the source only has a single minimum. By inspection, we see that the approximations for the linear coupling fail to capture the full behaviour of the field profiles, and this worsens as the sources become more compact.
One can see from these profiles that the qualitative behaviour is the same around spherically and cylindrically symmetric sources: for the quadratic coupling, the field at the centre of the source tends towards zero as the source density increases; for the linear coupling, the field at the centre of the source becomes more negative as the source becomes more compact. x s and the dimensionless source massm s are chosen so that we require a two-part (solid blue), three-part (solid orange) or four-part (solid green) solution. The quartic self-coupling was taken to be unity, i.e., α = 1. The dotted lines give the corresponding approximate solutions constructed as described in Secs. II A and II B. In the quadratically coupled case, the solid blue and dotted blue lines are indistinguishable. We remark that there is generally good agreement with our multipart solutions in the quadratically coupled case. On the other hand, we see that the standard approximate solutions deviate significantly from our multi-part solutions in the linearly coupled case and especially for cylindrically symmetric sources, i.e., in d = 2.
We begin to see indications of the behaviour of the field profiles around extremely compact sources. For the quadratic coupling, the field value at the origin decreases towards zero as the source becomes more compact. In three spatial dimensions, as the source is made more compact, the field profile becomes flatter in the exterior of the source, and the field rises more steeply in the vicinity of the surface of the source. This matches our conclusions in Sec. III that the gradient of the field should diverge at the origin, while the field itself should go toφ(ρ) = 0. This will be confirmed by our numerical solutions in Sec. V A. Also in agreement with our conclusions in Sec. III, we see in the linearly coupled model that the field evolves through φ = 0 and becomes negative for compact sources, and our calculations indicate the possibility that the value of the field at the origin diverges in the point-source limit. We will explore this further with our numerical field profiles in Sec. V A.

B. Numerical calculation of field profiles
The piecewise approximations to the field profiles, derived in the previous section, give us some intuition for how the quadratically and linearly coupled models behave around increasingly point-like sources. In this section, we compute numerical solutions to the full equations of motion for both the linearly and quadratically coupled models around cylindrical and spherical sources.
To obtain them, we have modified the ϕenics code 5 , a numerical code building on the FEniCS library 6 [59][60][61] and applying the finite element method to problems of screening.
The modified numerical code, as well as Jupyter notebooks for the main results presented in this paper, are available at https://github.com/benjthrussell. The numerical method underlying the ϕenics code has been extensively described in Ref. [62]. In the following, we limit ourselves to summarising its main qualitative aspects, while focusing on the details of the implementation specific to this work. We conclude by showing sample field profiles around highly compact sources.
For computational convenience, we work with dimensionless quantities: further details can be found in Appendix A of Ref. [62]. However, in this section, we will express all equations in their native form for readability.
The equations of motion can trivially be rewritten in terms of ψ = φ −φ ∞ , for which the Dirichlet boundary condition φ(r → ∞) = µ/ √ α becomes ψ(r → ∞) = 0. We have developed both solvers, i.e., one for the field φ and another one for the field ψ. We find the latter to be more accurate, thanks to its better characterisation of the exponential decay at large r.
We solve the non-linear Eqs. (52) and (53) using Newton's method, an iterative technique that aims to find successively improved approximations {φ (k) , k = 1, · · · , n iter } to the true solution φ, given some sufficiently close initial guess φ (0) , where n iter is the number of iterations. If φ (k) is sufficiently close to the true solution φ, we can expand the non-linear terms in the equation of motions to first order, so that they become linear equations in φ.
For the case of quadratic coupling, we can approximate Eq. (52) at some iteration k + 1 by and, for linear coupling, the solution to Eq. (53) can be found from We solve the discretised form of Eqs. (54) and (55), as described in the next subsection, to obtain an improved approximation φ (k+1) , which will be used for the subsequent iteration.
The {φ (k) } will not exactly satisfy the continuum Eqs. (54) or (55), but the residuals will decrease with the iterations until a convergence criterion has been met, or the solution has saturated the accuracy that the scheme is capable of achieving. The convergence criteria that we use are discussed in Ref. [62].
The finite element method works by solving a discrete approximation to the original continuous boundary value problem. We consider a finite interval of the radial domain r/r s ∈ [0,r max ], wherer max is taken to ber max = 10 4 /x s (so thatr max r s = 10 4 /µ, i.e., 10 4 times the field's Compton wavelength) when x s < 1, andr max = 100 when x s ≥ 1. This domain is discretised into a collection of N intervals {[r i ,r i+1 ], i ∈ 0, · · · , N }, withr i=0 = 0 andr i=N +1 =r max , which we will refer to as "cells". We will additionally refer to the collection of all cells as a "mesh" and to cell extrema as mesh "vertices". Details of the meshes that we use are given in Ref. [62].
A key aspect of the finite element method is that cells can be of varying size. In other words, mesh vertices do not need to be equally spaced. This is particularly helpful in addressing both the hierarchy between the Compton wavelength and the source radius, and In this work, we are interested in spherically or cylindrically symmetric functions, and we use the inner product: which is the Euclidean inner product in d dimensions except for a multiplicative geometrical factor. For the test space, we use the space of square-integrable functions on R d , whose gradients are also square-integrable on R d . This integral form is particularly adept at handling sharp transitions like the step source profile for which we intend to solve.
In this form, the Newton iteration for the quadratically coupled model in Eq. (54) becomes where we have used the Neumann boundary condition dφ(r → 0)/dr = 0 and integrated by parts the term ∇ 2 φ (k+1) v r d−1 dr to lower the order of the first differential operator. 7 It is not always desirable to employ globally continuous function spaces; see, e.g., Ref. [62].

Initial guess
In Sec. IV A, we showed how to approximate the scalar field profile around a compact source, which could serve as an initial guess for our iterative Newton's solver. We use the conditions on the matching radii, derived for the piecewise approximations of Sec. IV A to determine how many piecewise approximations we need to match together to find a good initial guess for a given set of model and source parameters. Whilst it would be possible to also use the analytic expressions previously determined as piecewise approximations in these regions, we find it numerically more efficient to use FEniCS to solve the linearised equations of motion in each region in order to find a good initial guess. The resulting initial guesses are an excellent match to our analytically derived piecewise profiles. The coefficients of the linearised equations of motion solved by FEniCS to find the required initial guess are given in Appendix B.
For the quadratically coupled model, whenρ = 1, the approximation described in Sec. IV A fails. In this case, we artificially detune the physical parameters away from the critical point ρ = (µM ) 2 for the purpose of obtaining an initial guess for the algorithm.
Further details are given in Appendix B. We stress that the initial guess only needs to be accurate enough to allow the solver to converge to the correct physical solution; it need not be an accurate reflection of the final profile, although a more accurate initial guess will ensure faster convergence to the true solution.

Field profiles
To ensure the validity of our numerical solutions, we perform all the convergence tests described in Appendix B of Ref. [62]. Moreover, we perform the Derrick's theorem tests of Eqs. (42a) and (42b). We have checked that our results do not depend on the specific source profile used for the density: a step function, a smoothed top-hat profile or a Gaussian profile. We present results for the step function in this article, since these can be compared with the analytic results in Eqs. (41a) and (41b), derived from the generalisation of Derrick's theorem.
For choices of parameters where we trust the validity of our piecewise initial guesses, such as those shown in Fig. 1 and the orange lines in Fig. 2, we find that the field profile found by our solver is always close to that initial guess. For the quadratically coupled model around spherically symmetric sources (or in three spatial dimensions), we find the limiting behaviour that φ(r s ) → 0 for more compact sources, as we might expect. For cylindrically symmetric sources (or in two spatial dimensions), we find the same limiting behaviour, but shallower gradients, with the field reaching zero more slowly, cf. Table II.
For the linearly coupled model around spherically symmetric (three-dimensional) sources, we find that φ(r = 0) ≈φ(ρ) = −φ ∞ρ 1/3 for small r s . For cylindrically symmetric (d = 2) sources, we find that that φ(r = 0) becomes very negative, although the particular scaling with the source parameters is more complicated, as we describe in Sec. V A. This confirms our hypothesis that the field at the origin diverges near a point source in the linearly coupled model. Examples supporting all the above inferences are shown in Fig. 2.

SOURCES
In this section, we describe the key implications of our results for the screening of dense objects whose spatial extent is significantly smaller than the Compton wavelength of the fifth-force mediator. We summarise these implications in two ways. First, we provide a set of scaling relationships, which capture the behaviour of the scalar field and its boundary conditions at the surface of the source. These allow our results to be extrapolated to even more compact sources. Second, we compare the screening factors obtained from our full numerical solutions with those obtained from the piecewise approximations usually applied in the literature. Where these screening factors disagree, our full numerical results provide refined screening conditions, which will find utility in phenomenological studies. Finally, and before offering our conclusions, we draw attention to the potential relevance of quantum corrections. in the linearly coupled case (blue), only the numerical solution is plotted, since the piecewise approximate solution overshoots the minimum to large negative values. In this case, the field value |φ(r s )| is far greater than the vevφ ∞ and is therefore outside the range of validity for the linearised potential in Eq. (47). Nevertheless, the numerical results agree with our expectations, such that the field reaches a minimum value with φ/φ ∞ = −ρ 1/3 at the centre of the source. This relationship holds for other parameter choices, allowing us to confirm that inside dense sources, the field generally reaches the minimum of the effective potential at the origin. Elsewhere, we see that there is excellent agreement between the full numerical profiles and the four-part piecewise constructed profiles.

Coupling
Dimensions  Table I). We reiterate that "two dimensional" versus "three dimensional" refers here

A. Scaling relationships
The piecewise approximations to the field profiles derived in Sec. IV, along with the inferences from applying our generalisation of Derrick's theorem, indicated the possible behaviour of the field profiles as the compactness of the sources increases. By complementing these results with full numerical solutions, we are able to analyse behaviour in regimes beyond the applicability of the piecewise approximations. By these means, we can extract a set of scaling relationships that characterise how the field value φ(r s ) and the field gradient dφ(r s )/dr at the surface of the source vary with source mass and radius. This allows us to see the scaling in both the point-like limit, where x s → 0 holdingm s fixed, and the compact limit, where x s → 0 andm s → ∞. These scaling relationships for compact sources are presented in Table II. Together, they provide a guide to the boundary conditions that should be applied at the surface of the source, allowing our results to be extrapolated to still more compact sources without the need to solve for the full behaviour of the field.
The key observations can be summarised as follows: Linear coupling, d = 3: Both the field and its gradient diverge at the surface of the source as the source becomes more compact. We note that the value of the field at the surface of the source scales with ρ 1/3 , as discussed in Section IV B and Fig. 2.
Linear coupling, d = 2: The field diverges at the surface of the source as the source becomes more compact. The behaviour of the gradient of the field at the surface of the source depends on the way in which the limit is taken. In the point-like limit, wherẽ m s is held fixed and x s → 0, the gradient of the field at the surface of the source tends to a fixed value. If x s is held fixed and the mass of the source is increased, the gradient of the field at the surface of the source diverges.
Quadratic coupling, d = 3: The field tends to zero at the surface of the source as the source becomes more compact. In the compact limit, the gradient of the field at the surface of the source diverges as 1/r s , in agreement with our expectations from Derrick's theorem.
Quadratic coupling, d = 2: The field tends to zero at the surface of the source as the source becomes more compact. The field gradient is independent of the source mass and scales as a function of the source radius only.
The above behaviours agree with the inferences drawn from Derrick's theorem in Sec. III in three spatial dimensions. Derrick's theorem does not allow similar inferences to be drawn in two spatial dimensions. We note that in the limit that the source becomes point-like, the field profiles may not be continuous functions, but instead be distributions, reflecting the nature of the delta-function source. We also note that, for d > 1, we find no indication of the breakdown phenomenon identified in Ref. [50] in d = 1 for the linearly coupled model (see also Sec. II B).

B. Screening factors
Having obtained full numerical solutions, we can determine screening factors without approximation. This allows us to validate existing calculations based on approximate analytic solutions, such as those presented in Sec. II.
In Fig. 3, we compare the analytic expressions for the screening factors that appear in Eqs. (13) and (28) with those obtained from Eqs. (12) and (27) using the field gradients from the full numerical solutions. We refer to the former as λ approx and the latter as λ full .
We note that λ full depends on the distance from the source but tends to a constant at large distances. We therefore compare the screening factors at a radius of 10/x s , which proves to be sufficiently distant for all cases. For positive-definite and perturbative values of the quartic self-coupling α, we expect results qualitatively similar to those shown in Fig. 3.
The comparisons of the screening factors can be summarised as follows: Quadratic coupling (Fig. 3a): We find that λ full ≥ λ approx , such that the approximate analytic solutions generally overestimate the degree of screening. This is because the full numerical field profiles have a smaller gradient in the vicinity of the source compared to the analytic approximation. As a result, the full numerical field profiles must have a larger gradient at greater distances from the source in order to reach the vev. We also find that the ratio λ full /λ approx depends mostly on the source radius x s in the ρ > µ 2 M 2 region (m s > Ω d x d s /d), tending to 1 as x s → 0, while the source massm s is only relevant for small masses. The disappearance of vertical contours for x s 1 (with the exception of the region around ρ = µ 2 M 2 ) results from the complete symmetry restoration at the centre of the source: once the symmetry is fully restored there, a further decrease in the radius of the source has no effect.
Most notable are the large deviations of λ full from λ approx that appear close to ρ = µ 2 M 2 . At this critical density, the effective mass of the scalar goes to zero inside the source. The potential is then dominated by the φ 4 term, and the analytic approximations become invalid. This leads to a significant enhancement in the ratio λ full /λ approx .
Linear coupling (Fig. 3b): While we observe contours along which λ approx = λ full (in them s > 1 and x s < 1 region of the parameter space), there is poor agreement between λ approx and λ full across the remainder of the parameter space. In contrast to the quadratically coupled case and for increasingly compact sources, λ approx > λ full , such that the degree of screening is underestimated by the analytic approximations.
In the spherically symmetric case (d = 3), there is a transition at smallm s , where the equation for r shell (29) has no real solutions and, consequently, the approximate screening factor becomes a constant λ approx = 1. For x s → 0, this occurs wheñ m s < 2 3/2 · 4π/3 ≈ 11.8. In the cylindrically symmetric case (d = 2), we see that the ratio of the screening factors becomes independent of the source radii for sufficiently small source mass (see Tab. II).
Overall, for the linearly coupled model, the analytic approximations provide a poor estimate of the degree of screening for compact sources. On the other hand, the analytic approximations provide a reasonable reflection of the true degree of screening for the quadratically coupled model. The important exception to this is the region of parameter space where the source nears the critical density ρ ∼ µ 2 M 2 and for which the effective mass of the scalar inside the source goes to zero. Interestingly, this is also the region corresponding to more diffuse sources, which may have phenomenological implications for astrophysical objects, such as nebulae. We remark that this more strongly non-linear region of parameter space has previously been identified as interesting in the context of disk galaxies [71][72][73]. However, we leave its detailed study to future work. This is not least because quantum corrections may play an important role here (see also Sec. V C), at the very least renormalising the condition that defines the critical density.

C. Quantum corrections
Throughout this work, we have studied the solutions to the classical equations of motion.
Before offering our conclusions, we assess the potential importance of quantum corrections.
Indeed, the field gradients become large in the vicinity of point-like sources, and we might expect quantum corrections to become important. More precisely, if the classical solution is to be valid, we require the quantum fluctuations to be subdominant once smoothed over scales smaller than the characteristic length scale of the classical solution, where the latter is proportional to the Compton wavelength of the field. The arguments that follow are based on those made in Ref. [74] and references therein.
Our aim is to find a smoothing length L that satisfies two criteria: i.e., the length L is smaller than the characteristic length scale of variations in the classical configuration. Assuming that the the field rolls to the vev in about one Compton wavelength gives ∼ 1.
i.e., the quantum fluctuations smoothed over an interval of length L are smaller than the change in the classical solution over the same interval.
Since we will only be considering fields with non-vanishing mass in vacuum, we can always choose L such that (i) is satisfied.
We note that, if the conditions (i) and (ii) are not met, it does not mean that the classical solution is invalid or uninteresting. Instead, this indicates only that one should look more closely at the predictivity of the classical approximation in this regime.
The change in the classical field configuration over the distance L is approximately where we have takenφ ∞ /( µ −1 ) to approximate the derivative in the vicinity of the source, viz. where the derivative is at its largest.
In order to estimate the quantum mechanical fluctuations, we first smear the quantum field over an interval of length L. This can be done by taking the average of the field over all space, weighted by a Gaussian of width L: Placing x at the origin without loss of generality, we find After performing the spatial integrals, we are left with Hence, in order to satisfy the condition (ii) above, we must be able to find an L such that √ α In other words, we must be able to choose L to be large enough that quantum fluctuations are under control. Combining this with the condition (i), we therefore require that (We note that α has dimensions of mass in two spatial dimensions.) From our analytic and numerical results, we have seen that 1 for small, dense sources, and our arguments from Derrick's theorem in Sec. III suggest that → 0 for point sources.
That is to say, although the classical solution varies rapidly in the vicinity of the point source and goes to zero, the quantum fluctuations are sufficiently large that they swamp the classical behaviour. Notice that this is not the case in d = 1, for which we require α µ, independent of , such that we can always find a coupling small enough that the quantum effects are suppressed. A comprehensive analysis of the quantum corrections to the scalar field profiles presented here is beyond the scope of this work.

VI. CONCLUSIONS
The aim of this work was to study the screening of fifth forces around extremely compact sources that are effectively point-like compared to the Compton wavelength of the scalar fifth-force mediator. To this end, we have derived a generalisation of Derrick's theorem in the presence of external sources, which allowed us to infer the behaviours of the scalar field as the compactness of these sources is increased. These behaviours were confirmed in two and three spatial dimensions by analysing the full numerical solutions for the scalar field in models with a spontaneous symmetry breaking potential and a linear or quadratic coupling to the matter source (at the level of the action). We found that fields with a quadratic coupling are driven to zero in the vicinity of point-like sources, whereas fields with a linear coupling diverge. In addition, the numerically obtained field profiles were compared to novel piecewise analytic solutions that extend the piecewise approximations typically applied in the literature.
By these means, we have been able to establish a set of scaling relationships from which the behaviour of the scalar fields and their boundary conditions at the surface of the source can be extrapolated to other sources. These are shown in Table II. In addition, we have undertaken a critical assessment of the validity of existing screening conditions extracted from so-called screening factors. For the linearly coupled model, this has provided a refined estimate of screening over a large region of parameter space. On the other hand, for the quadratically coupled model, we find reasonable agreement with existing screening estimates in the screened regime, except in the region of parameter space where the source has the critical density for which the field effectively become massless and the field equation becomes unavoidably non-linear. Here, we find that existing approximations significantly overestimate the amount of screening. Since this region of parameter space corresponds to diffuse sources, there may be potential implications for astrophysical systems, such as nebulae, which warrant further study. The comparisons between the screening factors extracted from our numerical results and those based on existing approximations are shown in Fig. 3.
We conclude by commenting on the possible future directions highlighted by this work.
We have laid the foundations for studying the conditions under which distributions of discrete sources can be treated approximately as continuous density distributions for the purpose of analysing screening. A natural first step towards addressing this question is the extension of this work to systems of two (or more) point-like sources, as were studied analytically in one spatial dimension in Ref. [50]. Additionally, the results of this work should be extended to curved spacetimes and time-dependent settings. For instance, in the linearly coupled model, we identified that the explicit breaking of the Z 2 symmetry (φ → −φ) by the matter coupling may give rise to an instability. As such, the increasing compactness of a collapsing matter source has the potential to induce a transition in the scalar field from one vacuum external to the source to another. Whether this transition is realised, and the phenomenological implications that it might have, warrants further investigation. Finally, we have identified the regions of parameter space in which it may be necessary to account for quantum corrections to the scalar field profiles in the vicinity of compact sources, since the latter induces large spatial gradients in the field. The treatment of these effects requires us to find the quantum corrections to spatially varying classical field profiles, and the relevant self-consistent procedures for dealing with such calculations have been developed in the context of vacuum decay [75,76]. and x 2 is the solution to (x s − x 2 )(x s + 2x 2 ) − 3(1 + √ 2x 2 ) + 15x s x s (x 2 − x s ) 2 + 3(x 2 − 1)(1 + √ 2x 2 ) + 15x s = tanh νx s νx s .
above are ignored. Similarly, when a three-part solution is the correct approximation, the column r s ≤ r < r 1 is ignored.
As we do not use the full form of Eq. (50), the linearly coupled sources withρ > √ 12/9 are split into high-density (ρ > 1) and medium-density (1 ≥ρ > √ 12/9) regimes. For the former, we make an approximation that the matter coupling dominates over the non-linear term in the potential, such that ∇ 2 φ = −µ 2 φ+ρ/M provides us with a suitable initial guess; for the latter, we take the matter coupling to only slightly displace the field from the vev, such that ∇ 2 φ = 2µ 2 (φ −φ ∞ ) + ρ/M is the equation that we use.
For the quadratically coupled model, whenρ = 1, the analytic approximations described in Sec. IV A fail: in this case, we setν = 10 −10 . To obtain the matching radii, we further set φ(ρ)/m r = 10 −10 , where m r is a rescaling mass that we use to make the field dimensionless (see Appendix A of Ref. [62]). This choice shifts the initial guess slightly away from the critical density ρ = (µM ) 2 , and we find it leads to fast convergence for all the parameter ranges examined in this paper.