If ultralight boson fields exist then vacuum misalignment populates them with nonzero relic abundance. For a broad range of particle mass m the field condenses into fuzzy cores in massive galaxies. We use numerical simulations to test this idea extending previous work Blum and Teodori 2021 and focusing on ultralight dark matter ULDM that makesup a subdominant fraction of the total dark matter density consistent with observational constraints. Our simulations mimic galactic halos and explore different initial conditions and levels of sophistication in the modeling of the halo potential. For msim1025 eV ULDM cores act as approximate internal mass sheets in strong gravitational lensing and could first be detected as an H0 bias in cosmography a scenario we dub AxionH0graphy. The mass sheet degeneracy is broken by finite core radius and by the dynamical displacement of cores from the halo center of mass which introduce imaging distortions and restrict the H0 bias limit of AxionH0graphy to mlesssim5times1025 eV. Cosmological simulations are called for to sharpen the predicted connection between the amplitude of ULDM galactic cores and the ULDM cosmological fraction.
['Astrophysics - Cosmology and Nongalactic Astrophysics', 'Astrophysics - Astrophysics of Galaxies', 'High Energy Physics - Phenomenology']
AxionH0graphy hunting for ultralight dark matter with cosmographic H0 bias
{'AxionH0graphy: hunting for ultralight dark matter with cosmographic H 0 bias': 'Kfir Blum 1, ∗ and Luca Teodori 1, †  \n1 Weizmann Institute of Science, Rehovot 7610001, Israel  \nIf ultralight boson fields exist, then vacuum misalignment populates them with nonzero relic abundance. For a broad range of particle mass m the field condenses into fuzzy cores in massive galaxies. We use numerical simulations to test this idea, extending previous work (Blum & Teodori 2021 [1]) and focusing on ultralight dark matter (ULDM) that makes-up a subdominant fraction of the total dark matter density, consistent with observational constraints. Our simulations mimic galactic halos and explore different initial conditions and levels of sophistication in the modeling of the halo potential. For m ∼ 10 -25 eV ULDM cores act as approximate internal mass sheets in strong gravitational lensing, and could first be detected as an H 0 bias in cosmography: a scenario we dub AxionH0graphy. The mass sheet degeneracy is broken by finite core radius and by the dynamical displacement of cores from the halo center of mass, that introduce imaging distortions and restrict the H 0 bias limit of AxionH0graphy to m ≲ 5 × 10 -25 eV. Cosmological simulations are called for to sharpen the predicted connection between the amplitude of ULDM galactic cores and the ULDM cosmological fraction.', 'II. ULDM cores and cosmographic H 0 bias': 'A. An aside: multiple-source systems', 'IV. Results': '- A. Static external CDM and stars\n- B. Dynamical CDM via m heavy ≫ m\n- C. Results summary for α χ = 0 . 15\n- D. Brief view of α χ < 0 . 15', 'VI. Summary': 'Acknowledgments  \nReferences', 'B. Shear and flexion terms': '13', 'I. INTRODUCTION': "UltraLight Dark Matter (ULDM) could arise via vacuum misalignment [2-4], and the scenario of particle mass as low as m < 10 -20 eV has been discussed extensively in the literature [5-8]. A plethora of axion-like particles spanning a range of masses may occur in string theory [911], motivating the possibility of several ULDM fields, each contributing a fraction of the DM budget.  \n1 2 3 3 4 5 5 6 8 8 10 10 10 13 ULDM differs from collisionless Cold DM (CDM) in several features, allowing to constrain the ULDM massfraction using gravitational probes alone, without relying on model-dependent assumptions. From a cosmological point of view, ULDM suppresses power on small scales, making it possible to constrain it via Cosmic Microwave Background (CMB) anisotropies and galaxy clustering [12], as well as the Lyman-alpha forest [13-16]. In galaxies, ULDM interference patterns characterized by the de Broglie wavelength λ dB ≈ 2 π/ ( mv ) , where v is the characteristic velocity in the system, are seen in simulations [17-28], and lead to detectable effects related to dynamical heating and dynamical friction [6, 29-33]. Another prominent feature of ULDM in galaxies is the formation of cored profiles in the inner part of the dark matter halo [6, 18, 19]. These cores, called 'solitons' in the ULDM literature, are the ground state solution of the Schrodinger-Poisson equations, obeyed by ULDM in the non-relativistic regime. They are consistently seen in simulations, and some of their key features are understood theoretically [24, 34-40]. For m ≲ 10 -21 eV, the ULDM core prediction is in tension with observations of the rotation curves of low-surface-brightness galaxies, if a single ULDM species makes up the full DM budget [41-43].  \nGravitational lensing offers another probe of ULDM. Magnification anomalies can be used to constrain the amount of substructure present in a galaxy [44-47], including substructure due to ULDM [48, 49]. Here we focus on a different effect, suggested in Ref. [1], that concerns the impact of ULDM cores on quasar time delays in strong gravitational lensing systems. Such cores affect the measurement of the Hubble parameter H 0 inferred from time delay cosmography [50], acting as an approximate internal Mass-Sheet Degeneracy (MSD) [51-53].  \nSome studies accounted for the MSD by direct implementation in the analysis pipeline, either via an effective ideal mass sheet [54], or via a particular core parametrization [55]. Such analyses are timely both because they could shed light on the H 0 tension [53, 56-59], and because current and upcoming surveys [60-63] are expected to discover many new time-delay systems [64-67], possibly  \nallowing a < 2% measurement of H 0 with improved stellar kinematics [68].  \nInterestingly, Ref. [54] found some positive preference (at the 2 σ level) for an internal mass sheet component, when combining time-delays and lens morphology prior based on a larger galaxy sample. It is amusing to entertain the possibility that this could be the first observational evidence of ULDM and thus, perhaps, a hint of ultra-high energy physics. More generally, the results of [1, 54] suggest that cosmography is sensitive to features in massive galaxy morphology on 10-100 kpc scales, that are difficult to detect via kinematics, weak lensing, or large-scale cosmological data. As noted in [53], this sensitivity can be significantly improved if one includes external prior constraints on the value of H 0 , coming from other data sets.  \nIn this paper we continue the exploration of ULDM cores in relation to gravitational lensing. Ref. [1] showed that for ULDM cores to be diffuse enough to behave as approximate mass sheets, and massive enough to yield detectable H 0 bias at or exceeding a few percent, requires 10 -25 ≲ m ≲ 10 -24 eV. In this range ULDM is constrained to comprise only a fraction α χ ≲ 0 . 2 of the total DM budget [12, 15]. For such small α χ , simple estimates of the dynamical relaxation time scale [37, 38] led Ref. [1] to identify the time required to form a core as a bottleneck for the lensing effect.  \nWe extend the analysis by using numerical simulations to study ULDM dynamics on the relevant galactic scales. We focus on 10 -25 ≲ m ≲ 10 -24 eV, α χ ≈ 0 . 15 consistent with all other observational constraints, and TDCOSMOlike galaxies [50]: massive ellipticals with M 200 ≈ 10 13 M ⊙ and lens and source redshift configurations typical in strong lensing.  \nOur simulations start from variations around a simple class of initial conditions, where the ULDM roughly follows the same distribution as the CDM component. For these initial conditions, the ULDM in the central region of the halo evolves due to wave pressure to form a flattened core of size ∼ λ dB . The system first expels excess mass from the core region (thus mostly erasing the dependence on the details of the initial conditions in the inner halo). Subsequently, the core density oscillates around a mean level that grows in qualitative agreement with previous related analyses [24].  \nFor most simulated cases, the ULDM spatial distribution for m ≲ 5 × 10 -25 eV is consistent with approximate internal MSD for galaxy lensing applications. We believe that this aspect of the results is robust to variations in initial conditions, and to the detailed modeling of the dominant CDM and stellar mass components of the halo. The equivalent H 0 bias in our simulations ranges between about 1 and 20 percent. For example, for m = 2 × 10 -25 eV (still with α χ = 0 . 15 ) the mean fractional bias across 6 simulations is δ H = δH 0 /H 0 ≈ +0 . 1 with standard deviation of roughly 0 . 06 . In most cases, our estimates suggest that MSD-breaking imaging distortions would be difficult to detect, even without considering degeneracies with other systematic effects like cosmological external shear, sub-halos or host group association. The magnitude  \nof the H 0 bias could, however, depend on the details of the initial conditions at the outskirts of the halo and may be different for cosmological initial conditions, which we are not yet able to perform in this work.  \nThe paper is structured as follows. We describe the approximate MSD in gravitational lensing and the relation to ULDM in Sec. II. In Sec. III we explain the simulations. The main results are presented in Sec. IV, with additional discussion regarding the formation and growth of cores given in Sec. V. We summarize in Sec. VI.  \nApp. A provides more details of the numerical code. App. B contains a discussion of shear and flexion terms.", 'II. ULDM CORES AND COSMOGRAPHIC H 0 BIAS': 'Consider a gravitational lensing system characterized by the lens convergence κ 0 ( ⃗ θ ) , with Einstein angle θ E . The equation relating the angular position of the source ⃗ β 0 with the observed image angle ⃗ θ reads  \n⃗ β 0 = ⃗ θ -⃗ α 0 ( ⃗ θ ) , (1)  \nwhere ⃗ α 0 is the displacement angle, related to κ 0 via 2 κ 0 ( ⃗ θ ) = ∇ θ · ⃗ α 0 ( ⃗ θ ) (and satisfying α 0 ( θ E ) = θ E for a spherically-symmetric system). Changing the model parameters via  \nκ ( ⃗ θ ) = (1 -κ c ( θ E )) κ 0 ( ⃗ θ ) + κ c ( ⃗ θ ) , ⃗ β = (1 -κ c ( θ E )) ⃗ β 0 , (2)  \nkeeps imaging data approximately invariant, and acts as an approximate MSD, provided that κ c ( θ ) is a sufficiently slowvarying function of θ for θ smaller than a few times θ E . (In particular, Eq. (2) expresses an exact MSD if κ c = Const.)  \nFollowing Ref. [53], we consider the case that κ c arises from a cored mass component associated with the main lens, and call κ c a core-MSD component. Following Ref. [1], the physical mechanism we consider to produce the core is ULDM.  \nH 0 bias: Lensing time delays are not invariant under the MSD. If the model of Eq. (2) is the truth model corresponding to the actual physical distribution of mass along the line of sight, while Eq. (1) is the model assumed by an inference analysis pipeline, one finds  \n∆ t truth -∆ t inferred ∆ t inferred ≈ κ c ( θ E ) . (3)  \nIn cosmography, this time delay bias translates into an H 0 inference bias  \nδ H = H inferred 0 -H truth 0 H inferred 0 ≈ κ c ( θ E ) . (4)  \nWe assume small κ c ( θ E ) ≪ 1 (our main discussion revolves around κ c ( θ E ) of order 0.1). It follows that a physical density core component that is not accounted for in the analysis pipeline would bias the inferred value of H 0 positively, that is, above the truth value of H 0 . Conversely,  \nan external prior on H 0 can allow cosmography analyses to detect density cores [53] down to fractional column densities of order ten percent, that are difficult to detect from kinematics.  \nEq. (4) equates the fractional H 0 bias to the core convergence at the Einstein angle. This may seem confusing, because κ c ( θ E ) is the ULDM convergence, not the fractional contribution of ULDM to the total convergence. In practice, however, the fractional H 0 bias is indeed proportional to fractional, rather than absolute, ULDM convergence: this is because the observable effect comes from the angular region containing imaging data, which is concentrated around the effective Einstein angle θ E . Focusing on spherical systems for simplicity, the Einstein angle of the null model is given by the solution of  \n2 ∫ 1 0 dy y κ 0 ( y θ E ) = 1 . (5)  \nThe integral is typically dominated by the region y = O (1) , implying κ 0 ( θ E ) = O (1) . For example, 3D density ρ 0 ( r ) ∝ r -γ leads to κ 0 ( θ E ) = (3 -γ ) / 2 , ranging between 1 and 0.5 for γ between 1 and 2. (The Navarro-Frenk-White (NFW [69]) halo falls in between γ = 1 and γ = 2 .) This means that up to a model-dependent O (1) factor that depends on the details of κ 0 ( θ ) , one actually has κ c ( θ E ) ∼ κ c ( θ E ) /κ 0 ( θ E ) , arising implicitly because lensing prescribes that we evaluate κ c ( θ ) specifically near θ E , determined in turn by κ 0 ( θ ) . Thus, up to a model-dependent O (1) factor, the fractional H 0 bias depends on the fractional column density of ULDM compared with the dominant column density of the halo near the projected Einstein angle.  \nImaging errors: Because the ULDM-induced convergence profile κ c ( θ ) is not exactly constant, the MSD is only approximate and imaging errors are produced. The magnitude of these errors can be roughly estimated by assuming spherical symmetry, and evaluating the change induced on θ E . The exact MSD would have the same value of θ E , satisfying α 0 ( θ E ) = θ E , regardless of the value of κ c . For finite core we parameterize the potentially-observable imaging distortion by the fractional change  \nα ( θ E ) := θ E (1 + δ E ) . (6)  \nFor this definition one finds [1]  \nδ E = 2 ∫ 1 0 dy y ( κ c ( y θ E ) -κ c ( θ E )) . (7)  \nMock data analysis [1, 54] suggests that current data for TDCOSMO systems may allow access to δ E ≳ 3% or so, although systematic modeling degeneracies may make this a difficult task (see, e.g. [70, 71]).  \nOur estimates for the H 0 bias δ H (Eq. (4)) and the imaging error δ E (Eq. (7)) are crude approximations because they are obtained in the spherical limit, while realistic systems exhibit deviation from sphericity, and also because they refer to the Einstein radius θ E while TDCOSMO systems have some imaging information also slightly off  \nof θ E . We explore this issue to some extent when we present our results, where, in addition to δ E , we report the maximal variation of κ c ( ⃗ θ ) over all | ⃗ θ | ≲ θ E . We also consider related effects, encoded in shear and flexion terms, in App. B. A more precise treatment of the lensing inference pipeline, likely using mock analyses along the lines of Ref. [1, 54], would become important once more comprehensive cosmological simulations including ULDM, CDM, and stars become available.', 'A. An aside: multiple-source systems': "The core-MSD effect could also, in principle, be detected via angular diameter distance ratios in multiple-source systems [70]. For N planes (ignoring lens-lens coupling for simplicity), the lens equation associated to the i th source is modified as  \n⃗ β i = ⃗ θ -C i ⃗ α ( ⃗ θ ) , C i := D S D LS i D S i D LS , (8)  \nwhere D S , L , LS are the angular diameter distances to source, lens, and lens-source, respectively, and z S , z L are the redshifts. A core-MSD term in the main lens is observationally equivalent to an effective rescaling of the factors  \nC ' i ≈ C i (1 + δκ S 1 i -δκ LS 1 i + κ c (1 -C i )) , (9)  \nwhere δκ S , LS 1 i are differential external convergence terms, and the equality holds at first order in external convergence and κ c . If one can constrain (e.g. from theory via cosmological simulations, or perhaps by combination with weak lensing data) the differential external convergence terms to be much smaller than the core factor κ c (1 -C i ) , and the redshifts and cosmological model parameters are measured well enough to determine the factors C i precisely (note that H 0 is not required for this purpose, because H 0 cancels out in C i ), then the core-MSD effect can be inferred from the difference between C ' i and C i . Preliminary analysis [70] suggests that this method may be useful for κ c ≳ 0 . 1 or so.", 'III. SIMULATION SETUP': "We consider 10 -25 ≤ m ≤ 10 -24 eV, and simulate galaxy halos with halo mass M halo 200 varying around 10 13 M ⊙ , and cosmological ULDM fraction M ULDM ≈ α χ M halo 200 , fixed to α χ = 0 . 15 in most of the analysis.  \nWe use a pseudo-spectral Schrodinger-Poisson solver, based on the method of Ref. [24], and implemented on a discretized grid with box length L ∼ 1 Mpc varying slightly between different simulations. This box length is chosen as compromise between our computational resources and the need to resolve the wave dynamics of ULDM in M halo 200 ≈ 10 13 M ⊙ halos, which have R 200 of a few hundred kpc. More details on the code are deferred to App. A.  \nSince ULDM makes-up only a subdominant fraction of the total DM, we need to implement also the dominant CDM component, as well as the stellar mass component in the halo (in TDCOSMO-like lens galaxies, the stellar mass is comparable to the DM mass near the projected Einstein radius). Our simulations cannot yet handle cosmological initial conditions, nor approximate realistic stellar physics. We therefore explore the sensitivity of the results to different approximate prescriptions for treating the galaxy halo:  \n- A. CDM and stars as static external potential: Our main results are obtained in a simple set-up, in which we approximate the dominant CDM and stellar mass distributions of the lens galaxy by fixed static external potentials. In this case, the numerical simulation includes only the ULDM dynamics, subject to the influence of the external potential.  \nWe approximate the CDM component, with mass M CDM = (1 -α χ ) M halo 200 , as an NFW potential. The NFW scale radius parameter r s in simulations is taken from the mass-concentration relations of [72]. The stellar component is taken to follow a Hernquist potential [73].  \n- B. CDM as a heavier field with m heavy ≫ m : Here we approximately implement CDM by adding to the code another dynamical ULDM field, with a larger particle mass m heavy ≫ m ; for simplicity, we continue to refer to this component as 'CDM'. The two fields interact with each other only via gravity. The approximation of CDM by a largem Schrodinger field relies on the Vlasov-Poisson correspondence [74]. We maintain sufficient numerical resolution to capture the de Broglie wavelength of the more massive field [75, 76]. The highest mass we can afford to simulate is m heavy = 10 -23 eV. Hence the ratio between m heavy and the lighter 'proper ULDM' m in our analysis varies between 10 ≤ m heavy /m ≤ 100 . In these simulations, the stellar component is neglected. Thus these simulations cannot be used directly to study gravitational lensing results. The goal of this analysis is to test the effect of backreaction of the subdominant ULDM field on the dominant CDM.  \nSimulating the background halo CDM with a dynamical field allows us to also test the impact of different initial conditions, notably the effect of galaxy mergers.  \nIn all simulations, the ULDM density distribution is initialized approximately as an α χ -scaled copy of the dominant CDM density profile. Initial conditions are defined via the Eddington procedure [31, 77], adapted to wave equations:  \nψ i ( x ) = (∆ v ) 3 / 2 ∑ ⃗v √ f ( E ( x, v ))e i m i ⃗x · ⃗v +i φ ⃗v , (10)  \nwhere ψ i is the Schrodinger field associated to m i (the index i refers to the proper ULDM field m and also to the  \nFIG. 1: Σ crit ( z L , z S ) in the ( z L , z S ) plane. TDCOSMO systems [50] are marked by dots.  \n<!-- image -->  \nm heavy implementation of 'CDM', where relevant), φ ⃗v is the random phase associated with the Fourier component of velocity ⃗v = ⃗ k/m i , ∆ v is the velocity resolution step available in the simulation, and  \nf ( E ) = 2 √ 8 π 2 ∫ √ E 0 d Q d 2 ρ dΨ 2 ( Q ) , (11)  \nE = Ψ( r ) -v 2 2 , Ψ = -Φ+Φ( r max ) , (12)  \nwhere ρ is the mass density of the simulated field, Φ is the total gravitational potential and Q = √ E Ψ . The random phases φ ⃗v cause interference patterns that make runs with different random seeds differ from each other on scales smaller than λ dB , even if they are designed to give the same density profile on average.  \nUnless stated otherwise, simulations begin at t = 0 and final results are read out at the Hubble time t = t H = 1 /H 0 .", 'IV. RESULTS': 'Our simulations produce results in terms of mass density profiles. However, gravitational lensing depends not only on the intrinsic density profile of the lens galaxy, but also on projection effects related to the observer-lens-source configuration. In particular, mapping mass column density into lensing convergence requires the critical density Σ crit = D S / (4 πGD LS D L ) . Fig. 1 shows Σ crit in the ( z L , z S ) plane (with the constraint z S > z L +0 . 3 ). The seven lens systems from TDCOSMO [50] are marked by white dots. Their Σ crit ranges from 1.8 to 3.4, with a mean of 2.4 × 10 9 M ⊙ / kpc 2 . In presenting results that require convergence, we either use a reference value of Σ crit = 2 × 10 9 M ⊙ / kpc 2 , or refer to specific examples from the TDCOSMO sample.', 'A. Static external CDM and stars': "We begin with simulations in which the dominant CDM and stellar components of the lens galaxy are approximated using static external NFW+Hernquist potentials. The dynamical ULDM has a cosmological fraction of α χ = 0 . 15 . The CDM NFW profile has M halo 200 ≈ 10 13 M ⊙ with scale radius r s ≈ 70 kpc. The stellar Hernquist profile has a total mass of ≈ 7 . 5 × 10 11 M ⊙ and scale radius ≈ 7 kpc. At r = 10 kpc the enclosed stellar mass is roughly twice the enclosed DM mass, compatible with estimates for typical TDCOSMO systems.  \nFig. 2 shows the radially-averaged 3D mass density for m = 2 × 10 -25 eV . We show H 0 bias δ H (top number  \nFIG. 2: Radially-averaged density. Here m = 2 × 10 -25 eV , α χ = 0 . 15 . Background halo implemented as static NFW+Hernquist potentials. We show H 0 bias δ H (top number in blue) and imaging error δ E (bottom number in red) for three redshift configurations ( z L , z S ) , corresponding to real TDCOSMO systems. Each redshift configuration maps to slightly different projected Einstein radius D L θ E . Colormap inset shows ULDM convergence κ c using Σ crit = 2 × 10 9 M ⊙ / kpc 2 (grid units: kpc). Yellow spot marks the lens halo center of mass.  \n<!-- image -->  \nin blue) and imaging error δ E (bottom number in red) for three redshift configurations ( z S , z L ) , corresponding to TDCOSMO systems 1 . Each redshift configuration maps to different projected Einstein radius D L θ E , indicated on the plot by thin vertical lines connected to the corresponding δ H and δ E . Colormap inset shows ULDM convergence κ c , using a typical reference value of Σ crit = 2 × 10 9 M ⊙ / kpc 2 to translate column density to lensing convergence (grid units: kpc). Yellow spot marks the halo center of mass.  \nFig. 2 showcases the AxionH0graphy scenario. ULDM is  \npredicted to produce cosmographic H 0 bias of a few percent, that would be essentially invisible in imaging and very challenging for stellar kinematics. Since the ULDM core falls off and joins the dominant CDM halo morphology after a few tens of kpc, it would also be difficult to detect in weak lensing, while cosmological constraints allow α χ = 0 . 15 . Thus, if cosmography delivers ≲ 2% precision in inferred H 0 [68] when the MSD is not modeled explicitly, then Fig. 2 demonstrates a physically-motivated scenario in which an ultralight axion may be first discovered from H 0 bias.  \nFig. 3 extends the results to m = (1 , 2 , 5 , 10) × 10 -25 eV (this contains again the result of Fig. 2 for comparison).  \nFor m = 10 -25 eV (top-left panel) the core extends far beyond the Einstein radius, making the MSD very precise. The ULDM-induced convergence is small - about a factor of 5 smaller than α χ - causing a small H 0 bias δ H ≲ 3% that is still, however, somewhat larger than cosmography precision expected in forecasts. It is interesting to also investigate the regime m < 10 -25 eV, but our simulations do not extend there because of our technical limitation in simulating a sufficiently large box.  \nm ≈ 2 × 10 -25 eV (top-right panel, or Fig. 2) is a 'sweetspot' for AxionH0graphy with H 0 bias δ H ∼ 10% .  \nFor m = 5 × 10 -25 eV (bottom-left) the spatial scale of interference patterns in the ULDM field; the displacement of the ULDM core from the center of mass of the dominant CDM+stars halo [75, 78]; and the core radius itself; which are all of the order of λ dB , approach the scale of the projected Einstein radius for a TDCOSMO-like galaxy. This causes deviations from the MSD approximation, seen in the bottom-left panel of Fig. 3 with δ E ≳ 4% . In this regime the spherical approximation used in Eqs. (4) and (7) becomes less accurate. Mock data analysis is needed to investigate this regime quantitatively. As a tentative rough estimate of the impact on imaging, noting that Eq. (7) is the spherical average of the difference between κ c ( ⃗ θ ) and its value at θ E , we can compare δ E to the maximal variation in κ c ( ⃗ θ ) with respect to κ c ( θ E ) , taken across the region | ⃗ θ | ≲ θ E :  \n∆ κ c := max | ⃗ θ | <θ E ( | κ c ( ⃗ θ ) -κ c ( θ E ) | ) , (13)  \nwhere the coordinates are centered on the CDM halo center of mass. Choosing D L θ E ≈ 10 kpc for the boundary of the region, and using Σ crit = 2 × 10 9 M ⊙ / kpc for reference, we find ∆ κ c ≈ 0 . 09 for the m = 5 × 10 -25 eV simulation. Thus ULDM interference patterns and core displacement cause imaging errors of comparable magnitude to those caused by the finite-core effect.  \nFor m = 10 -24 eV (bottom-right), the core radius is smaller than D L θ E and the convergence map shows that the MSD approximation is poor, violating Eqs. (4) and (7) at the O (1) level.", 'B. Dynamical CDM via m heavy ≫ m': "The middle panel of Fig. 4 shows simulation results for m = 5 × 10 -25 eV, where the CDM component of the halo  \nFIG. 3: Same set-up as Fig. 2 (reproduced in the top-right panel), extended to m = (1 , 5 , 10) × 10 -25 eV (top-left, bottom-left, and bottom-right panels, respectively).  \n<!-- image -->  \nis simulated by a dynamical field with m heavy = 10 -23 eV. Initial conditions are shown by dashed lines. Final time results are shown by thick bands that contain the extreme values of the density field during the time t = t H ± 1 Gyr. The thickness of these bands demonstrates the amplitude of core oscillations and dynamical interference patterns. (Although we did not highlight this fact, comparable core oscillations exist in the simulations of Fig. 3.)  \nFor reference, the left panel of Fig. 4 shows results obtained with a fixed static NFW background, which approximately matches the CDM initial conditions in the dynamical case. These simulations do not include a stellar mass component: the goal is to investigate the impact of the dynamical background on the ULDM core in a simple toy model.  \nTo test the sensitivity of our results to initial conditions, notably major mergers, we also perform runs initialized as two NFW halos of similar initial mass profiles that merge dynamically. The two halos are initially separated by ≈ 100 kpc , at rest w.r.t. each other. The result is shown in the right panel of Fig. 4.  \nFig. 4 suggests that ULDM cores are not very sensitive to the details of CDM backreaction or halo initial conditions. In all three panels, the ULDM core central density is similar to within about a factor of two, with a similar core radius of  \nabout 20 kpc. Note that our simulations can only roughly resolve the detailed dynamics of the m heavy = 10 -23 eV 'CDM' component, which exhibits intereference patters with typical wavelength of order 1 kpc. In particular, the right panel in Fig. 4 shows that the m heavy field is forming a soliton of its own in the inner few kpc region of the halo. This could be interesting as it is expected that the lighter ULDM species, despite being subdominant in the total mass budget, could accelerate dynamical relaxation in the dominant CDM component [42]. For AxionH0graphy, however, it is sufficient to notice that the ULDM core is not strongly dependent on the details of the main CDM component.", 'C. Results summary for α χ = 0 . 15': "Tab. I summarizes the main results of our simulations for α χ = 0 . 15 . For each of m = (1 , 2 , 5 , 10) × 10 -25 eV, we simulate a set of scenarios: single halo evolution with static CDM NFW potential; single halo with static CDM+stars; single halo with dynamical CDM implemented via m heavy =  \n<!-- image -->  \n<!-- image -->  \nFIG. 4: Simulation results for m = 5 × 10 -25 eV , α χ = 0 . 15 . Panels differ in the implementation of the CDM halo component, and in initial conditions. Left panel: CDM implemented as static NFW potential. Middle panel: CDM implemented by another dynamical field with m heavy = 10 -23 eV. Right panel: CDM implemented as m heavy = 10 -23 eV, initialized as two halos that merge dynamically. Final time results are shown by thick bands that contain the extreme values of the density field during the time t = t H ± 1 Gyr, capturing the effect of core oscillations. Initial conditions are shown by dashed lines.  \n<!-- image -->  \nTABLE I: Summary of simulation results for α χ = 0 . 15 . From left to right: (1) ULDM particle mass, (2) M 200 (final snapshot), (3) soliton mass from Eq. (14), (4) core radius, vis ρ ULDM ( r c ) = 0 . 5 ρ ULDM (0) , (5) half de-Broglie wavelength, using Eq. (16), (6) projected Einstein radius r E = D l θ E , (7) δ E from Eq. (7) (in %), (8) maximal convergence variation from Eq. (13) (in %), (9) H 0 bias, (10) simulated configuration.  \n10 -23 eV 2 ; two-halo merger with dynamical m heavy CDM; and three-halo merger with dynamical m heavy CDM. All simulations are run until t = t H . For each simulation we quote the CDM halo mass M 200 (extracted at t = t H when CDMis dynamical). We now discuss each additional column  \nof Tab. I in turn.  \nSoliton mass: Ref. [1] parametrized the ULDM soliton profile in the presence of external potential with the form  \nρ sol ( r ) = λ 4 (1 + a 2 λ 2 r 2 ) 2 b m 2 4 πG , (14)  \nwhere a, b are fit parameters, and λ a scaling variable. For simplicity, we use this parameterization even though the field configurations we find in the simulations deviate from  \nspherical symmetry. Fitting the core with Eq. (14), we find the soliton mass as  \nM sol = π 3 / 2 λ 4 Γ( -3 / 2 + 2 b ) a 3 λ 3 Γ(2 b ) . (15)  \nCore scale radius r c , and λ dB : We denote the scale length of the core in simulations by r c , defined as the radius at which the ULDM density decreases by a factor of 2 from it's peak value. r c is comparable with λ dB / 2 , estimated by  \nλ dB ≈ 20 kpc ( 180 kms -1 v )( 5 × 10 -25 eV m ) , (16)  \nwhere v is the 1D velocity dispersion  \nv ≈ √ GM 200 3 R 200 ≈ 180 kms -1 ( M 200 10 13 M ⊙ ) 1 / 3 , (17)  \nusing R 200 ≈ 450 kpc ( M 200 10 13 M ⊙ ) 1 / 3 . 3  \nEinstein radius r E : It is useful to compute the projected Einstein radius r E = D L θ E for comparison with r c . We expect precise MSD when r E ≪ r c . In simulations that do not include a stellar mass component, we assign a value for θ E by adding to the convergence computation a contribution from a Hernquist profile, scaling the stellar mass linearly with M 200 . The error bars we quote for r E are intended to represent projection effects. They relate to the smallest and largest θ E values obtained across the seven TDCOSMO systems (see Fig. 1, and compare also to the vertical thin lines in Figs. 2 and 3).  \nImaging error estimates δ E and ∆ κ c : We quote both δ E (Eq. (7)) and ∆ κ c (Eq. (13)) as estimates of imaging errors. Error bars refer to maximum and minimum values obtained from the redshift configurations of the TDCOSMO sample.  \nWe highlight in red δ E and ∆ κ c entries that exceed 0 . 03 , as these halos may exhibit detectable MSD-breaking imaging errors. Note again that our the error estimates do not take into account degeneracies involving external shear and flexion [70, 71]. We comment on this point in App. B.  \nH 0 bias: We use Eq. (4) to estimate H 0 bias from the radial-averaged convergence profile. The error bars we show come from the maximum and minimum values of κ c ( θ E ) obtained by the redshift configurations of TDCOSMO systems.  \nFinally, it is important to note that the ULDM mass in our simulated halos is matched to the cosmological fraction α χ . If cosmological evolution causes significant suppression of the ULDM fraction on distance scales of ∼ Mpc, the actual ULDM mass in the inner region of galaxies may be lower than α χ . Our results could then exaggerate the H 0 bias. This adds motivation to extend the analysis to full cosmological simulations in future work.  \nFIG. 5: Comparison of results for α χ = 0 . 01 (thick solid lines) and α χ = 0 . 15 (thin dot-dashed lines), done for m = 5 × 10 -25 eV, with dynamical CDM modeled via m heavy field.  \n<!-- image -->", 'D. Brief view of α χ < 0 . 15': 'Our work focused on α χ = 0 . 15 , which is roughly the maximal ULDM fraction consistent with cosmological data [12, 15]. Here we take a brief glance at smaller α χ .  \nFig. 5 compares simulations for α χ = 0 . 01 (thick solid lines) and α χ = 0 . 15 (thin dot-dashed lines), done for m = 5 × 10 -25 eV. Results for α χ = 0 . 15 are the same as in the middle panel of Fig. 4. The comparison highlights the soliton core that begins to form in the α χ = 0 . 15 case, but not in the α χ = 0 . 01 case, in qualitative agreement with results in [27, 79].', 'V. TIMESCALE OF CORE FORMATION': "So far we focused on simulation snapshots at t = t H . We now comment on the time evolution leading to these snapshots.  \nFig. 6 shows the time evolution of the peak ULDM density of two simulations, containing m = 2 × 10 -25 eV (left panel) and 5 × 10 -25 eV (right panel) at α χ = 0 . 15 with dynamical m heavy CDM. We can identify three regimes in the dynamics of the core:  \n- 1. Wave pressure expels mass outwards, forming flattened region. As explained in Sec. III, in the (simplistic) initial conditions of our simulations the ULDM halo is an approximate copy of the dominant CDM NFW halo, scaled-down by a factor of α χ . The Eddington prescription used to initialize the field does not take into account wave mechanics pressure. This omission is analogous to neglecting two-body stellar encounters when assigning initial conditions in standard N-body simulations. Once we start the run, ULDM in the dense central region of the halo evolves quickly due to wave pressure to form a flattened core,  \n<!-- image -->  \nFIG. 6: ULDM core density as a function of time. Left (right) panel: m = 2 (5) × 10 -25 eV.  \n<!-- image -->  \non a wave-mechanics time scale t dB given by  \nt dB = 1 mv 2 ≈ 0 . 1 ( 200 km / s v ) 2 ( 5 × 10 -25 eV m ) Gyr . (18)  \nt dB is comparable but probably somewhat shorter than the ( m -independent) free-fall time,  \nt freefall = √ 3 πR 3 200 32 GM 200 ≈ 0 . 8 Gyr . (19)  \nWe do not investigate the details of the initial shorttime dynamics, because these necessarily depend on the details of our rather ad-hoc initial conditions. The only comment we make in this regard is that the later evolution of the core is relatively independent on the details of the initial conditions, as those are quickly 'ironed out' by the wave pressure.  \n- 2. Dynamical relaxation and emergence of a soliton. After the initial density drop, the core density fluctuates for a while without substantial mean growth.  \nThe state of the field during this stage may not be very different from the random noise state used in Ref. [24, 37, 38] to investigate soliton condensation in a statistically homogeneous background. The relaxation timescale in that case was found to be  \nt rel ≈ √ 2 12 π 3 m 3 v 6 G 2 ρ 2 χ Λ , (20)  \nwhere ρ χ refers to the ULDM density and Λ is a Coulomb logarithm, of order unity near the core. The simulations of soliton condensation from random noise, done in Ref. [24], show self-gravitating solitons emerging from noise around t ≈ t rel , after which time the soliton mass grows approximately as M sol ∝ t 1 2 .  \nTo compare this expectation to our simulations, we can read ρ χ directly from Fig. 6. However, t rel depends very sensitively on the local velocity dispersion v in Eq. (20), which is not a straightforward quantity to define precisely in our situation, where dispersion sourced by the background potential adds to the motion of the ULDM core around the halo center of mass. We make a crude estimate as follows. From the simulation, we read the core radius r c of the ULDM; we find r c ≈ 26 (18) kpc for m = 2 (5) × 10 -25 eV, with enclosed background halo mass M halo ( r c ) ≈ 0 . 35 (0 . 2) × 10 12 M ⊙ . We then estimate the circular velocity at r c as v 2 circ ( r c ) = GM halo ( r c ) /r c , and use v ≈ √ 3 × v circ in Eq. (20). This crude estimate gives v circ ( r c ) ≈ 245 (222) km/s, and t rel ≈ 0 . 9 (8) Gyr, for m = 2 (5) × 10 -25 eV, corresponding to the left (right) panel of Fig. 6. This estimate is not very far from the behavior seen in Fig. 6. However, to appreciate how sensitive these estimates are to numerical details, note that using r c → 2 × r c in our estimate would multiply the resulting t rel by a factor of ∼ 8 . We can only conclude that our results in Fig. 6 are not clearly inconsistent with Eq. (20).  \n- 3. Soliton growth. In the third stage of the evolution the central density in our simulations shows a clear trend of growth, consistent with ρ c ∝ t 1 2 . To guide the eye, we superimpose a red line with this slope in both panels.  \nNote that we see ρ c ∝ t 1 2 , while Ref. [24] reported that the mass of the soliton, M sol , grows as ∝ t 1 2 . This scaling is consistent if the soliton is not selfgravitating, but rather dominated by the velocity dispersion due to the background potential, in which case the soliton radius only weakly depends on the soliton mass and one finds M sol ∝ ρ c . This is indeed the relevant regime in our simulations.  \nWe track the build-up of the soliton to t ≫ t H .  \nWe do not see a clear break of the soliton growth, although we cannot exclude the beginning of a break at late time. Saturation of soliton growth is expected to occur when the velocity dispersion due to the soliton's self-gravity becomes comparable to the mean velocity dispersion in the halo (see discussion and references in [41, 42]). The solitons in our simulations do not reach this point at t = t H ; even the relatively prominent ULDM cores seen in the middle and right panels of Fig. 4, where the ULDM density is comparable to the CDM density near the core radius, still induce a contribution to v circ that is smaller by a factor of few compared to the CDMdominated peak halo v circ . This result is consistent with expectations in Ref. [1] that identified soliton growth as a bottleneck for AxionH0graphy.", 'VI. SUMMARY': "Even a subdominant contribution of ultralight dark matter (ULDM), of order ten percent of the total DM budget, can affect strong gravitational lensing at a level accessible to current observations and makes a promising target for near-future surveys. This happens because ULDM condenses into cores around massive galaxies. We reported results of a series of numerical simulations investigating this effect, extending the original analysis of Ref. [1].  \nFor ULDM particle mass m ∼ 10 -25 eV, the scale size of the cores is a few tens of kpc in massive elliptical galaxies of the type used for cosmographic determination of H 0 . The main effect of a subdominant component of ULDM would be to cause a subtle positive bias in the inferred value of H 0 . Adding an external prior on H 0 from another data set (e.g. CMB, or the supernovae distance ladder), the lensing H 0  \n- [1] K. Blum and L. Teodori, 'Gravitational lensing H0 tension from ultralight axion galactic cores,' Phys. Rev. D 104 no. 12, (2021) 123011, arXiv:2105.10873 [astro-ph.CO] .\n- [2] J. Preskill, M. B. Wise, and F. Wilczek, 'Cosmology of the Invisible Axion,' Phys. Lett. B 120 (1983) 127-132.\n- [3] L. F. Abbott and P. Sikivie, 'A Cosmological Bound on the Invisible Axion,' Phys. Lett. B 120 (1983) 133-136.\n- [4] M. Dine and W. Fischler, 'The Not So Harmless Axion,' Phys. Lett. B 120 (1983) 137-141.\n- [5] W. Hu, R. Barkana, and A. Gruzinov, 'Cold and fuzzy dark matter,' Phys. Rev. Lett. 85 (2000) 1158-1161, arXiv:astro-ph/0003365 [astro-ph] .\n- [6] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, 'Ultralight scalars as cosmological dark matter,' Phys. Rev. D95 no. 4, (2017) 043541, arXiv:1610.08297 [astro-ph.CO] .\n- [7] L. Hui, 'Wave Dark Matter,' Ann. Rev. Astron. Astrophys. 59 (2021) 247-289, arXiv:2101.11735 [astro-ph.CO] .\n- [8] E. G. M. Ferreira, 'Ultra-light dark matter,' Astron.  \nbias can be detected. We dub this method AxionH0graphy. It provides concrete theoretical motivation for combining cosmography with an external H 0 prior to investigate the structure of massive galaxies, demonstrating the potential for fundamental discovery.  \n- m ≳ 5 × 10 -25 eV ULDM would introduce imaging distortions in addition, if not instead of the H 0 bias. This regime for m is also interesting but the emphasis shifts to detailed lens modeling, rather than an H 0 prior.  \nWe believe our simulations provide robust results regarding the formation, general scale, and dynamics of ULDM cores. However, the set-up we explore is still simplistic in important aspects: first, we cannot yet simulate fully self-consistently the dominant component of cold DM (CDM) and stars, that are both important for strong lensing in massive elliptical galaxies. Our simulations either use static external potentials to replace these dominant backgrounds, or a wave implementation of CDM that cannot resolve the dynamics below ∼ 1 kpc scales. Second, our simulations start from rather ad-hoc initial conditions, and may miss some features of cosmological structure formation. With these caveats in mind, our predictions for the amplitude of the H 0 bias must be understood as crude approximations, motivating more comprehensive analysis.", 'Acknowledgments': "We are grateful to E. Hardy and M. Gorghetto for useful discussion and guidance in the use of the pseudo-spectral solver, originally developed by E. Hardy. Different versions of the code used for this work were also used in [80-82]. We thank M. Nori and M. Baldi for collaboration in the early stage of this work. KB and LT were supported by ISF grant 1784/20 and by MINERVA grant 714123  \n- Astrophys. Rev. 29 no. 1, (2021) 7, arXiv:2005.03254 [astro-ph.CO] .\n- [9] P. Svrcek and E. Witten, 'Axions In String Theory,' JHEP 06 (2006) 051, arXiv:hep-th/0605206 [hep-th] .\n- [10] A. Arvanitaki, S. Dimopoulos, S. Dubovsky, N. Kaloper, and J. March-Russell, 'String Axiverse,' Phys. Rev. D81 (2010) 123530, arXiv:0905.4720 [hep-th] .\n- [11] D. J. E. Marsh, 'Axion Cosmology,' Phys. Rept. 643 (2016) 1-79, arXiv:1510.07633 [astro-ph.CO] .\n- [12] A. Lague, J. R. Bond, R. Hloˇzek, K. K. Rogers, D. J. E. Marsh, and D. Grin, 'Constraining ultralight axions with galaxy surveys,' JCAP 01 no. 01, (2022) 049, arXiv:2104.07802 [astro-ph.CO] .\n- [13] V. Irˇsiˇc, M. Viel, M. G. Haehnelt, J. S. Bolton, and G. D. Becker, 'First constraints on fuzzy dark matter from Lymanα forest data and hydrodynamical simulations,' Phys. Rev. Lett. 119 no. 3, (2017) 031302, arXiv:1703.04683 [astro-ph.CO] .\n- [14] E. Armengaud, N. Palanque-Delabrouille, D. J. E. Marsh, J. Baur, and C. Y¨ı¿œche, 'Constraining the mass of light bosonic dark matter using SDSS Lymanα forest,' Mon.  \nNot. Roy. Astron. Soc. 471 no. 4, (2017) 4606-4614, arXiv:1703.09126 [astro-ph.CO] .  \n- [15] T. Kobayashi, R. Murgia, A. De Simone, V. Irˇsiˇc, and M. Viel, 'Lymanα constraints on ultralight scalar dark matter: Implications for the early and late universe,' Phys. Rev. D 96 no. 12, (2017) 123514, arXiv:1708.00015 [astro-ph.CO] .\n- [16] K.-H. Leong, H.-Y. Schive, U.-H. Zhang, and T. Chiueh, 'Testing extreme-axion wave dark matter using the BOSS Lyman-Alpha forest data,' Mon. Not. Roy. Astron. Soc. 484 (2019) 4273, arXiv:1810.05930 [astro-ph.CO] .\n- [17] F. S. Guzman and L. A. Urena-Lopez, 'Evolution of the Schrodinger-Newton system for a selfgravitating scalar field,' Phys. Rev. D69 (2004) 124033, arXiv:gr-qc/0404014 [gr-qc] .\n- [18] H.-Y. Schive, T. Chiueh, and T. Broadhurst, 'Cosmic Structure as the Quantum Interference of a Coherent Dark Wave,' Nature Phys. 10 (2014) 496-499, arXiv:1406.6586 [astro-ph.GA] .\n- [19] H.-Y. Schive, M.-H. Liao, T.-P. Woo, S.-K. Wong, T. Chiueh, T. Broadhurst, and W. Y. P. Hwang, 'Understanding the Core-Halo Relation of Quantum Wave Dark Matter from 3D Simulations,' Phys. Rev. Lett. 113 no. 26, (2014) 261302, arXiv:1407.7762 [astro-ph.GA] .\n- [20] B. Schwabe, J. C. Niemeyer, and J. F. Engels, 'Simulations of solitonic core mergers in ultralight axion dark matter cosmologies,' Phys. Rev. D94 no. 4, (2016) 043513, arXiv:1606.05151 [astro-ph.CO] .\n- [21] J. Veltmaat and J. C. Niemeyer, 'Cosmological particle-in-cell simulations with ultralight axion dark matter,' Phys. Rev. D94 no. 12, (2016) 123523, arXiv:1608.00802 [astro-ph.CO] .\n- [22] P. Mocz, M. Vogelsberger, V. H. Robles, J. Zavala, M. Boylan-Kolchin, A. Fialkov, and L. Hernquist, 'Galaxy formation with BECDM: I. Turbulence and relaxation of idealized haloes,' Mon. Not. Roy. Astron. Soc. 471 no. 4, (2017) 4559-4570, arXiv:1705.05845 [astro-ph.CO] .\n- [23] J. Veltmaat, J. C. Niemeyer, and B. Schwabe, 'Formation and structure of ultralight bosonic dark matter halos,' Phys. Rev. D98 no. 4, (2018) 043509, arXiv:1804.09647 [astro-ph.CO] .\n- [24] D. G. Levkov, A. G. Panin, and I. I. Tkachev, 'Gravitational Bose-Einstein condensation in the kinetic regime,' Phys. Rev. Lett. 121 no. 15, (2018) 151301, arXiv:1804.05857 [astro-ph.CO] .\n- [25] B. Eggemeier and J. C. Niemeyer, 'Formation and mass growth of axion stars in axion miniclusters,' Phys. Rev. D 100 no. 6, (2019) 063528, arXiv:1906.01348 [astro-ph.CO] .\n- [26] J. Chen, X. Du, E. W. Lentz, D. J. E. Marsh, and J. C. Niemeyer, 'New insights into the formation and growth of boson stars in dark matter halos,' arXiv:2011.01333 [astro-ph.CO] .\n- [27] B. Schwabe, M. Gosenca, C. Behrens, J. C. Niemeyer, and R. Easther, 'Simulating mixed fuzzy and cold dark matter,' Phys. Rev. D 102 no. 8, (2020) 083518, arXiv:2007.08256 [astro-ph.CO] .\n- [28] J. Zhang, H. Liu, and M.-C. Chu, 'Cosmological Simulation for Fuzzy Dark Matter Model,' Front. Astron. Space Sci. 5 (2019) 48, arXiv:1809.09848 [astro-ph.CO] .\n- [29] N. C. Amorisco and A. Loeb, 'First constraints on Fuzzy Dark Matter from the dynamics of stellar streams in the  \n- Milky Way,' arXiv:1808.00464 [astro-ph.GA] .\n- [30] B. V. Church, J. P. Ostriker, and P. Mocz, 'Heating of Milky Way disc Stars by Dark Matter Fluctuations in Cold Dark Matter and Fuzzy Dark Matter Paradigms,' arXiv:1809.04744 [astro-ph.GA] .\n- [31] L. Lancaster, C. Giovanetti, P. Mocz, Y. Kahn, M. Lisanti, and D. N. Spergel, 'Dynamical Friction in a Fuzzy Dark Matter Universe,' JCAP 01 (2020) 001, arXiv:1909.06381 [astro-ph.CO] .\n- [32] N. Bar, D. Blas, K. Blum, and H. Kim, 'Assessing the Fornax globular cluster timing problem in different models of dark matter,' Phys. Rev. D 104 no. 4, (2021) 043021, arXiv:2102.11522 [astro-ph.GA] .\n- [33] D. Dutta Chowdhury, F. C. van den Bosch, P. van Dokkum, V. H. Robles, H.-Y. Schive, and T. Chiueh, 'On the Dynamical Heating of Dwarf Galaxies in a Fuzzy Dark Matter Halo,' Astrophys. J. 949 no. 2, (2023) 68, arXiv:2303.08846 [astro-ph.GA] .\n- [34] P.-H. Chavanis, 'Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: I. Analytical results,' Phys. Rev. D84 (2011) 043531, arXiv:1103.2050 [astro-ph.CO] .\n- [35] P. H. Chavanis and L. Delfini, 'Mass-radius relation of Newtonian self-gravitating Bose-Einstein condensates with short-range interactions: II. Numerical results,' Phys. Rev. D84 (2011) 043532, arXiv:1103.2054 [astro-ph.CO] .\n- [36] A. H. Guth, M. P. Hertzberg, and C. Prescod-Weinstein, 'Do Dark Matter Axions Form a Condensate with Long-Range Correlation?,' Phys. Rev. D 92 no. 10, (2015) 103513, arXiv:1412.5930 [astro-ph.CO] .\n- [37] B. Bar-Or, J.-B. Fouvry, and S. Tremaine, 'Relaxation in a Fuzzy Dark Matter Halo,' Astrophys. J. 871 (2019) 28, arXiv:1809.07673 [astro-ph.GA] .\n- [38] B. Bar-Or, J.-B. Fouvry, and S. Tremaine, 'Relaxation in a Fuzzy Dark Matter Halo. II. Self-consistent Kinetic Equations,' Astrophys. J. 915 no. 1, (2021) 27, arXiv:2010.10212 [astro-ph.GA] .\n- [39] J. H.-H. Chan, S. Sibiryakov, and W. Xue, 'Condensation and evaporation of boson stars,' JHEP 01 (2024) 071, arXiv:2207.04057 [astro-ph.CO] .\n- [40] J. H.-H. Chan, S. Sibiryakov, and W. Xue, 'Boson star normal modes,' JHEP 08 (2023) 045, arXiv:2304.13054 [astro-ph.CO] .\n- [41] N. Bar, D. Blas, K. Blum, and S. Sibiryakov, 'Galactic rotation curves versus ultralight dark matter: Implications of the soliton-host halo relation,' Phys. Rev. D98 no. 8, (2018) 083027, arXiv:1805.00122 [astro-ph.CO] .\n- [42] N. Bar, K. Blum, and C. Sun, 'Galactic rotation curves versus ultralight dark matter: A systematic comparison with SPARC data,' Phys. Rev. D 105 no. 8, (2022) 083015, arXiv:2111.03070 [hep-ph] .\n- [43] T. Bernal, L. M. Fern'andez-Hern'andez, T. Matos, and M. A. Rodr'ıguez-Meza, 'Rotation curves of high-resolution LSB and SPARC galaxies with fuzzy and multistate (ultralight boson) scalar field dark matter,' Mon. Not. Roy. Astron. Soc. 475 no. 2, (2018) 1447-1468, arXiv:1701.00912 [astro-ph.GA] .\n- [44] L. V. E. Koopmans, 'Gravitational-mass imaging of CDM substructure,' Mon. Not. Roy. Astron. Soc. 363 (2005) 1136, arXiv:astro-ph/0501324 .\n- [45] D. D. Xu, D. Sluse, L. Gao, J. Wang, C. Frenk, S. Mao, and P. Schneider, 'How well can cold-dark-matter substructures account for the observed lensing flux-ratio anomalies?,' 7, 2013. arXiv:1307.4220 [astro-ph.CO] .\n- [46] S. Vegetti and L. V. E. Koopmans, 'Bayesian Strong Gravitational-Lens Modelling on Adaptive Grids: Objective Detection of Mass Substructure in Galaxies,' Mon. Not. Roy. Astron. Soc. 392 (2009) 945, arXiv:0805.0201 [astro-ph] .\n- [47] S. Vegetti et al. , 'Strong gravitational lensing as a probe of dark matter,' arXiv:2306.11781 [astro-ph.CO] .\n- [48] A. Laroche, D. Gilman, X. Li, J. Bovy, and X. Du, 'Quantum fluctuations masquerade as halos: Bounds on ultra-light dark matter from quadruply-imaged quasars,' arXiv:2206.11269 [astro-ph.CO] .\n- [49] D. M. Powell, S. Vegetti, J. P. McKean, S. D. M. White, E. G. M. Ferreira, S. May, and C. Spingola, 'A lensed radio jet at milli-arcsecond resolution II: Constraints on fuzzy dark matter from an extended gravitational arc,' arXiv:2302.10941 [astro-ph.CO] .\n- [50] M. Millon et al. , 'TDCOSMO. I. An exploration of systematic uncertainties in the inference of H 0 from time-delay cosmography,' Astron. Astrophys. 639 (2020) A101, arXiv:1912.08027 [astro-ph.CO] .\n- [51] E. E. Falco, M. V. Gorenstein, and I. I. Shapiro, 'On model-dependent bounds on H 0 from gravitational images : application to Q 0957+561 A, B.,' apjl 289 (Feb, 1985) L1-L4.\n- [52] P. Schneider and D. Sluse, 'Mass-sheet degeneracy, power-law models and external convergence: Impact on the determination of the Hubble constant from gravitational lensing,' Astron. Astrophys. 559 (2013) A37, arXiv:1306.0901 [astro-ph.CO] .\n- [53] K. Blum, E. Castorina, and M. Simonovi'c, 'Could Quasar Lensing Time Delays Hint to a Core Component in Halos, Instead of H 0 Tension?,' Astrophys. J. Lett. 892 no. 2, (2020) L27, arXiv:2001.07182 [astro-ph.CO] .\n- [54] S. Birrer et al. , 'TDCOSMO - IV. Hierarchical time-delay cosmography joint inference of the Hubble constant and galaxy density profiles,' Astron. Astrophys. 643 (2020) A165, arXiv:2007.02941 [astro-ph.CO] .\n- [55] A. J. Shajib et al. , 'TDCOSMO. XII. Improved Hubble constant measurement from lensing time delays using spatially resolved stellar kinematics of the lens galaxy,' Astron. Astrophys. 673 (2023) A9, arXiv:2301.02656 [astro-ph.CO] .\n- [56] L. Verde, T. Treu, and A. G. Riess, 'Tensions between the Early and the Late Universe,' Nature Astron. 3 (7, 2019) 891, arXiv:1907.10625 [astro-ph.CO] .\n- [57] E. Di Valentino, O. Mena, S. Pan, L. Visinelli, W. Yang, A. Melchiorri, D. F. Mota, A. G. Riess, and J. Silk, 'In the realm of the Hubble tension-a review of solutions,' Class. Quant. Grav. 38 no. 15, (2021) 153001, arXiv:2103.01183 [astro-ph.CO] .\n- [58] W. L. Freedman and B. F. Madore, 'Progress in direct measurements of the Hubble constant,' JCAP 11 (2023) 050, arXiv:2309.05618 [astro-ph.CO] .\n- [59] L. Verde, N. Schoneberg, and H. Gil-Mar'ın, 'A tale of many H 0 ,' arXiv:2311.13305 [astro-ph.CO] .\n- [60] DES Collaboration, T. Treu et al. , 'The STRong lensing Insights into the Dark Energy Survey (STRIDES) 2016 follow-up campaign I. Overview and classification of candidates selected by two techniques,' Mon. Not. Roy. Astron. Soc. 481 no. 1, (2018) 1041-1054, arXiv:1808.04838 [astro-ph.CO] .\n- [61] SKA Collaboration, D. J. Bacon et al. , 'Cosmology with Phase 1 of the Square Kilometre Array: Red Book 2018: Technical specifications and performance forecasts,' Publ.  \n- Astron. Soc. Austral. 37 (2020) e007, arXiv:1811.02743 [astro-ph.CO] .\n- [62] Euclid Collaboration, A. Blanchard et al. , 'Euclid preparation: VII. Forecast validation for Euclid cosmological probes,' Astron. Astrophys. 642 (2020) A191, arXiv:1910.09273 [astro-ph.CO] .\n- [63] LSST Science, LSST Project Collaboration, P. A. Abell et al. , 'LSST Science Book, Version 2.0,' arXiv:0912.0201 [astro-ph.IM] .\n- [64] M. Oguri and P. J. Marshall, 'Gravitationally lensed quasars and supernovae in future wide-field optical imaging surveys,' Mon. Not. Roy. Astron. Soc. 405 (2010) 2579-2593, arXiv:1001.2037 [astro-ph.CO] .\n- [65] B. M. Dobke, L. J. King, C. D. Fassnacht, and M. W. Auger, 'Estimating cosmological parameters from future gravitational lens surveys,' Mon. Not. Roy. Astron. Soc. 397 (2009) 311, arXiv:0904.1437 [astro-ph.CO] .\n- [66] R. Gavazzi, T. Treu, L. V. E. Koopmans, A. S. Bolton, L. A. Moustakas, S. Burles, and P. J. Marshall, 'The Sloan Lens ACS Survey. VI: Discovery and analysis of a double Einstein ring,' Astrophys. J. 677 (2008) 1046, arXiv:0801.1555 [astro-ph] .\n- [67] K. Liao et al. , 'Strong Lens Time Delay Challenge: II. Results of TDC1,' Astrophys. J. 800 no. 1, (2015) 11, arXiv:1409.1254 [astro-ph.IM] .\n- [68] S. Birrer and T. Treu, 'TDCOSMO - V. Strategies for precise and accurate measurements of the Hubble constant with strong lensing,' Astron. Astrophys. 649 (2021) A61, arXiv:2008.06157 [astro-ph.CO] .\n- [69] J. F. Navarro, C. S. Frenk, and S. D. M. White, 'A Universal density profile from hierarchical clustering,' Astrophys. J. 490 (1997) 493-508, arXiv:astro-ph/9611107 [astro-ph] .\n- [70] L. Teodori, K. Blum, E. Castorina, M. Simonovi'c, and Y. Soreq, 'Comments on the mass sheet degeneracy in cosmography analyses,' JCAP 07 no. 07, (2022) 027, arXiv:2201.05111 [astro-ph.CO] .\n- [71] L. Teodori and K. Blum, 'Host group degeneracy in gravitational lensing time delay determination of H 0 ,' JCAP 08 (2023) 065, arXiv:2305.19151 [astro-ph.CO] .\n- [72] A. V. Maccio', A. A. Dutton, and F. C. v. d. Bosch, 'Concentration, Spin and Shape of Dark Matter Haloes as a Function of the Cosmological Model: WMAP1, WMAP3 and WMAP5 results,' Mon. Not. Roy. Astron. Soc. 391 (2008) 1940-1954, arXiv:0805.1926 [astro-ph] .\n- [73] L. Hernquist, 'An Analytical Model for Spherical Galaxies and Bulges,' ApJ 356 (June, 1990) 359.\n- [74] P. Mocz, L. Lancaster, A. Fialkov, F. Becerra, and P.-H. Chavanis, 'Schrodinger-Poisson-Vlasov-Poisson correspondence,' Phys. Rev. D 97 no. 8, (2018) 083519, arXiv:1801.03507 [astro-ph.CO] .\n- [75] X. Li, L. Hui, and G. L. Bryan, 'Numerical and Perturbative Computations of the Fuzzy Dark Matter Model,' Phys. Rev. D 99 no. 6, (2019) 063509, arXiv:1810.01915 [astro-ph.CO] .\n- [76] S. May and V. Springel, 'Structure formation in large-volume cosmological simulations of fuzzy dark matter: impact of the non-linear dynamics,' Mon. Not. Roy. Astron. Soc. 506 no. 2, (2021) 2603-2618, arXiv:2101.01828 [astro-ph.CO] .\n- [77] L. M. Widrow and N. Kaiser, 'Using the Schrodinger equation to simulate collisionless matter,' Astrophys. J. Lett. 416 (1993) L71-L74.  \n[78] X. Li, L. Hui, and T. D. Yavetz, 'Oscillations and Random Walk of the Soliton Core in a Fuzzy Dark Matter Halo,' Phys. Rev. D 103 no. 2, (2021) 023508, arXiv:2011.11416 [astro-ph.CO] .  \n[79] A. Lague, B. Schwabe, R. Hloˇzek, D. J. E. Marsh, and K. K. Rogers, 'Cosmological simulations of mixed ultralight dark matter,' Phys. Rev. D 109 no. 4, (2024) 043507, arXiv:2310.20000 [astro-ph.CO] .  \n[80] M. Gorghetto, E. Hardy, J. March-Russell, N. Song, and S. M. West, 'Dark photon stars: formation and role as dark matter substructure,' JCAP 08 no. 08, (2022) 018, arXiv:2203.10100 [hep-ph] .  \n- [81] D. Budker, J. Eby, M. Gorghetto, M. Jiang, and G. Perez, 'A generic formation mechanism of ultralight dark matter solar halos,' JCAP 12 (2023) 021, arXiv:2306.12477 [hep-ph] .\n- [82] M. Gorghetto, E. Hardy, and G. Villadoro, 'More axion stars from strings,' JHEP 08 (2024) 126, arXiv:2405.19389 [hep-ph] .  \n[83] M. Frigo, 'A fast fourier transform compiler,' SIGPLAN Not. 34 no. 5, (May, 1999) 169-180.  \nhttps://doi.org/10.1145/301631.301661 .  \n[84] D. J. Bacon, D. M. Goldberg, B. T. P. Rowe, and A. N. Taylor, 'Weak gravitational flexion,' Mon. Not. Roy. Astron. Soc. 365 (2006) 414-428, arXiv:astro-ph/0504478 .", 'Appendix A: Schrodinger-Poisson solver': 'We use a pseudo-spectral 3D Schrodinger-Poisson solver, which closely resembles the one described in [24]. In the following, we review the basics of its implementations, together with some details of the particular current setup.  \nThe Schrodinger-Poisson Equations for n fields ψ i with mass m i read (in our conventions, ψ has the dimensions of a mass squared, in particular the density is ρ = | ψ | 2 )  \ni ∂ψ i ∂t = -1 2 m i ∇ 2 ψ i + m i (Φ + Φ ext ) ψ i , (A1)  \n∇ 2 Φ = 4 πG ∑ i ( | ψ i | 2 -⟨| ˜ ψ i | 2 |⟩ ) , (A2)  \nwhere Φ ext is a possible source of external gravitational potential, not sourced by ψ i , for which back-reaction is neglected.  \nTo convert to adimensional variables, define  \n˜ ψ i = 1 β √ 4 πG m 1 ψ , ˜ x = √ βm 1 x , ˜ t = βm 1 t , ˜ Φ = Φ /β , (A3)  \nwith β an adimensional rescaling parameter. This gives  \ni ∂ ˜ ψ i ∂ ˜ t = -1 2 r i ∇ 2 ˜ x ˜ ψ i + r i ( ˜ Φ+ ˜ Φ ext ) ˜ ψ i , (A4)  \n∇ 2 ˜ x ˜ Φ = ∑ i | ˜ ψ i | 2 -⟨ ∑ i | ˜ ψ i | 2 |⟩ , (A5)  \nwhere we defined  \nr i := m i m 1 . (A6)  \nThe energy associated to each field component reads  \n˜ E i = ∫ d 3 ˜ x ( | ∇ ˜ ψ i | 2 2 r 2 i + 1 2 | ˜ ψ i | 2 ˜ Φ+ | ˜ ψ i | 2 ˜ Φ ext ) . (A7)  \nThe field is evolved in time via the unitary operator  \n˜ ψ i ( ˜ t +d ˜ t ) = ∏ α e -i r i d α d ˜ t ( ˜ Φ α + ˜ Φ ext ) e -i c α d ˜ t ( -i ∇ ) 2 2 r i ˜ ψ i ( ˜ t ) .  \n(A8)  \nThe equation is meant to be read from right to left, i.e., one first applies the kinetic operator  \ne -i c α d ˜ t ( -i ∇ ) 2 2 r i ˜ ψ i ( ˜ t ) =: ˜ ψ ( α ) i ( ˜ t +d ˜ t ) , (A9)  \nand then the potential operator, where  \n∇ 2 ˜ Φ α = ∑ i | ˜ ψ ( α ) i | 2 . (A10)  \nThe constants c α , d α can be found in [24]. Eqs. (A9),(A10) are solved via Fast Fourier Transforms (FFT), using the FFTW library [83].  \nWe use adaptive time-steps, ensuring conservation of energy ∆ E/E ≲ 10 -5 between time-steps, and overall conservation of energy | E tot final -E tot initial | / | E tot final + E tot initial | ≲ 10 -3 .  \nThe code used for the results presented in this work is available upon request.', 'Appendix B: Shear and flexion terms': 'MSD-breaking effects for realistic ULDM cores, which have a finite radius and dynamical deviations from spherical symmetry due to time-varying interference patterns, generate shear and flexion terms that could, in principle, be sought after in lensing analyses unless they are smaller than or comparable to modeling uncertainties. Consider the lensing potential Ψ( ⃗ θ ) , defined from the convergence as  \n∇ 2 Ψ = 2 κ . (B1)  \nFor θ ≪ θ c , where θ c ≈ λ dB / (2 D L ) is an approximation for the core radius in angular variables, we can expand in ⃗y := ⃗ θ/θ c ,  \nΨ( ⃗y ) ≈ Ψ(0) + ∂ Ψ(0) ∂y i y i + 1 2 ∂ 2 Ψ(0) ∂y i ∂y j y i y j + 1 6 ∂ 3 Ψ(0) ∂y i ∂y j ∂y k y i y j y k + . . . . (B2)  \nFor systems at D L ≈ 1 Gpc and the core radii of the ULDM we are considering, y ≲ 0 . 2 in the region where images form.  \nThe first term in the expansion is an unobservable constant, whereas the second term is degenerate with the unobservable source position. The lowest order term which  \nhas an impact on inference analysis is the third, which amounts to convergence and shear  \nγ 1 ( ⃗ θ ) = 1 2 θ 2 c ( ∂ Ψ( ⃗y ) ∂y 1 ∂y 1 -∂ Ψ( ⃗y ) ∂y 2 ∂y 2 ) , γ 2 ( ⃗ θ ) = 1 θ 2 c ∂ Ψ( ⃗y ) ∂y 1 ∂y 2 . (B3)  \nSimilarly, the third derivatives of Ψ correspond to flexion terms F, G (for definitions, see e.g. [71, 84]). The expansion in Eq. (B2) converges fast enough if (we remind that Ψ has the dimensions of a squared angle)  \nθ c ∂ 3 Ψ(0) ∂θ i ∂θ j ∂θ k ≲ 1 = ⇒ F, G ≲ 1 θ c . (B4)  \nOur simulation outputs convergence on an equally-spaced grid, so it is natural to compute Ψ and its derivatives via FFT. From Eq. (B1) (hatted quantities are in Fourier space),  \nˆ Ψ ⃗ k = -2 k 2 ˆ κ ⃗ k . (B5)  \nIn a grid of length L with N points, with ⃗ j = ( j 1 , j 2 ) , j 1 , j 2 = 0 , . . . , N -1 labelling the grid point, we can write  \nΨ( ⃗ θ ) = -2 (2 π ) 2 ∫ d 2 k e i ⃗ k · ⃗ θ k 2 ˆ κ ⃗ k ≈ -2 (2 π ) 2 N/ 2 -1 ∑ l 1 ,l 2 = -N/ 2 e 2 π i ⃗ l · ⃗ j/N ˆ κ ⃗ k l 2 , (B6)  \nwhere we discretized the integral via θ i = j i L/N , k i = 2 πl i /L , l i = 0 , . . . , N -1 . With the discrete Fourier transform definitions 4  \nˆ κ ⃗ l = N/ 2 -1 ∑ j 1 ,j 2 = -N/ 2 κ ⃗ j e -2 π i ⃗ j · ⃗ l/N =: F N ( κ ) ⃗ l , (B7)  \nκ ⃗ j = 1 N 2 N/ 2 -1 ∑ l 1 ,l 2 = -N/ 2 ˆ κ ⃗ l e 2 π i ⃗ j · ⃗ l/N := F -1 N (ˆ κ ) ⃗ j , (B8)  \nfrom which  \nˆ κ ⃗ k = ( L N ) 2 F N ( κ ) ⃗ k , (B9)  \nwe can write  \nΨ = -2 ( L 2 π ) 2 F -1 N ( F N ( κ ) ⃗ l l 2 ) . (B10) Analogously, we can write, for Ψ derivatives,  \n∂ i Ψ = -i L π F -1 N ( F N ( κ ) ⃗ l l i l 2 ) , (B11)  \n∂ i ∂ j Ψ = 2 F -1 N ( F N ( κ ) ⃗ l l i l j l 2 ) , (B12)  \n∂ i ∂ j ∂ m Ψ = i 4 π L F -1 N ( F N ( κ ) ⃗ l l i l m l j l 2 ) . (B13)  \nIn Fig 7, we show γ 1 and F 1 as proxy for shear and flexion magnitude, computed using Eqs. (B12)-(B13). These results refer to the system showed in the bottom-left panel of Fig. 3, which has an high ∆ κ c ≳ 0 . 09 . With a large ∆ κ c , the system under consideration can yield possible detectable imaging errors. However, the magnitude of the shear near the center of the core is comparable to (and degenerate with) typical cosmological external shear, which is poorly constrained theoretically. Similarly, the flexion magnitude is not too far off from the flexion caused by a host group of the lens galaxy [71] (notice that criterium Eq. (B4) is satisfied, justifying the relevance of the expansion in Eq. (B2)). This calls for the possibility that MSD-breaking effects from cores with ∆ κ c ≳ 0 . 03 could still be reabsorbed by other modeling parameters. Once the features of these ULDM cores are better constrained theoretically (e.g. with cosmological simulations), one can use mock systems to explore the extent of such degeneracies.  \n<!-- image -->  \nFIG. 7: Plots of γ 1 and F 1 on the left and right respectively, referring to the convergence field in the bottom-left panel of Fig. 3.  \n<!-- image -->'}