The Symbiosis of Vainshtein and Chameleon Screening

Introduction

While recent cosmological observations seem to provide compelling motivation to modify General Relativity (GR) on large scales [1, 2], stringent Solar System tests demand that any viable theory dynamically hide its new degrees of freedom in local dense environments [3, 4, 5, 6, 7, 8, 9]. This process, known as screening, allows models to naturally transition from cosmological-scale deviations to recovering GR locally. The handful of known screening mechanisms can be broadly divided into two categories: i) derivative screening (e.g. Vainshtein [10] and K-mouflage [11]), and ii) environmental screening (e.g. Chameleon [12] and Symmetron [13]).¹ Each one relies on specific nonlinear operators in the Lagrangian to suppress scalar interactions in distinct ways.

Traditionally, these mechanisms have been studied in isolation. However, treating them as strictly independent is an artificial simplification. Generic modified gravity models can simultaneously activate multiple nonlinear screening operators. In particular, dark energy models able to accommodate the phantom-divide crossing favoured by recent data feature either or both operator classes [15, 16, 17, 18, 19, 20], naturally providing the framework for mixed screening regimes. Moreover, combining environmental and derivative operators has been argued to be theoretically advantageous for a consistent ultraviolet completion [21], for stable early-Universe background evolution [22], and for admitting viable dark-energy models [23].² Spherically symmetric equations containing both operator classes have appeared before [25, 26], but their non-linear spatial interplay, and its subsequent impact on observable small-scale structure, has largely been ignored. An initial step in this direction was taken in [14] by deriving a master equation for screening simultaneously including multiple mechanisms.

In this paper, building on that approach, we address this gap by studying the symbiosis between Vainshtein screening as representatives of the two broad screening categories. We introduce an illustrative symbiotic theory, present analytical and numerical solutions for the nonlinear regime, and discuss three phenomenological results: i) the resolution of the void pathology in pure Galileon models, ii) the reduction of the Vainshtein screening radius, and iii) the de-screening (shrinking) of the Chameleon core. Crucially, we demonstrate that this interplay reshapes the spatial fifth-force profiles to be modelled in $N$-body simulations, and thus fundamentally alters the standard observable signatures targeted by upcoming Stage IV surveys [27, 28].

Mixed screening equation

Vainshtein and Chameleon screening are arguably the most well-studied mechanisms for suppressing fifth forces on non-linear scales. The former operates via nonlinear derivative self-interactions, while the latter dynamically increases the effective mass of the scalar through (non-derivative) self-interactions of the field, both of which naturally activate in local, high-density environments. To investigate their interplay, we introduce a minimal Jordan-frame scalar-tensor theory:³

where $\phi$ is the scalar field, $S_m$ is the matter action, $R$ is the Ricci scalar, $g$ is the determinant of the metric, $M_{\mathrm{Pl}}$ is the reduced Planck mass, and $\Box\equiv\nabla_\mu\nabla^\mu$ is the d’Alembertian operator. The physical mass scales are matched to the cosmological background as $\Lambda_3^3\equiv M_{\mathrm{Pl}}H_0^2$ and $\Lambda_2^4\equiv M_{\mathrm{Pl}}^2H_0^2$, where $H_0$ is the value of the Hubble parameter today. These operators span a vast hierarchy ($H_0$–$\Lambda_3$–$\Lambda_2$–$M_{\rm Pl}$). The classical treatment is nonetheless self-consistent, since the sources lie well below the Galileon strong-coupling scale $\Lambda_3$ [38, 39], redressed upward by the Vainshtein mechanism [40, 41]. While the model’s radiative stability is not guaranteed (the Chameleon breaks the shift symmetry protecting the Galileon), our focus is the classical phenomenology of the interplay, itself a concrete instance of the broader regulating role additional operators can play [42]. The dimensionless parameters $\alpha$, $c_2$, and $c_3$ control the strengths of the conformal coupling, Chameleon potential, and Vainshtein kinetic braiding, respectively. This representative model thus combines simultaneous realisations of both Vainshtein (through the cubic Galileon term) and Chameleon screening (through the potential and non-minimal coupling), and is expected to capture the generic phenomenology of their nonlinear interplay, regardless of the specific Chameleon potential or higher-order Vainshtein operators.

To obtain the local spatial dynamics for the scalar, we evaluate the covariant equations of motion on a cosmological Friedmann-Lemaître-Robertson-Walker (FLRW) background. Applying the standard weak-field [43] and quasistatic [44, 45, 46, 47] approximations, we can express the Ricci scalar in terms of the local matter density and the scalar field’s spatial Laplacian (see derivation in Appendix A). Substituting this back into the scalar equation yields the nonlinear field equation:

where $\nabla_{ij}\equiv\nabla_i\nabla_j$ and $\nabla^2\equiv\nabla_i\nabla_i$ are spatial derivatives and $a$ is the scale factor. Here, the linear gradient has absorbed the conformal coupling, shifting the canonical kinetic coefficient. Unlike the perturbative treatment of [14], which expands the potential about a fixed background, Eq. (2) retains the full scalar, so the field is free to leave its background minimum.

Assuming spherical symmetry, we scale the radial coordinate to the source radius ($x\equiv r/R$) and the scalar field to a fixed reference scale ($y\equiv\phi/\phi_\ast$), and normalise the equation of motion so that the potential’s contribution carries unit coefficient ($-1/y^2$). Here, $\phi_\ast$ is a fiducial normalisation (the field value at the effective potential minimum for the reference density $\rho_\ast=M_{\rm Pl}\Lambda_5^2c_2/(\alpha\phi_\ast^2)$). We model the local density with a top-hat profile, but confirm in Appendix § 1.12 the robustness of results to more realistic profiles. Eq. (2) then reduces to the dimensionless ordinary differential equation:⁴

where primes denote differentiation with respect to $x$. For a specified density profile $\widehat{\rho}(x)\equiv\rho_m(x)/\rho_\ast$ (with central and background values denoted $\widehat{\rho}_{\rm in}$ and $\widehat{\rho}_{\rm bg}$), the screening dynamics is completely governed by three independent dimensionless parameters: i) $\varepsilon\equiv2(\lambda_c/R)^2$, the squared ratio of the reference Chameleon Compton wavelength ($\lambda_\ast^2\equiv\phi_\ast^2/(2\Lambda_5^2c_2)$, evaluated at the reference amplitude $\phi_\ast$) to the source radius; note that the general Compton wavelength is density-dependent ($\lambda\propto\widehat{\rho}^{-3/4}$), shorter inside a dense source and longer inside the dilute background ($\lambda_{\rm in}\ll\lambda_{\rm bg}$); ii) the Vainshtein kinetic strength $\beta\equiv c_3\phi_\ast^4/(a^4\Lambda_3^5\Lambda_5^2c_2R^4)$; and iii) the conformal coupling $\alpha$, which modifies the linear gradient term via $\widetilde{\alpha}\equiv(1+6\alpha^2)$. A fourth non-independent parameter, $\xi\equiv2\alpha\phi_\ast/M_{\rm Pl}$, is fixed by $\alpha$ and the reference field amplitude. Since $\xi_{\rm bg}=2\alpha\phi_{\rm bg}/M_{\rm Pl}\ll1$ for the sub-Planckian fields of interest, we set $\xi=0$ throughout our numerical solutions, leaving $\alpha$ to enter solely through $\widetilde{\alpha}$. Throughout the numerical solutions we adopt a representative conformal coupling $\alpha=1/\sqrt{6}$, matched to the canonical $f(R)$ value [48], for which the kinetic renormalisation is exactly $\widetilde{\alpha}=2$. However, note that different $\widetilde{\alpha}$ choices do not affect the resulting phenomenology, only rescale the respective thresholds.⁵

Void pathology resolution

Pure Galileon theories, where Vainshtein is the only active screening mechanism, suffer from a critical breakdown in underdense environments. While the mechanism merely becomes inefficient in mild cosmic voids [49, 50, 51], deep voids trigger an imaginary fifth-force pathology that persists even beyond the quasi-static approximation [52]. Within N-body simulations, this pathology is often bypassed by implementing ad-hoc constant-field cutoffs to prevent numerical breakdown [53, 54]. The requirement that a theory remain stable within realistic cosmic voids therefore seemingly places severe constraints on viable parameter spaces [19, 20, 55, 56].

Here, we demonstrate that the presence of an effective potential admitting a stable, density-tracing minimum naturally resolves this instability. In other words, the symbiosis with the Chameleon cures the pure Vainshtein pathology. To see this, we isolate the spatial gradient in Eq. (3) by integrating from the void centre ($x=0$, $y'(0)=0$) to a radius $x$:

As a quadratic equation for $y'$, we have that⁶

and the reality of the fifth force then relies strictly on the condition $\Delta\ge0$. Substituting the enclosed charge [Eq. (4)] back, the discriminant takes the form of a perfect square,

This identity is general, so whether a void is pathological is decided entirely by what sources $\mathcal{J}$. In a pure Vainshtein model, the scalar is sourced by the density perturbation, $\mathcal{J}_{\mathrm{V}}(x)=\int_0^x s^2[\hat{\rho}(s)-\hat{\rho}_{\mathrm{bg}}]\,\mathrm{d}s$, where $\hat{\rho}_{\mathrm{bg}}=\hat{\rho}(x\to\infty)=\bar{\rho}_{m}/\rho_*$ is the far-field value the profile approaches, set by the cosmological mean density $\bar{\rho}_{m}$. Inside a void the core is underdense, $\hat{\rho}<\hat{\rho}_{\mathrm{bg}}$ (a central density contrast $\delta\equiv\hat{\rho}_{\mathrm{in}}/\hat{\rho}_{\mathrm{bg}}-1<0$), so $\mathcal{J}_{\mathrm{V}}(x)<0$. For a sufficiently deep void this drives $\Delta<0$ and the gradient $y'$ becomes imaginary, hence resulting in a pathology.

The Chameleon changes what sources $\mathcal{J}$. It drives the field towards its local density-tracking minimum $y_{\mathrm{v}}(x)=\hat{\rho}^{-1/2}$, defined exactly where the integrand in Eq. (4) vanishes. The condition $\Delta\ge0$ is exactly the requirement that a real gradient exist, so the void pathology is precisely the loss of such a real branch. By holding the field on this minimum, the Chameleon keeps $y'$ real throughout the void, so $\Delta\ge0$ and the model smoothly bridges the underdensity (Figure 1). This requires only the existence of a real, density-dependent minimum, not the specific form of the potential, so the resolution is generic: any environmentally-coupled potential admitting such a minimum (e.g. a simple mass term) similarly cures the pathology.

This resolution is not unconditional in $\beta$, however. A sufficiently strong Vainshtein term drives the central field far enough from the interior minimum that the real branch is lost and the pathology re-emerges. To locate this threshold, we evaluate the discriminant near the void centre, where the density is $\hat{\rho}_{\mathrm{in}}$ and the field is approximately constant (regular at the origin, $y'\to0$), with central value $y_0\equiv y(x\to0)$. Hence, the integrand of Eq. (4) is constant and the enclosed charge is $\mathcal{J}(x)\simeq(\hat{\rho}_{\mathrm{in}}-1/y_0^2)x^3/3$. Inserting this into Eq. (6) gives $\Delta(x)\simeq x^2[(\epsilon\tilde{\alpha})^2-\frac{8}{3}\beta(1/y_0^2-\hat{\rho}_{\mathrm{in}})]$. In a void, the central field lies below the interior minimum $(1/y_0^2>1/y_{\mathrm{in}}^2=\hat{\rho}_{\mathrm{in}})$, so raising $\beta$ eventually turns the bracket negative, vanishing at the boundary

The same estimate for a pure (shift-symmetric) Vainshtein Galileon gives $\beta_{\mathrm{V}}^{\mathrm{crit}}=3\epsilon^2/[8(\hat{\rho}_{\mathrm{bg}}-\hat{\rho}_{\mathrm{in}})]$. The Chameleon raises the threshold in two ways: (i) the conformal coupling stiffens the kinetic term, contributing a factor $\tilde{\alpha}^2$, and (ii) it replaces the background density $\hat{\rho}_{\mathrm{bg}}$ with the smaller central potential term $1/y_0^2$, raising the threshold by a total factor $\tilde{\alpha}^2(\hat{\rho}_{\mathrm{bg}}-\hat{\rho}_{\mathrm{in}})/(1/y_0^2-\hat{\rho}_{\mathrm{in}})$. Numerical solutions in Figure 1 confirm this boundary, widening the pathology-free window in $\beta$ by about an order of magnitude.

IV Reducing the Vainshtein radius

Having established that the Chameleon stabilises the theory in underdensities, we now turn to dense sources to see how the two mechanisms interact together. This interplay is displayed in Figure 2, which maps the numerical screening factor $\mathcal{S}\equiv|y'/y'_{\mathrm{lin}}|$, i.e. the scalar-field gradient normalised to its unscreened linear value (generalising the analytic Vainshtein screening factor of [54]).

The Vainshtein radius marks where the nonlinear Galileon term overtakes the linear kinetic one in Eq. (6),

where $x_{\mathrm{V+Ch}}$ is the Vainshtein radius of the symbiotic model, different from $x_{\mathrm{V}}^3\approx8\beta\hat{\rho}_{\mathrm{in}}/(3\epsilon^2)$ of the pure Galileon. The net suppression

separates cleanly into a fixed ceiling $\tilde{\alpha}^{-2/3}\le1$ and a dynamical charge reduction $(\mathcal{J}_{\mathrm{V+Ch}}/\mathcal{J}_{\mathrm{V}})^{1/3}\le1$. The Chameleon freezes the screened interior at the minimum $y_{\mathrm{in}}=\hat{\rho}_{\mathrm{in}}^{-1/2}$, which sources no gradient, so the field only has to climb the excursion $\Delta y \equiv y_{\rm bg}-y_{\rm in}$ to reach the background. Two charges then compete to set the Vainshtein radius: the excursion charge $\mathcal{I}_{\rm exc}=\varepsilon\tilde{\alpha}\Delta y$ this climb demands (Appendix C), and the enclosed mass $\mathcal{I}_{\rm V}=\hat{\rho}_{\rm in}/3$ physically available. Since the field sources only the charge needed to reach the background, but cannot source more than its contained mass, the effective charge is the smaller of the two, $\mathcal{I}_{\rm V+Ch}=\min(\mathcal{I}_{\rm exc},\mathcal{I}_{\rm V})$. When the background is dilute the excursion is the binding constraint ($\mathcal{I}_{\rm exc}<\mathcal{I}_{\rm V}$), and the charge, hence the sphere, drops below its pure-Galileon size. The outcome falls into three regimes (Appendix C), shown in the left panel of Figure 3: (i) a reduced but still macroscopic sphere; (ii) a saturation ceiling set by $\tilde{\alpha}$; and (iii) complete elimination. In (i), the still macroscopic sphere ($x_{\rm V+Ch}>1$) shrinks below its pure-Galileon size to

In (ii), once the excursion charge reaches the full mass ($\mathcal{I}_{\rm exc}\ge\mathcal{I}_{\rm V}$), the charge reduction switches off and the ratio Eq. (9) is saturated at the ceiling $\tilde{\alpha}^{-2/3}$. In (iii), when the coupling is too weak to push the sphere past the source it is drawn inside it ($x_{\rm V+Ch}<1$) and Vainshtein screening disappears entirely. Setting $x_{\rm V+Ch}=1$ gives a two-branch boundary [Eq. (23)] (confirmed numerically in Figure 3): an excursion-limited $8\beta\Delta y\lesssim\varepsilon\tilde{\alpha}$ with slope $\propto\varepsilon$, crossing over to a steeper mass-capped branch $\propto\varepsilon^{-2}$ at higher $\varepsilon$. Since $\varepsilon\propto R_*^{-3/2}$, at fixed $\rho_*$ the $\varepsilon$ axis tracks source size; changing $\rho_*$ rescales the $(\varepsilon,\beta)$ axes, not the physics.

V De-screening the Chameleon core

Having seen how the Chameleon reshapes the Vainshtein sector, we now examine the reciprocal effect. In a pure Chameleon model, a screened body has a sharp screening transition: the field is pinned at the high-density minimum $y_{\rm in}$ throughout the interior (where the Compton wavelength $\lambda_{\rm in}\ll R$), and climbs to the background $y_{\rm bg}$ within a few source radii outside it, following a $1/x$ Coulomb tail. We show that the inclusion of the Vainshtein term de-screens (shrinks) that core, which we quantify through the de-screening function

where recall that $y_0$ (the actual central field value) and $y_{\mathrm{in}}$ (the interior potential minimum it would occupy if perfectly screened) are distinct quantities in general. Hence, we have $\mathcal{D}\simeq0$ for a fully screened core sitting at its interior minimum, and $\mathcal{D}>0$ once the central field is lifted from the potential minimum, i.e. the Chameleon fails to activate.

Within a macroscopic Vainshtein sphere, the nonlinear kinetic term caps the spatial gradient [Eq. (5)], flattening the field across $x\lesssim x_{\mathrm{V+Ch}}$, so the interior can no longer reach $y_{\mathrm{in}}$ and $\mathcal{D}$ grows. The right panel of Figure 3 maps $\mathcal{D}$ across the $(\varepsilon,\beta)$ plane, showing that at fixed $\varepsilon$, increasing $\beta$ drives the core from screened ($\mathcal{D}\simeq0$) to eliminated ($\mathcal{D}\to1$).

To derive an expression for the onset of this elimination, we assume a screened core solution exists and match it to the solutions in the transition region and the exterior Coulomb tail. This gives a screened-core radius $x_{\mathrm{s}}^2 = 1 - 2\varepsilon\tilde{\alpha}\Delta y/\hat{\rho}_{\mathrm{in}}$, which shrinks to zero at (Appendix D)

with the numerical value associated with the density contrast in Figure 3. The Vainshtein term relaxes this bound: a larger $\beta$ de-screens the core at smaller $\varepsilon$, so the elimination boundary (blue dotted, $\mathcal{D}=0.1$) bends leftwards in Figure 3. Finally, note that where the core remains screened, the same gradient cap also widens the screening transition, stretching its radial width towards the Vainshtein scale $x_{\mathrm{V+Ch}}$.

Conclusions

Screening mechanisms, essential for the viability of dark energy models, have historically been treated as isolated phenomena. While this has enabled their respective understanding and modelling, such isolation represents an unnatural simplification in general scalar-tensor theories capable of realising multiple mechanisms simultaneously, as is true for many recently proposed dynamical dark energy models. Here, we have derived a unified spatial master equation to analyse the nonlinear symbiosis between Vainshtein and Chameleon screening. Our analytical and numerical solutions demonstrate that this interplay fundamentally alters the kinematics of the scalar field. Specifically, we have shown that the inclusion of a Chameleon potential resolves the imaginary fifth-force pathology inherent to pure Galileon models. By dynamically tracking the effective minimum, the scalar field maintains strict theoretical stability deep within cosmic voids. Furthermore, we have shown a concrete interplay in dense environments. On one hand, the Chameleon mechanism shields the inner core of massive sources, pulling the active Vainshtein radius inwards, eventually eliminating it completely. On the other hand, the Vainshtein operator acts as a kinematic upper bound on the allowable spatial gradient, de-screening the Chameleon core and, where the core survives, broadening its screening transition.

These reciprocal effects carry profound implications for the next generation of cosmological probes. Precision tests of gravity using Stage-IV surveys rely heavily on accurate models and $N$-body simulations. Our results indicate that assuming purely isolated screening mechanisms will systematically mischaracterise the local fifth-force profile in theories that can realise multiple screenings. These results can similarly impact laboratory bounds. Since the Vainshtein term lifts the Chameleon core, exclusion regions mapped for isolated Chameleons [57] need not transfer to theories realising both mechanisms. Ultimately, by rescuing Galileon stability in voids while uniquely altering the structural signatures around dark matter halos, this symbiotic relation opens up previously discarded regions of parameter space, establishing a theoretically robust and observationally distinct phenomenological target for future surveys.

Acknowledgments

We thank Claudia de Rham and Philippe Brax for useful discussions. SS, TB and KN are supported by ERC Starting Grant SHADE (grant no. StG 949572). TB is further supported by a Royal Society University Research Fellowship (grant no. URF\R\231006). CB is supported in part by STFC Consolidated Grants (grant no. ST/T000732/1 and ST/C004235/1). JH is supported by a PhD studentship from UKRI-STFC. We have used escut [58] and xAlpha [59], where we include all the files to reproduce the results of this paper. The authors are grateful for computational resources on the Sciama High Performance Computing (HPC) cluster, which is supported by the Institute of Cosmology and Gravitation, the South-East Physics Network (SEPNet) and the University of Portsmouth. For the purpose of open access, the authors have applied a Creative Commons Attribution (CC BY) licence to any Author Accepted Manuscript version arising from this work.

A Covariant Equations and the Spatial Master Equation

Varying the action in Eq. (1) with respect to the metric $g_{\mu\nu}$ yields the covariant Einstein equation. Taking the trace of this equation and assuming a non-relativistic matter source ($T^{\mu}_{\ \mu} \approx-\rho_m$) gives

In the sub-horizon quasi-static limit, the scalar’s own stress-energy trace is negligible. Being sub-Planckian ($\phi\ll M_{\mathrm{Pl}}$), its $\mathcal{O}((\nabla\phi)^2)$ contributions are dimensionally suppressed relative to the $\sim M_{\mathrm{Pl}}\Box\phi$ term. Approximating the d’Alembertian by its spatial gradients ($\Box\approx\nabla^2/a^2$), the Ricci scalar evaluates to

Varying the action with respect to $\phi$ yields the covariant scalar equation of motion

B Perturbative stability

We analyse linear fluctuations about the static, spherically symmetric solutions of Eq. (3). In the weak-field limit it suffices to perturb the scalar, $y\rightarrow y+\pi$, in the effective Lagrangian whose equation of motion is Eq. (3),

At second order, $\mathcal{L}^{(2)}=-\frac{1}{2}\varepsilon\tilde{\varepsilon}(\partial\pi)^2-\frac{\beta}{2}(\partial\pi)^2\Box y-\beta(\partial y\cdot\partial\pi)\Box\pi$. Integrating the last term by parts, the fluctuation propagates on the effective inverse metric

the cubic Galileon supplying the background-dependent term. For the static background ($\Box y=(x^2y')'/x^2$) its components give the no-ghost condition and the radial and angular sound speeds,

The solutions are stable where all three conditions hold: $\mathcal{G}^{tt}>0$, $c_r^2\geq0$, and $c_\theta^2\geq0.$ Radial stability is governed by the very discriminant $\Delta$ of Eq. (6), so resolving the void pathology ($\Delta\geq0$) simultaneously restores radial stability. Evaluating Eqs. (B3)–(B5) on the numerical solutions across the parameter space (Figure 4), all three conditions are satisfied wherever the Chameleon is active. Only in the deep-Galileon corner (small $\varepsilon$, large $\beta$), where the model reduces to the pure Galileon, do the ghost and the angular instability re-emerge, recovering the known pure-Galileon pathologies [60, 61]. Both $c_r^2$ and $c_\theta^2$ reduce to unity in the linear regime, as required.

C Reduction and elimination of the Vainshtein radius

We derive here the results quoted in Sec. IV. For a deeply-screened top-hat halo, the Chameleon pins the field at $y_{\mathrm{in}}$ throughout the interior, so the entire excursion to the background occurs outside the source. Integrating the exact gradient Eq. (5) from the surface outward fixes the excursion $\Delta y\equiv y_{\rm bg}-y_{\rm in}$ to $^{8}$

where we have used

and $u\equiv x/x_{\rm V+Ch}$ is the radius in units of the Vainshtein radius. The integrand is the Vainshtein enhancement of the gradient over the linear one, i.e. significant only inside the sphere ($u<1$) and vanishing outside. The lower limit $a=1/x_{\rm V+Ch}=R/r_{\rm V+Ch}$ measures how deeply the source sits inside its own Vainshtein sphere, and its two limits correspond to different phenomenology.

Reduction ($a\to0$, $x_{\rm V+Ch}\gg1$). The source sits deep inside a large Vainshtein sphere ($a\to0$), so the full nonlinear region out to $x_{\rm V+Ch}$ contributes to the excursion and $G\to G(0)=\Gamma(1/6)\Gamma(1/3)/(4\sqrt{\pi})\equiv G_0\simeq2.10$, with $\Gamma$ the gamma function. Solving Eq. (21) then gives Eq. (10), showing that the radius grows as $\beta^{1/2}$ yet remains inside the pure-Galileon radius ($\propto\beta^{1/3}$) throughout, saturating at the geometric ceiling $\tilde{\alpha}^{-2/3}$.

Elimination ($a\to1$, $x_{\rm V+Ch}\to1$). The Vainshtein sphere has contracted to the source surface, so the exterior is essentially linear and $G\to1/(2a)$, its linear-tail value (already accurate at $a=1$, where $G(1)\simeq1/2$). Equation (C1) then reduces to

Integrating from the surface gives $y_{\rm bg}-y(1)$; the field departs from $y_{\rm in}$ by only $\mathcal{O}(\varepsilon)$ across the surface layer, so $y(1)\simeq y_{\rm in}$.

$x_{\rm V+Ch}^{3}\simeq8\beta\Delta y/(\varepsilon\tilde{\alpha})$, so the mechanism is eliminated once $8\beta\Delta y\lesssim\varepsilon\tilde{\alpha}$ holds.

Finally, the enclosed charge cannot exceed the full mass, so $\mathcal{I}_{\rm V+Ch}=\min(\mathcal{I}_{\rm exc},\mathcal{I}_{\rm V})$, with $\mathcal{I}_{\rm exc}=\varepsilon\tilde{\alpha}\Delta y$ the excursion charge and $\mathcal{I}_{\rm V}=\hat{\rho}_{\rm in}/3$ the mass. Setting $x_{\rm V+Ch}=1$ in Eq. (8) therefore yields a two-branch elimination boundary,

whose excursion-limited slope $\propto\varepsilon$ and mass-capped slope $\propto\varepsilon^2$ cross where $\varepsilon\tilde{\alpha}\Delta y=\hat{\rho}_{\rm in}/3$. This is the grey boundary in the left panel of Figure 3.

The Chameleon core de-screening threshold

We derive the pure-Chameleon elimination threshold quoted in Eq. (12). Setting $\beta=0$ and $\xi=0$, the master equation reduces to $\varepsilon\tilde{\alpha}(x^2y')'/x^2=\hat{\rho}(x)-1/x^2$, with a top-hat source ($\hat{\rho}=\hat{\rho}_{\rm in}$ for $x<1$, $\hat{\rho}_{\rm bg}$ for $x>1$) and minima $y_{\rm in}=\hat{\rho}_{\rm in}^{-1/2}$, $y_{\rm bg}=\hat{\rho}_{\rm bg}^{-1/2}$. A screened body admits a thin-shell solution in three regions: a pinned core ($y=y_{\rm in}$, $y'=0$) for $x<x_s$; a rolling shell for $x_s<x<1$, where the potential term is negligible against the matter term so that $\varepsilon\tilde{\alpha}(x^2y')'/x^2\simeq\hat{\rho}_{\rm in}$ with general solution $y=\hat{\rho}_{\rm in}^{1/2}(x^2/6+a/x+b)$; and an exterior Coulomb tail $y=y_{\rm bg}+c/x$ for $x>1$ (valid since $\lambda_{\rm bg}\gg1$). Regularity at the core edge, $y'(x_s)=0$ and $y(x_s)=y_{\rm in}$, fixes $a=\hat{\rho}_{\rm in}x_s^3/(3\varepsilon\tilde{\alpha})$ and $b=y_{\rm in}-\hat{\rho}_{\rm in}x_s^2/(2\varepsilon\tilde{\alpha})$; continuity of $y'$ at the surface fixes $c=\hat{\rho}_{\rm in}(x_s^3-1)/(3\varepsilon\tilde{\alpha})$. Continuity of $y$ at $x=1$ then gives

i.e. $x_s^2 = 1 - 2\varepsilon\tilde{\alpha}\Delta y/\hat{\rho}_{\rm in}$. The screened core shrinks as $\varepsilon$ grows and disappears ($x_s \to0$) at

using $\Delta y \simeq y_{\rm bg} = \hat{\rho}_{\rm bg}^{-1/2}$. For $\hat{\rho}_{\rm in}=300$, $\hat{\rho}_{\rm bg}=3\times10^{-4}$, $\tilde{\alpha}=2$ this gives $\varepsilon_{\rm crit}\approx1.3$. Restoring a finite $\beta$ flattens the interior gradient and lifts the core, lowering $\varepsilon_{\rm crit}$, as traced by the boundary in Figure 3.

Numerical method

We solve the dimensionless master equation [Eq. (3)] using the escut [58] solver. We compute this as a first-order boundary-value problem in the pair $(y,\mathcal{J})$, where $\mathcal{J}(x)$ is the enclosed charge [Eq. (4)]. Differentiating yields

the second relation being [Eq. (5)]. This form avoids the singular denominator of the second-order equation for $y$ and keeps the smooth, physical charge $\mathcal{J}$ as a dependent variable. Regularity at the centre imposes $\mathcal{J}(x_{\min})=0$ (equivalently $y'(x_{\min})=0$); at the outer edge we impose a Robin (Yukawa-decay) condition, $y'(x_{\max})+(m_{\rm bg}+1/x_{\max})(y-y_{\rm bg})=0$ with $m_{\rm bg}^2=2/y_{\rm bg}^3$ the background Chameleon mass, so the field relaxes to the background.

We integrate the system with the collocation solver solve_bvp [62], seeded by a minimum-tracking initial guess $y(x)=y_v(\hat{\rho}(x))$ with the consistent $\mathcal{J}$, and continued in the Vainshtein strength: $\beta$ is ramped from 0 (the readily solved pure-Chameleon limit) to its target over several homotopy steps, each seeded by the previous solution. Deeply screened cores and deep voids, where the field spans several decades, are solved in the logarithmic variable $u=\ln y$.

Robustness to the density profile

The main-text results use a smoothed top-hat source. Here, we confirm that the three effects persist for realistic density profiles, matched to the same central and background densities.

Void resolution (Figure 5). We compare the top-hat with a compensated void [63], whose overdense ridge at the edge renders the enclosed charge non-monotonic. For both profiles the pure-Vainshtein discriminant turns negative in the core (an imaginary fifth force), whereas the Vainshtein+Chameleon solution keeps $\Delta\geq0$ throughout. The resolution is therefore generic: it requires only a real density-tracking minimum, not the top-hat shape.

For haloes we test two realistic shapes alongside the top-hat: an NFW-like profile [64], $\rho\propto[(r/r_s)(1+r/r_s)^2]^{-1}$, the cuspy form predicted by cold-dark-matter N-body simulations (regularised with a small inner core so that $\rho(0)$ is finite); and a cored Burkert profile [65], $\rho\propto[(1+r/r_s)(1+(r/r_s)^2)]^{-1}$, which better fits dwarf-galaxy rotation curves. For more details on the cusp/core contrast we refer to [66].

Vainshtein-radius reduction (Figure 6, top-right panel). For all three shapes the symbiotic Vainshtein radius lies below its pure-Galileon value throughout. It equals the geometric value $\tilde{\alpha}^{-2/3}$ at low density and deepens as the growing core charge is increasingly screened (most for the flat-core profiles, least for the NFW-like cusp).

Chameleon-core de-screening (Figure 6, bottom-right panel). The de-screening $\mathcal{D}$ [Eq. (11)] grows with the Vainshtein strength $\beta$ for all profiles, confirming the effect is generic, though its onset depends on the core. The flat-core profiles (top-hat, cored) retain the thin-shell threshold (the core stays pinned ($\mathcal{D}=0$) until $\beta$ exceeds a critical value) whereas the NFW-like cusp, lacking a flat core to pin to, de-screens continuously. The idealised top-hat is thus the conservative case: realistic cuspy haloes de-screen more readily.