Wake forces in a background of quadratically coupled mediators

Two particles can exert forces on each other when embedded in a sea of weakly-coupled particles. These"wake forces'' occur whenever the source and target particles have quadratic interactions with the mediating particles; they are proportional to the ambient energy density, and typically have a range of order the characteristic de Broglie wavelength of the background. The effect can be understood as source particles causing a disturbance in the background waves -- a wake -- which subsequently interacts with the target particles. Wake forces can be mediated by bosons or fermions, can have spin dependence, may be attractive or repulsive, and have a generally anisotropic spatial profile and range that depends on the phase-space distribution of the ambient particles. In this work, I investigate the application of wake forces to dark matter searches, recast existing limits on short-range forces into leading constraints on dark matter with quadratic couplings, and sketch out potential experimental modifications to optimize sensitivity. Wake forces occur in the Standard Model: the presence of the cosmic neutrino background induces a millimeter-range force about 22 orders of magnitude weaker than gravity. Wake forces may also be relevant in condensed-matter and atomic physics.


I. INTRODUCTION
The birth of the modern scientific method can arguably be traced to Galileo's gravitational experiments with inclined planes.His determination that the free-fall acceleration of objects is independent of their mass and composition, now known as the (weak) equivalence principle, is a cornerstone of the theory of gravitational interactions.Since then, swaths of experiments have been performed to test the equivalence principle, and the inversesquare-radius scaling of the two long-range forces of nature, gravity and electromagnetism.
These efforts were further boosted by the realization that motivated theories beyond the Standard Model (SM) of particle physics-the QCD axion [1][2][3][4][5][6][7] and large extra dimensions [8,9] in particular-can exhibit small deviations from these predictions or even qualitatively new forces [10].These forces are typically mediated by low-mass particles-dilatons, moduli, axions, dark photons, or gravitons-which couple linearly but weakly to regular matter.Single-particle-exchange forces generally have a range equal to the Compton wavelength of the mediator; any such particle with gravitational-strength, spin-independent couplings to matter must be heavier than 5.11 meV ≈ (38.6 µm) −1 from tests of the gravitational inverse-square law [11].
These light, weakly-coupled fields can also make up the dark matter (DM), as they are naturally long-lived and have generic production mechanisms in the early universe [12][13][14].Linearly-coupled DM fields generally give rise to temporally oscillating phenomena at a frequency equal to the mass.Light scalar DM causes oscillations of fundamental "constants" [15], which in turn lead to time-varying energies [15,16], length scales [17], and forces [15,18], the latter through gradients in the field.Pseudoscalar DM can excite electromagnetic fields [19], spin resonances [20], and acoustic modes [21].At higher DM masses, current and near-future single-quantum detection techniques are sufficiently sensitive to search for DM absorption [22,23] and conversion [24,25].
Quadratic interactions between DM and SM operators O SM are also possible, and are the leading interactions if the linear coupling is forbidden by a symmetry of the DM field.Examples include a Z 2 parity symmetry ϕ → −ϕ for a real scalar field, forbidding ϕO SM but allowing ϕ 2 O SM .Likewise, a U(1) symmetry Φ → e iα Φ of a complex scalar field implies a leading interaction of |Φ| 2 O SM in the effective field theory (EFT).Fermions beyond the SM (other than sterile neutrinos) necessarily have quadratic couplings to SM operators because of fermion number symmetry.Notably, the QCD axion a has an irreducible quadratic coupling to nucleons N = (p, n) of the form L ⊃ ( σ/2)(a/f a ) 2 N N , with σ ≈ 15 MeV and f a the axion decay constant [26].Finally, SM neutrinos couple quadratically to matter at low energies through the four-Fermi interaction: where the first and second lines contain the relevant neutral-and charged-current interactions, respectively, with g V ψ and g A ψ the vector and axial-vector couplings of the SM fermion ψ to the Z boson, and U the unitary matrix that diagonalizes the neutrino mass matrix.(In this convention, g V p = 1/2 − 2s 2 w , g A p = 1/2, g V n = −1/2, arXiv:2401.08745v2[hep-ph] 9 Sep 2024 , where s w is the sine of the weak mixing angle.) Many of the oscillatory phenomena that can occur for linearly-coupled DM carry over to quadratic interactions through simple rescalings of the coupling and the oscillation frequency (now twice the mass), see e.g.Refs.[27,28].At higher masses, the leading effect becomes scattering in low-threshold targets rather than (double) absorption [29].When the field has no ambient density, there are still static forces from a virtual exchange of two particles, but they decrease with distance as 1/r 4 for two-scalar exchange rather than 1/r 2 for single-scalar exchange, and is exponentially suppressed beyond half the Compton wavelength of the scalars [30].Famously, the feeble two-neutrino exchange forces scale as G 2 F /r 6 [19,31,32] at short distances.In this paper, I calculate the density-dependent forces between two SM particles or macroscopic bodies.If two SM particles are embedded in a sea of waves interacting with them via quadratic couplings, there exists a static force that scales linearly with the energy density, as the square of the quadratic couplings, and with a range of order the typical spatial (de Broglie) wavelength of the ambient waves, unless they have a high degree of coherence.The spatial profile of the resulting potentials is not generally spherically symmetric, and is determined by the couplings and the phase-space distribution of the mediating particles.I call this phenomenon a wake force, as it is analogous to the wake of a stationary boat (source particle) in a mild ocean swell (the ambient waves), which is felt by the nearby surrounding boats (the target particles).Curiously, wake forces generically violate Newton's third law-the vector sum of the wake force from source to target does not cancel that from target to source, since there is some momentum transferred (which depends on the source-target separation vector) to the ambient medium.
In Sec.II, I present the general formalism for calculating (scalar, non-derivative) wake forces both classically (Sec.II A) and quantum mechanically (Sec.II B), with analogous derivations for other scalar interactions and fermions relegated to Apps.A & B. I calculate the range and spatial profile of wake forces in Sec.II C. In Sec.III, I discuss applications of wake forces to DM searches (Sec.III A) and neutrino detection (Sec.III B).Sec.IV contains parametric comparisons of wake forces and potentials to other effects that necessarily occur for quadratically-coupled fields: double-exchange forces (Sec.IV A), elastic scattering (Sec.IV B), in-medium potentials and forces (Sec.IV C), and screening (Sec.IV D).A validation against nonperturbative numerical simulations is provided in App. C.
Certain aspects of this work have appeared in prior literature in different guises, and were developed independently of the material presented here.I comment on similarities and differences with Refs.[27,30,[33][34][35][36][37][38][39][40] in Sec.V, where I also discuss future directions and open questions.

II. THEORY
In this section, I develop the general framework for calculating the wake force between two particles, which can then be pair-wise integrated to obtain the force between two macroscopic bodies.The basic effect can be understood in simple classical terms, at the level of equations of motion and solutions thereof, presented in Sec.II A. Alternatively, Sec.II B contains the exactly equivalent quantized treatment using tree-level scattering amplitudes, which may be illuminating vis-à-vis standard derivations from single-particle virtual exchanges, and is more powerful for spin-dependent interactions.
The most minimal example is that of a real scalar field ϕ with a quadratic coupling to the number density current of a SM fermion ψ = {p, n, e, . . .}: The normalization of the interaction term proportional to both the ϕ mass m and the coupling constant G (in analogy to Fermi's constant) of inverse-mass-squared dimension is chosen for later notational convenience and to make contact with quadratically-coupled fermions.Ultraviolet (UV) completions of the EFT in Eq. 3 are discussed in Sec.III A. The calculations in the rest of this section are repeated for other quadratic scalar interactions in App.A (e.g.|Φ| 2 ψψ for a complex scalar Φ and ϕ 2 ψiγ 5 ψ) and for fermions in App.B. They follow the same steps as shown here.

A. Classical Description
Consider a static "source" ψ particle at the origin x = 0, and a "target" ψ particle at position x with r ≡ |x|.The equation of motion for ϕ in the presence of this single source particle, with a number density of ψψ = δ (3) (x), is: Assume that the ϕ medium is composed of a background wave ϕ 0 (t, x) with amplitude φ 0 and wavenumber k 0 , and a scattered wave δϕ(t, x): The angular frequency is ω = m 2 + k 2 0 in vacuum; ψdensity-dependent corrections are discussed in Sec.IV D. Sec.II B will generalize the simple background plane wave of Eq. 6 to the more phenomenologically relevant case of a random superposition of different k modes.
Eq. 4 can be solved perturbatively in G. Using the Green's function corresponding to the retarded propagator, one finds as the leading-order solution in G for the perturbed spherical wave δϕ.
A "target" ψ particle at x experiences the nonrelativistic potential: The first term is the in-medium potential from the background wave only, and is independent of the source particle's position.The second term is the leading-order correction from the scattered wave.Its time-averaged expectation value is the wake potential: The wake force is the gradient of this potential, and is universally attractive for relative separations smaller than the de Broglie wavelength r ≪ 1/|k 0 |.The parametric form of this wake potential, is generic: it persists for other types of (spinindependent) quadratic interactions, including with complex scalars and fermions.The form factor is defined to go to unity in the small-radius and nonrelativistic limit.For phase-space distributions beyond that of a single k mode in Eq. 9, the form factor falls off significantlypolynomially or exponentially-at distances larger than the typical de Broglie wavelength.The form factor is also generally aspherical (at large distances), unless the phase-space distribution of ϕ particles is spherically symmetric itself.Sec.II C will cover the form factor in greater detail.
Fig. 1 depicts the wake potential (bottom panel) for a monochromatic unidirectional wave moving in the zdirection k = (0, 0, 2π), as the multiplicative cross-term FIG. 1. Time-averaged wake potential (bottom panel) of Eq. 9 for a target test particle at position x produced by a monochromatic unidirectional background wave ϕ0(t = 0, x) of Eq. 6 moving in the +z direction (top left panel) and the scattered spherical wave δϕ(t = 0, x) of Eq. 7 (top right panel) from a perturbing source particle at the origin.For illustrative purposes, units are such that G = m = φ0 = 1 and k = (0, 0, 2π).
(the second term on the RHS of Eq. 8) of the background wave (top left panel) and the perturbed spherical wave (top right panel).For the monochromatic unidirectional background wave of Eq. 6, the wake force effectively has an infinite range in the forward direction x ≃ k, and is spatially oscillatory with wavenumber magnitude O(|k|) at large angles (2|k| in the backward direction).The generalization to a random superposition of background waves of varying wavenumber magnitudes and directions is crucial for phenomenological purposes, and is more easily derived in the quantized treatment below.

B. Quantum Description
The classical calculation from the previous section can be repeated from the perspective of tree-level scattering amplitudes in quantum field theory, where the generalization to an ensemble of many k modes and other types of interactions (Apps.A & B) is more transparent.The real scalar field of Eq. 3 is quantized in the usual way: with The environment, i.e. the ensemble of background waves, is taken to be a mixed state of k modes, with a momentum distribution f (k) normalized as d ¯3k f (k) = 1.In this background, the expectation value of the product of creation and annihilation operators is: with n the number density of ϕ particles.Using the freeparticle Hamiltonian density H = [ φ2 +(∇ϕ) 2 +m 2 ϕ 2 ]/2, the expectation value of the energy density is The Feynman diagrams contributing to the wake force are depicted in Fig. 2. Taking the nonrelativistic limit of distinguishable SM fermions ψ 1 and ψ 2 with masses m ψ , the matrix element M for this scattering amplitude is in the background of Eq. 12.The four-momentum exchanged between ψ 1 and ψ 2 is denoted by q with direction as indicated in Fig. 2. Use of the standard relation V (q) = −M/(2m ψ ) 2 between the matrix element and the 3D Fourier transform of the potential V (q) = d 3 x e −iq•x V (x) yields the result: where, once again, the iε prescription corresponds to the use of the retarded propagator to retain causality.
Fourier transforming back to position space gives the wake potential: In the second line, the integration variable is shifted as q + k → q.The result in the third line is equivalent to the classical result of Eq. 9 for f (k) = δ(3) (k − k 0 ) and with the identification n/E k0 = φ 2 0 /2.The long range of the wake force (of order the de Broglie wavelength) is explained by the nearly on-shell internal line: the square of the four-momentum k carried by the on-shell external state cancels the pole in the internal propagator in Fig. 2 and Eq. 13.In contrast, the short range of Yukawa-type forces (of order the Compton wavelength) from single-particle exchange and of double-particle-exchange forces can be attributed to their off-shell internal lines.
The parametric form promised in Eq. 10, is recovered in the nonrelativistic limit, where lim x→0 F(x) ≃ 1 and ρ ≃ mn.This is the relevant limit for the DM searches in Sec.III A, and will be used to calculate the range of the wake force in the next section.

C. Range
The range and profile of the wake force is governed by the momentum distribution f (k) of the background particles, which determines the form factor in Eq. 16: Isotropic distributions-For an isotropic Maxwell-Boltzmann (MB) distribution, the form factor takes on a simple analytic form (for E k /m ≃ 1): falling off as a (spherical) Gaussian beyond the characteristic de Broglie wavelength σ −1 k .It is depicted as the thick black curve in the top panel of Fig. 5.For an isotropic Fermi-Dirac (FD) distribution with temperature T : where ψ (1) is the first derivative of the digamma function.
In this case, the form factor is quartically suppressed for r ≫ 1/T , shown as the thick black curve in the bottom panel of Fig. 5.
For a general isotropic momentum distribution f (k) = f (|k|), the form factor integral reduces to: proving at least 1/r suppression at large r, and even more suppression for smooth phase-space distributions due to the oscillatory nature of the integrand.For a (physically less motivated) top-hat momentum distribution, , the form factor decreases in an oscillatory fashion with an inverse quadratic envelope, according to the parametric form F ∝ cos(2k 0 r)/(k 0 r) 2 for r ≫ 1/k 0 .The top-hat distribution shows that despite "hard edges" in f (k), there is still a significant suppression of the form factor beyond the typical wavelength k −1 0 .(These edges may come from incomplete phasespace equilibration at velocities below the solar system's escape velocity, or from the screening effects discussed in Sec.IV D.) Anisotropic distributions-The momentum distributions for many background fields of phenomenological interest, such as the DM and CνB particles, are strongly anisotropic in the laboratory frame due to our motion relative to that of the MW DM halo or the cosmic rest frame.
The DM halo in Earth's frame can be modeled to first approximation by a truncated, boosted Maxwell-Boltzmann (BMB) momentum distribution: with coordinates chosen such that the DM wind points towards the +z-direction (circular velocity v circ = −v circ ẑ).The distribution is truncated above the escape velocity v esc , which affects the normalization via N esc = erf(z)−2ze −z 2 / √ π, where z ≡ mv esc / √ 2σ k .The corresponding form factor F BMB cannot be calculated analytically.However, its numerical evaluation is shown in the top panels of Fig. 3 and Fig. 5 for fiducial parameters for the velocity distribution at the Sun's location in the Milky Way (MW): σ v ≡ σ k /m ≈ 165 km/s ≈ √ 2v circ and v esc ≈ 550.7 km/s [41].At distances larger than those shown in Fig. 3, the form factor falls off quadratically in all directions.The distance dependence in the forward (backward) direction is indicated by the red (gold) curve in Fig. 5.The exponential suppression for the isotropic MB distribution is thus lifted, but it still constitutes a dramatic suppression relative to the forward limit from Eq. 9-compare the bottom panels of Figs. 1 & 3.
In a standard cosmology, the CνB is expected to be at rest relative to the cosmic microwave background (CMB) at a temperature T ν ≃ (4/11) 1/3 T γ ≈ 1.95 K, where T γ ≈ 2.725 K is the CMB temperature [42].The CMB dipole anisotropy thus implies that the CνB is moving with re- spect to the Earth at a velocity of v CMB ≈ 370 km/s in the direction of the CMB dipole [43].The corresponding momentum distribution is a boosted Fermi-Dirac (BFD) distribution: Fig. 4 displays the form factor F BFD for the CνB in Earth's frame for m = m 3 = 60 meV, a fiducial value for the largest of the three neutrino masses, in agreement with current constraints on the sum of neutrino masses [44] and not far from the lowest possible value allowed by neutrino oscillation experiments [45,Ch. 14].The strength of the anisotropy is quantified by mv CMB /T ν , which equals about 0.44 for the assumed mass m 3 , and is significantly smaller for the lighter mass eigenstates.The form factor is again suppressed relative to the forward limit from Eq. 9. Like for the BMB distribution, the boost mitigates the distance suppression of the form factor to F BFD ∝ 1/r 2 at large radii, which is less severe than the 1/r 4 scaling of the unboosted F FDcfr. the bottom panel of Fig. 5.
The reason for the 1/r 3 scaling of the wake potential is analogous to the 1/r 3 effective potential of a charged particle in a Vlasov plasma with either a relative velocity or anisotropy [46].In that case too, the boost and/or anisotropy lift the exponential Debye screening of the potential to a polynomial (cubic) one.
The anisotropic momentum distributions f BMB and f BFD are but rudimentary approximations to the actual phase-space distributions of DM and the CνB, respectively.Effects that would certainly alter the DM form factor include: MW substructure such as streams and subhalos, the annual modulation from Earth's orbit around the Sun, environmental screening (Sec.IV D), and corrections to the phase-space distribution at low velocities (|k|/m ≲ 10 −4 ) due to the gravitational potential of the Sun.Similarly, the momentum distribution of the heaviest neutrino mass eigenstate likely deviates from a BFD distribution due to its gravitational clustering within the Virgo/Laniakea Supercluster, and, to a lesser extent, the MW halo and solar system.A detailed study of these corrections is left to future work.
Multiple fields-The form factor can parametrically change for multiple wake fields, possibly extending the range of the wake force through a phenomenon akin to neutrino oscillations.Suppose the Lagrangian of Eq. 3 is modified to an interaction with two real scalar fields ϕ 1 and ϕ 2 , with masses m 1 and m 2 : For concreteness, assume there is only an ambient ϕ 1 field with number density n and momentum distribution f (k), and that ϕ 2 is in its vacuum state.If m 1 = m 2 = m, one would get exactly the same results as for a single real scalar field of mass m, from the same diagrams as in Fig. 2, except with ϕ 2 on the internal line, and ϕ 1 on the external legs.
Qualitatively new dynamics occur for a mass splitting: If |∆| is greater than the typical momentum |k 0 | in f (k), then the wake force has a shorter range, of order 1/|∆| instead of 1/|k 0 |, because the propagator is further off shell.More interesting is the opposite limit: 0 < |∆ 2 | ≪ k 2 0 .First, the appropriate generalization of the form factor is: At distances r ≲ 1/|k 0 |, the form factor is almost the same as if ∆ 2 were to vanish exactly.
To illustrate the oscillation phenomenon most simply, consider an isotropic momentum distribution f (k) = f (|k|): with the approximation in the second line valid for |k 0 |r ≫ 1 and 0 < |∆|/|k 0 | ≪ 1, with |k 0 | again the typical momentum in f (|k|).A new length scale appears in the form factor: which is the inverse of the spatial "beat" wavenumber of the two oscillatory factors in the first line of Eq. 28, and is reminiscent of neutrino oscillations.As in that case, the oscillation length scale arises because of the mismatch of interaction and mass eigenstates, and is parametrically equal to the characteristic wavenumber divided by the mass-squared splitting.Thus, instead of being severely suppressed for r ≫ 1/|k 0 |, the form factor generally takes on a size largely independent of the precise shape of the momentum distribution.In a sense, the form factor is just linearly suppressed at r ∼ λ osc , since F(λ osc , ∆ 2 ) ∼ sgn(∆ 2 )/(|k 0 |λ osc ), as visually indicated by the gray dashed lines in Fig. 5.Note that the sign of the form factor-and thus the attractive or repulsive nature of the wake force-is commensurate with the sign of ∆ 2 over this distance range.The behavior for various ratios of |∆|/σ k and |∆|/T for ∆ 2 > 0 is plotted in Fig. 5; other than a sign change at large distances, the negative-∆ 2 case is similar.

III. APPLICATIONS A. Dark Matter Searches
In this section, I explore the application of wake forces to searches for DM.As a strawman example for a first case study, I focus on the quadratic scalar coupling to nucleons N = (p, n): though the generalization to the analogous electron coupling −(G s,e /2)mϕ 2 ēe is straightforward.(In particular, given that the number density of nucleons is roughly twice that of electrons in standard matter, one can rescale the resulting wake force from the nucleon coupling to the total coupling G s,tot ∼ = G s,N + G s,e /2 to first approximation.) The low-energy EFT may also include a parityviolating interaction with the pseudoscalar current −(G p,N /2)mϕ 2 N iγ 5 N , analogous (derivative) couplings to the pseudovector current, or similar parity-violating interactions with electron currents.These will lead to spin-dependent dipole-dipole wake forces proportional to G 2 p , or in combination with the quadratic scalarcurrent couplings of Eq.31, to spin-dependent monopoledipole wake forces proportional to G s G p , as calculated in App. A. Their employment in the search for DM is fruit for further studies.
As shown in App.B, quadratic fermionic couplings of the form −(G s,N /2) χχ N N lead to precisely the same wake force as Eq. 31 in the nonrelativistic limit, so in principle wake forces can be used to search for fermionic DM particles χ as well.However, the relatively short range of the wake force (while at least 10 3 times longer than that of the double-exchange force of Sec.IV A in vacuum) precludes its practical use in searches for realistic fermionic DM candidates, due to the lower bound on the mass: m χ ≳ keV [47].This confines the (unsuppressed) wake force range to the sub-micron regime, where current experiments are less sensitive.
Models-The simplest UV completion for the quadratic scalar coupling of Eq. 31 is a quadratic coupling to the Higgs doublet H: with the linear term ϕH † H forbidden by a Z 2 symmetry, which also guarantees the stability of ϕ as a DM candidate (by forbidding decays such as ϕ → γγ).At low energies, where the Higgs can be safely integrated out, effective quadratic couplings to nucleons and electrons are generated: , where m N , m e , m h are the nucleon, electron, and Higgs masses.
The Higgs portal coupling g H is constrained by the observation that the invisible branching ratio of the Higgs, to which the process h → ϕϕ contributes, is less than 0.18 [50] (falling pink line in Fig. 6).The g H coupling also leads to quantum corrections to the ϕ quartic L ⊃ −λ ϕ ϕ 4 /4! of order δλ ϕ ∼ g 2 H /16π 2 , and thus a minimum characteric single-particle self-interaction crosssection σ ϕ ≳ δλ 2 ϕ /(128πm 2 ) [51].Astrophysical observations of merging clusters and dwarf galaxies require that (1+f occ )σ ϕ /m ≲ 1 cm 2 /g [52], where f occ ≃ (2π) 3/2 n/σ 3 k is a Bose enhancement factor [53]. (This coherent enhancement of f occ has not been discussed in previous literature on self-interacting DM, where the DM constituents are usually assumed to be sufficiently heavy such that f occ ≪ 1.) Self-interactions thus indirectly bound g H from above as shown by the other pink solid line in Fig. 6, for which I take conservative values of ρ DM = 5 × 10 6 M ⊙ /kpc 3 and σ v = 3,000 km/s for the Bullet Cluster.A conservative estimated bound based on large-scale structure formation for a natural quartic of |δλ ϕ | ≲ 5 × 10 −8 (m/eV) 4 [15, Sec.VII] is shown as a dotted pink line.Ref. [49] reports Big Bang nucleosynthesis (BBN) constraints on quadratically-coupled dark matter; the corresponding limit on the nuclear coupling is shown by the dashed pink line, but would be dramatically less stringent for the Higgs UV completion, which also includes couplings to electrons and photons which induce a high effective mass for ϕ in the early universe.Similarly, ϕ receives a mass-squared quantum correction δm 2 ∼ g H Λ 2 /16π 2 with Λ the UV cutoff, which implies that this model has a severe hierarchy problem.As for the Higgs, one may resort to anthropic selection in a landscape of vacua to explain the smallness of m 2 .
Several of these observational and naturalness constraints can be alleviated (at the expense of introducing others) by introducing a combination of linear and cubic interactions.For example, the super-renormalizable Lagrangian leads to the same low-energy wake force coupling as Eq.32 with the identification g H ∼ = A H A ϕ /m 2 .Its super-renormalizability means that the model of Eq. 33 is UV safe: in particular, the invisible branching ratio of the Higgs is parametrically smaller than for the "hard" g H coupling, and ϕ does not suffer from a hierarchy problem.An analogous self-interaction cross-section bound applies to A 2 ϕ , and there are significant singleexchange (Yukawa) force and stellar cooling constraints on A H .One must also address irreducible cosmological production channels, decays of ϕ proportional to A 2 H , and the phenomenology of the electroweak phase transition within these classes of models, all left to future work.
Other pathways to generate quadratic interactions of the form −(G s,ψ m/2)ϕ 2 ψψ and −(G p,ψ m/2)ϕ 2 ψiγ 5 ψ are via a light mediator a with mass m a , which couples to ϕ with coupling g s ϕ and to ψ with couplings g s ψ and g p ψ : (An alternative and more minimal possibility is, as in Eq. 33, through the combination of a cubic selfinteraction ϕ 3 and linear couplings ϕ ψψ and ϕ ψiγ 5 ψ.) Integrating out the mediator a, valid for momentum exchange max{|q|, |k|} ≪ m a yields low-energy quadratic couplings of G s,ψ ≃ g s ϕ g s ψ /m 2 a and G p,ψ ≃ g s ϕ g p ψ /m 2 a .A full charting of the parameter space of couplings g s ϕ , g s ψ , g p ψ , and masses m and m a is beyond the scope of this work, but see e.g.Refs.[28,29] in this direction.
Arguably the best-motivated quadratic scalar coupling to nucleons is that of the QCD axion a, which irreducibly has (in the parametrization −(G s,N /2)m a a 2 N N of Eq. 31): with σ ≈ 15 MeV [26] and f a ≈ 10 12 GeV (5.7 µeV/m a ) as a function of the QCD axion mass m a [54].This coupling is depicted by the green line in Fig. 6 for f a ≳ 10 8 GeV, the approximate lower bound on the decay constant.FIG. 6. Parameter space of quadratic scalar coupling Gs,N to nucleons as a function of ϕ mass m, in the interaction L ⊃ −(Gs,N m/2)ϕ 2 N N .In the blue region, the wake force is within reach of existing monopole-monopole force experiments if ϕ makes up all the DM.The other blue lines indicate wake forces detectable by a torsion balance (TB) with 1 cm-radius tungsten spheres as source and target masses, angular frequency ω = 2π mHz, quality factor Q = 10 8 , and noise temperatures T = 10 mK (solid) and T = ω/2 (dashed, SQL).In the light brown region, the wake force is screened by matter at standard densities.The characteristic range 1/σ k is indicated by the top axis.Below the upper black dashed line, the DM wake force between nucleons is weaker than gravity at short distances; the lower black dashed line is the reference point Gs,N = GF for electroweak-strength couplings.The yellow dashed curve corresponds to a single scattering event per kg-yr exposure in superfluid helium [48].The pink solid lines are upper bounds on the effective quadratic coupling in the Higgs model of Eq. 32 from DM self-interactions (ϕϕ → ϕϕ) and from invisible Higgs decays (h → ϕϕ).The dotted pink line is the structure formation exclusion from Ref. [15], and the dashed pink line the BBN bound [49] on the nuclear coupling only but is less stringent in the quadratic Higgs UV completion.The green line delineates the quadratic nuclear coupling of the QCD axion for fa > 10 8 GeV.
Approximate recasted limits-Dating back to the experiments of Cavendish, there has been a tremendous effort to search for deviations from Newton's gravitational potential, usually parametrized by the Yukawa potential with strength α (relative to gravity) and range λ: as would be generated from a single virtual exchange of a particle with mass m = 1/λ and linear coupling √ 4πG N α.The form of the wake potential (Eq.16) is similar at short distances, but the exponential suppression is replaced by the anisotropic form factor F(x).The existing limits on α lim (λ) [11,[55][56][57][58][59][60][61][62][63][64][65][66] can be (approximately) recast into a value of the quadratic coupling G lim s,N above which the wake force would have been detectable in past experimental searches-had efforts been made to look for it.I take this recasting map to be: The m dependence of G lim s,N (m) is implicit in the form factor F BMB in the forward direction of the boosted Maxwell-Boltzmann distribution with σ v = σ k /m ≈ 165 km/second and v circ ≃ √ 2σ v , plotted as the red line in the top panel of Fig. 5.At large m (and fixed r ∼ λ), the form factor decouples as F ∝ 1/m 2 , so that G lim s,N (m) ∝ 1/m, with the experimental limit from Ref. [66] dominating this tentative "constraint".At small m, the form factor is constant, so that G lim s,N tends to a constant.In the latter regime, the wake potential scales as 1/r and is only distinguishable from gravity via its generic violation of the weak equivalence principle (i.e.dependence on chemical composition) and the apparent mismatch of the effective value of G N at short range and on astronomical scales.
The map of Eq. 37, shown as the blue filled region in Fig. 6, cannot be deemed a limit with a quantifiable confidence level, because the geometry and orientation of the experiment relative to the DM wind must be taken into account.The experiments with leading limits on α(λ) did not explicitly search for composition dependence (only deviations from the inverse-square-radius behavior of the force), although they can be (and have been) modified to do so [67].Nevertheless, the recasting via Eq.37 should be an accurate (if conservative) guide to the approximate sensitivity of existing experimental setups.
Sensitivity projections-The sensitivity of dedicated wake force experiments can be estimated faithfully.As a representative case, consider the wake force exerted on a "target" tungsten sphere of radius R = 1 cm and mass density ρ N = 19.28g/cm 3 positioned in the forward direction of the wake emanating from an identical "source" sphere, with a center-to-center separation of 2.2R.This arrangement should be close to the optimal geometry given a target of size R: making the source much larger than the target does not parametrically increase the signal in the regime of large mσ k R in view of the largeradius scaling of the wake force |F| ∝ F BMB /r 2 ∝ 1/r 4 .For this geometry, the wake force on the target points towards the source, and has a parametric size: with M = (4π/3)ρ N R 3 ≈ 81 g the mass of the spheres.Suppose the target sphere is confined to its position in a (generalized) harmonic oscillator potential well of angular frequency ω = 2π mHz and losses corresponding to a quality factor Q = 10 8 , as one may obtain with a future torsion-balance (TB) setup [39].In this rather generic case, the minimal variance on the force measurement after an integration time t int = 3 yr is given by the fluctuation-dissipation theorem [68]: In the above, E(ω, T ) is the noise temperature, which is typically (at least) the thermodynamic temperature T of the environment; in thermal equilibrium, E(ω, T ) = ω 1/2 + 1/(e ω/T − 1) [69, Ch. 17.2].With active cooling methods and a sufficiently decoupled oscillator, one could in principle achieve the standard quantum limit (SQL) of E(ω, T ) = ω/2 ≈ 2 × 10 −14 K, though presentday capabilities are far from these values.Eq. 39 represents the minimum contribution from thermal/quantum noise, and does not include readout/back-action and environmental noise sources.
This corresponding idealized coupling sensitivity, at unit signal-to-noise ratio |F|/σ F = 1, of this representative oscillator is depicted in Fig. 6 as blue lines, at the thermal noise limit E = T = 10 mK (solid) and at the SQL E = ω/2 (dashed).For this forecast, the size of the wake force between the spheres is computed numerically, as opposed to the parametric estimate of Eq. 38, and the target-source separation is assumed to remain aligned with v circ of Eq. 23.Any realized sensitivity will be degraded by an order-unity factor based on the desired discovery or exclusion threshold, and for a non-corotating experimental setup.Especially the thermal-noise-limited sensitivity is a benchmark that current experimental setups may strive towards; attaining it would have significant discovery potential for the DM wake force of e.g. the quadratic Higgs coupling in Eq. 32.Dedicated wake force experiments may mitigate systematics through various control knobs, such as its weak-equivalence-principleviolating nature, distance dependence, and most strikingly, its asphericity and angular variation with respect to the DM wind direction, which changes on a diurnal basis in the lab frame.It would be interesting to explore the optimal experimental geometry and orientation to maximize the wake force signal and minimize more prosaic backgrounds.For example, short-range anisotropic corrections to gravity can be sensed by modified torsionbalance setups [70].
Other DM experiments-Traditional DM experiments searching for isolated scattering events (which produce low-energy quasiparticles in the scattering target) lose sensitivity at low DM masses, where the recoil energy is suppressed.The effective coupling at which one expects a single scattering event for a kg-yr exposure of the low-threshold superfluid helium proposal from Refs.[48,71] is shown in Fig. 6 as the dashed yellow line.Wake force experiments instead become more sensitive at low masses due to the increased range, and may thus probe effective single-particle cross-sections far below 10 −40 cm 2 and even 10 −50 cm 2 (gray dashed lines).In this sense, DM wake force searches can bridge the gap between traditional DM scattering experiments and coherent interactions of bosonic DM with macroscopic targets (haloscopes).

B. Neutrino Detection
Cosmic neutrino background (CνB)-The presentday CνB is, in a standard cosmology, expected to have a temperature: which also sets the characteristic momentum and the corresponding coherence length T −1 ν of the neutrino background (and wake force).The number densities per mass eigenstate are n νi = nνi ≈ 56 cm −3 .Given the observed mass-squared splittings of the neutrino mass eigenstates, at least two (and possibly all three) are nonrelativistic.As shown in App.B and in Refs.[33,34,[36][37][38][39][40], wake forces can be mediated by fermions too.
In App.B 2, I calculate the full wake potential for a single Dirac-neutrino mass eigenstate interacting through neutral-current interactions with SM fermions ψ with vector and axial-vector couplings g V ψ and g A ψ , respectively.Those results are easily generalized to include three Dirac neutrinos and the charged-current interactions.However, the wake potential is largest for the highest-mass neutrino eigenstate, so in the "normal" mass hierarchy case with strong ordering m 1 ≲ m 2 ≪ m 3 , it is a fine approximation up to O(m 2 /m 3 ) corrections to take into account only the heaviest neutrino ν 3 .Neutrino oscillations from the charged-current interactions can similarly be neglected, as the mass-squared splittings between the neutrinos are higher than the square of the neutrino temperature: |∆ 2 | ≫ T 2 in the language of the discussion around Eqs. 26-29.With those assumptions, the monopole-monopole wake potential mediated by the CνB between two SM fermions ψ 1 and ψ 2 can be translated from Eq. B20 (for Dirac neutrinos) to: where δ ψe is the effective vector "charge" of the fermion ψ as seen by ν 3 , including both the neutral-current vector coupling g V ψ and the chargedcurrent vector coupling U 3e U † e3 (only to electrons).The form factor F V V (x) is given in Eq.B21 (where one should take m = m 3 ), which in the nonrelativistic limit matches Eq. 17 up to O(T 2 /m 2 ) corrections and thus that of the top panel in Fig. 4 for a BFD distribution.
The characteristic strength of the CνB wake force can be expressed in terms of its relative size α CνB relative to the gravitational force between two macroscopic bodies 1 and 2: at short distances r ≪ T −1 ν ≈ 1.17 mm, where Q 3,V i is the average effective vector charge per nucleon of the heaviest neutrino mass eigenstate in the two bodies i = 1, 2. The present best constraint on α(T −1 ν ) is about 10 −3 [63], far worse than Eq.43.An improvement by 19 orders of magnitude is not in reach of any conceivable force experiment; even the SQL-limited sensitivity of the dashed line in Fig. 6 would "only" constitute a 10 13 improvement in α(T −1 ν ).The monopole-dipole wake potential mediated by the CνB between two SM fermions ψ 1 and ψ 2 can similarly be translated from Eqs. B22 and B23.The precise expression of the form factor is not important for a back-ofthe-envelope estimate of its size, which is O(T /m).For a spherical tungsten source mass 1 of radius R ∼ T −1 (larger sources are not relevant given the short range), the integrated monopole-dipole wake potential for a polarized target with (nuclear or electronic) spin σ 2 is of the parametric form V (σ 2 ) ∼ Ω • σ 2 , with an angular precession frequency: In comparison, the CνB wind effect linear in G F from Ref. [72] is of the form The precession induced from the wake potential in Eq. 44 is larger than that of Eq. 46 if the neutrino asymmetry is as small as the baryon asymmetry, i.e. n ν − nν ∼ 10 −10 (n ν + n ν ).For neutrino asymmetries larger than about 10 −5 , the effect from Eq. 46 would exceed that of the CνB wake potential in any geometry.However, either effect appears to be far below the quantum projection noise limits for nuclear and electron spins [73].Solar and reactor neutrinos-There exist more intense sources of relativistic neutrinos, for which the heuristic form of the wake potential in Eq. 10 still holds, now with ρ ≃ |k 0 |n instead of ρ ≃ mn, as can be seen from taking the relativistic limits of Eqs.B21 and B23.The solar neutrino flux and thus number density is about n ν,⊙ ≈ 7 × 10 10 cm −2 s −1 ≈ 2 cm −3 , but with typical energies of E k0 ∼ MeV.Solar neutrinos lead to a fractional correction to gravity on the order of α ν,⊙ ∼ 2G 2 F E k0 n ⊙,ν /(4πG N m 2 N ) ∼ 10 −16 at distances of k −1 0 ∼ 10 −13 m, again far too small to be detected, even after taking into account the neutrino oscillation effects described in Sec.II C. Fluxes from reactor neutrinos, whose individual coherent scattering events have been detected already [74], can be a couple of orders of magnitude larger, but not enough to close the gap towards practical observability of neutrino wake forces with present-day technology.(The same conclusion holds for neutrinos from other sources, such as supernovae, radioactive samples, and spallation.)

IV. COMPARISONS
In this section, I compare the size of wake potential effects to other phenomena that necessarily occur in the same theories with quadratic interactions to matter.Other than in Sec.IV A, the focus will be on the interaction of a light field ϕ coupled to nucleons as in Eq.31, since some comparisons for neutrino wake forces have already been made in Sec.III B.

A. Double-Exchange Forces
The diagrams in Fig. 2 can be regarded as "cut" versions of the two-particle-exchange diagram in Fig. 7, with the ϕ propagator split into two external legs which have nonzero occupation number.The resulting wake potential is a purely classical wave effect, as evident from Sec. II A. The double-exchange potential from Fig. 7 is a quantum effect arising from vacuum fluctuations, as can be shown by loop-or ℏ-counting arguments.For the coupling in Eq. 3, the double-scalar-exchange potential between two ψ particles is [30,75,76]: where K 1 is a modified Bessel function of the second kind.Firstly, the range of the force from the exchange of two virtual scalars is roughly half the Compton wavelength of ϕ, whereas wake forces have a range of order the de Broglie wavelength of the ambient ϕ particles.Secondly, for a force experiment carried out at r = 1/2m, the ratio of the double-exchange and wake potentials is: The wake force is enhanced at large number densities, and, unlike the double-exchange force, is not suppressed by a loop factor.A range of 1 mm or above, where force experiments achieve their best sensitivity, corresponds to m ≲ 1 meV to avoid exponential suppression in the double-exchange force.If ϕ makes up all of the DM, then this implies n/m 3 ≳ 10 6 and thus a much larger wake force than double-exchange force.The occupation number decreases quickly for sub-millimeter Compton wavelengths, but so do the sensitivities of force experiments, so that the longer range of the wake force by a factor of O(σ −1 v ) is a major advantage.Thus, if ϕ is a nonnegligible fraction of DM, then its wake force is always the dominant effect in force experiments, and the double-exchange force is negligible.(I do not attempt to recast α lim (λ) in terms of the double-exchange force sensitivity to G s,N , because current experimental setups are not sensitive to couplings in the perturbative, unscreened regime below the brown region in Fig. 6.) Similar arguments apply to double-fermion-exchange potentials.For SM neutrinos at distances r ≪ (2m ν ) −1 , they are of the form V 2ν ∼ G 2 F /(4π 3 r 5 ) and contain both spin-independent and spin-dependent (dipoledipole) components [31,32,77].This phenomenon is far too small to be observable with current technology, except possibly at very short distances with future experimental and theoretical efforts in muonium spectroscopy [78].For fermionic DM, one must have n/m 3 ≲ σ 3 v by the Pauli exclusion principle, which is the primary reason for why fermionic DM wake forces are harder to detect.

B. Elastic scattering
For m ≪ m ψ and for momentum transfers where the nonrenormalizable operator of Eq. 3 is valid, ϕ particles can scatter elastically off SM ψ particles, with a crosssection: In the limit of m ≫ m ψ or where the momentum transfer in the interaction is above the effective cutoff of the interaction (as would be possible in e.g. the light-mediator UV completion of Eq. 34), there are generally further suppression factors, which only strengthen the point below.
Individual scattering events occur stochastically, and are in practice only measurable if suitable quasiparticles with sufficiently low energy threshold can be excited and detected.The ultimate reach for single-phonon excitation in superfluid helium [48,71], which incurs additional suppression factors beyond Eq.49, is shown as the dashed yellow line in Fig. 6.An alternative approach is to look for the collective, integrated recoil of many elastic scattering events on a macroscopic target over time; each scattering event will, on average, impart some momentum to the target in the direction of the DM wind [79].
In a bath of nonrelativistic ϕ particles with number density n and typical momentum magnitude |k 0 |, the scattering rate per ψ particle is Γ ∼ nσ el |k 0 |/m.If the average momentum transfer per scattering event is k 0 , the time-averaged force due to elastic scattering, on a single ψ particle, is: The wake force F = −∇V (x) induced by the potential in Eq. 16 with e.g. the MB form factor of Eq. 19 (for simplicity) is equal to: Comparison of Eqs.50 & 51 shows that the elastic scattering and wake forces on a single ψ particle are roughly equal in magnitude when the two ψ particles separated by a distance r 1/σ k and if |k 0 ||k 0 | ∼ σ 2 k .However, the wake force is larger than the elastic scattering force for r ≲ 1/σ k or if the average momentum transfer |k 0 | is much smaller than the typical momentum, which is the case for the CνB and likely other types of primordial dark radiation.
To maximize detection prospects, one generally wants to measure either force on as many ψ particles as possible.The elastic scattering force is coherent when the target size is smaller than the coherence length of ϕ, R target ≲ 1/σ k , leading to a scaling of (n ψ,target R 3 target ) 2 with the number of particles in that regime.The wake force enjoys a similar coherent enhancement as (n ψ,target R 3 target ) × (n ψ,source R 3 source ) for R target , R source ≲ 1/σ k .However, in certain cases of practical relevance (e.g.torsion balances, optically-levitated dielectrics, atom interferometry), maximum acceleration sensitivity may be achieved on a very small target.Having the option to use a larger source of ψ particles over which the wake force is coherent can be advantageous in such cases.Furthermore, the wake force is generally easier to control and manipulate, as it depends on the density, size, distance, and (generally) chemical composition of the source, allowing for a variety of experimental knobs to be turned.Both coherent elastic scattering forces and wake forces enjoy a cosmic reference direction, which should help mitigate systematics.However, like any pair-wise force, wake forces can be made to vary temporally-by employing rotary stages [11,58,60,[64][65][66], driven mechanical oscillators [56], or linear motions of source [55,57,[61][62][63] or target masses [59]-allowing them to be probed by AC experiments as opposed to the daunting experimental challenge of searching for an uncontrollable DC force.Ref. [79] computes elastic scattering and the associated DC forces in the regime where screening is important, and discusses detection prospects of these forces with satellite tests of the equivalence principle and torsion balance experiments.

C. In-Medium Potentials and Forces
In a sea of ϕ particles, individual ψ particles experience an in-medium potential and force at linear order in G, represented by the right three diagrams of Fig. 7, independent of any neighboring ψ particles.For the quadratic scalar field interaction of Eq. 3, this in-medium potential is simply: with the second equality for a simple plane wave of ϕ, for which the in-medium force is: Diagrammatically, the "constant" term in Eq. 52 arises from the same diagram as elastic scattering (second from left in Fig. 7, with both positive-and negative frequency factors), whereas the oscillatory terms at 2ω in Eq. 52 and the force of Eq. 53 arise from the rightmost two diagrams in Fig. 7.For a random superposition of waves, n itself is stochastically varying with coherence times and lengths of order m/σ 2 k and 1/σ k , respectively, generating additional forces from spatial gradients in n, which are of the same order of magnitude as the oscillatory force in Eq. 53.
In the perturbative regime for G, these in-medium potentials and forces are naively larger than their corresponding wake counterparts.In particular, V wake /V IM ∼ Gm/2πr ≲ 1 for a single source particle due to perturbativity; for many source particles over a de Broglie volume of O(k −3 0 ), V wake /V IM ∼ Gmn ψ /k 2 0 ≲ 1 for a cosmic field when screening effects are negligible (Sec.IV D).However, in the "high-mass" regime of Fig. 6, the fractional nucleon mass variation is ), well out of reach of atomicclock sensitivity or that of other metrological experiments.The force from Eq. 53 oscillates at a frequency of m/2π ∼ 242 THz (m/eV), too rapidly to be detected by torsion balances or other force sensors over the mass range in Fig. 6.Similarly, the stochastic acceleration from gradients in n is about 10 −12 m/s 2 (G s,N /GeV −2 ) independent of mass, but also has a prohibitively short coherence time of 2π/(mσ 2 v /2) ≈ 273 µs (eV/m), after which it averages to zero.
The considerations for in-medium potentials from fermions are similar; the effect from Ref. [72] and Eq. 45 is one such spin-dependent incarnation for the CνB, and the only known nonstochastic effect linear in G F for the CνB [80].

D. Screening
In this work, I have so far treated the quadratic ϕ interactions perturbatively, in which the separation between a cosmic background wave and a linear superposition of scattered waves from all perturbing sources is justified.However, at large G and/or small momentum, this perturbative expansion can break down [27,38,81].For an arbitrary but static source distribution with number density n ψ (x), the full ϕ field obeys the Klein-Gordon equation with a spatially varying effective mass-squared m 2 eff = m 2 + Gmn ψ (x), which one can model as a corresponding index of refraction [38,39].
To illustrate the breakdown of the perturbative treatment at large G and/or small momenta, consider the 1D problem in the z-direction with a step function in number density, which vanishes for z < 0 and is a constant n ψ for z > 0. Suppose there is an incoming wave of the form cos(ωt − k 0 z) for z < 0. A wave will be reflected from the step at z = 0: R c cos(ωt + k 0 z) + R s sin(ωt + k 0 z).wise, a wave is transmitted inside the medium (z > 0) with wavenumber: For real k ψ , write the transmitted solution as T + c cos(ωt− k ψ z) + T + s sin(ωt − k ψ z), and for imaginary k ψ , as e −|k ψ |z [T − c cos(ωt) + T − s sin(ωt)].This simple boundary problem can be solved for the reflection and transmission coefficients: for real ( + ) and imaginary ( − ) k ψ , respectively.The fraction of transmitted power is 4k 0 k π /(k 0 + k ψ ) 2 for Gmn ψ < k 2 0 and vanishes for Gmn ψ > k 2 0 .Either way, most or all of the incoming wave is reflected when |G|mn ψ > k 2 0 .(Disastrous effects due to tachyonic instabilities could theoretically occur for large negative G, which I will not consider here.) This "screening" effect has important implications for detection of the wake force from a cosmic field.For sufficiently large coupling, most of the "dark" field could be reflected from the environment, and not make it to the laboratory setup.In general, the perturbative expansion can be expected to hold at large distances if where n ψ is the typical density of the environment (Earth's crust or atmosphere, laboratory walls).In Fig. 6, the regime where Eq. 59 is not satisfied is shaded in light brown, for n ψ corresponding to the nucleon density of aluminum.Even when Eq. 59 is satisfied, some small fraction of incoming waves will be shielded by the environment [38], but when the inequality is strong, the perturbative expansion of this work should be valid to a good approximation.
A full calculation of screening effects would be highly specific to the geometry of the system under study, and is left to future work, though the approach of Ref. [27] should be valid for DM masses m ≪ 10 −11 eV, i.e. for vacuum de Broglie wavelengths larger than Earth's radius.In any case, given self-interaction and structure formation bounds on the model space discussed in Sec.III A, there is a significant theoretical bias (for technically natural self-couplings) towards low values of |G| that satisfy Eq. 59. Close to the boundary of Eq. 59, the calculation of the wake force to O(G 3 ) may be merited.Based on the simulations of App.C, I speculate that the next-to-leading order wake potential will generate a fractional correction to the leading O(G 2 ) wake force of O(Gmn ψ /k 2 0 ).Shielding and screening effects are already negligible at the leading edge of present-day sensitivity (lower boundary of blue region in Fig. 6), and will be even more so for future (wake) force experiments.

V. DISCUSSION
In this work, I have introduced a general formalism for the perturbative calculation of wake forces and potentials, which can be understood as classical wave phenomena that generically occur in a "sea" of particles quadratically coupled to nonrelativistic source and target particles.I have computed the analytic form of wake forces from a variety of quadratic interactions of both (pseudo)scalars (Sec.II A, Sec.II B, and App.A) and fermions (App.B), as well as their spatial profiles and range (Sec.II C).Before outlining potential future applications and directions, I will comment on the overlap and differences with previous literature (of which I am aware).
Firstly, Ref. [27] is a study of both linearly-and quadratically-coupled light scalar DM, wherein the inmedium potentials of Sec.IV C and some screening phenomena (Sec.IV D) were discussed for spherical sources (Earth in particular).The screening effects are nothing but the nonperturbative manifestation of the wake potential.The approach of Ref. [27] is valid in the regime where the de Broglie wavelength of the mediating particles is much larger than the size of the source, but finite-wavenumber effects and thus the short-range nature of the wake force were neglected, so their constraints should be quantitatively reconsidered for m ≳ 10 −11 eV.Separately, Ref. [35] found an attractive 1/r 2 force on galactic scales between test particles in a Bose-Einstein condensate (BEC) of quadratically-coupled scalar field DM with self-interactions.In the limit of vanishing selfinteractions, their results are equivalent to the perturbative wake force formalism of Sec.II A at zero background momentum.
Wake forces in a neutrino medium have also been posited in Refs.[33,34,[36][37][38][39], and in a bosonic medium in Refs.[30,39], mostly independent of this work.Refs.[30,36,37] regarded the medium-dependent forces as a version of the two-particle-exchange diagram, with the effective propagators modified by the (thermal) background.Refs.[36,37] also found an additional 1/(T r) 4 suppression of the G 2 F ρ CνB /r potential at large distances for an isotropic neutrino background.Ref. [30] primarily discussed the case of a Bose-Einstein condensate background.The background-modified propagator approach of Refs.[33,34,36,37] is equivalent to my approach, though the kinetic-theory formalism of Sec.II B is more efficient for wake force calculations for arbitrary, anisotropic phase-space distributions.(Eqs.B20, B21, and 21 match Eqs.3.24 and 3.25 of Ref. [37] with n ν = n ν = 3ζ(3)T 3 ν /4π 2 per generation.)Refs.[38,[82][83][84] compute the density modulation of the neutrino background near Earth's surface at the level of the nonrelativistic field equations with a spatially varying index of refraction, which should again be equivalent the formalism presented here in the respective applicable regimes.App.C constitutes partial numerical confirmation of both the perturbative results of the main text and of the conclusions of Refs.[83,84].While this work was being completed, Ref. [40] came to my attention; their wavepacket approach for neutrinos more closely resembles the kinetic-theory formalism of Sec.II B.
The primary novelty of this work is its more illuminating treatment of wake forces as perturbative classical wave phenomena in Sec.II A and within the framework quantized kinetic theory in Sec.II B. The methods presented here should agree with those of the abovementioned works in their regimes of validity.However, the formalism of this work allows for the efficient computation of the spatial profile of wake forces (Sec.II C) for phenomenologically realistic background distributions (Figs. 3, 4, and 5), including novel oscillation phenomena.Apps.A and B show that wake forces arise and can be computed for any quadratically-coupled field.App.C describes the implementation and results of nonperturbative numerical simulations of wake forces, which validate the perturbative methods in this work in 2 + 1 and 3 + 1 dimensions.
I have shown in Sec.III A that wake forces can form the basis for a new class of DM searches with precisionfrontier force experiments, in a mass range that is challenging to probe with other methods.Force experiments can thus bridge the gap between traditional DM scattering experiments on the one hand, and DM searches based on coherent interactions to first order in the coupling on the other hand.The optimization of sensors and setups for wake forces should be explored further, and the parameter space of relevant DM models charted out more extensively.Spin-dependent DM wake force searches are an obvious next step.A formal proof of the convergence of the perturbative series, including corrections to the wake potential of O(G 3 ), is left to future work.The exquisite accuracy of the leading-order wake potential relative to the nonperturbative simulations in App.C strongly suggests that the expansion parameter is no larger than ϵ ∼ Gmn ψ /k 2 0 , so that the perturbative series converges whenever the no-screening condition of Eq. 59 is satisfied.I speculate that the agreement with the leading-order wake potential prediction, despite potentially large phase accumulation across large sources, is due a combination of its phase insensitivity and its short-ranged nature in two or more spatial dimensions.
It would be amusing to obtain a positive measurement of wake forces mediated by known elementary particles.Unfortunately, there are not too many suitable light bosons to go around in the SM.The simplest possibility may be force measurements between (electrically neutral) dielectric source and target samples, to which the photon couples quadratically.Wake forces and potentials from neutrinos appear to be prohibitively tiny, no mat-ter the neutrino source (Sec.III B).An ensemble of ultracold neutrons or atoms/molecules in a cold beam may be more promising fermionic mediators of wake forces.Finally, the phenomena calculated in this work for elementary particles undoubtedly carry over to quasiparticles in condensed matter systems-indeed, they may form the basis for certain types of exotic superconductors [85], but it would be interesting to study phonons, magnons, and other collective excitations as mediators of wake forces.

ACKNOWLEDGMENTS
I am grateful for collaboration with Sergei Dubovsky in the initial stages of this project in 2018-2019; he contributed conceptual understanding to Sec.II and modelbuilding aspects of Sec.III A. Riccardo Rattazzi's questions about the treatment in Sec.II B encouraged me to come up with the simpler description in Sec.II A. I enjoyed discussing this work and Ref. [38] with Asimina Arvanitaki and Savas Dimopoulos, even though none of us fully appreciated the extent of the overlap of our work until recently.They raised concerns about the possible breakdown of the perturbative expansion, and pushed me to produce App.C. I thank Mitrajyoti Ghosh, Yuval Grossman, Walter Tangarife, Xun-Jie Xu, and Bingrong Yu for comments, encouragement, and discussions of Ref. [37].Hannah Day, Da Liu, Markus Luty, and Yue Zhao gracefully coordinated a companion paper [79] to this work that expounds on the elastic scattering and screening effects from Secs.IV B & IV D in much greater detail; they had also started considering the concepts of wake forces.Credit to Olivier Simon for pointing out the connections to Ref. [46].I also benefited from discussions with Masha Baryakhtar, Kyle Cranmer, David Dunsky, Marco Farina, Simon Foreman, Marat Freytsis, Marios Galanis, Isabel Garcia Garcia, Andrei Gruzinov, Matthew Kleban, Simon Knapen, Marius Kongsøre, Paul Langacker, Matthew Low, Duccio Pappadopulo, Olivier Simon, Philip Sørensen, Neal Weiner, Zachary Weiner, Tien-Tien Yu, and especially Junwu Huang, who joined me on (failed) attempts to make the wake force from high-energy monochromatic neutrinos observable.In this appendix, I present the derivation of the wake force for two interactions other than that of ϕ 2 ψψ in 3. The form of the monopole-monopole wake force by a complex scalar Φ is identical to that of a real scalar ϕ.This extension is (almost) trivial since a complex scalar can be decomposed into two real scalars, but it shows that the wake force does not depend on particle-antiparticle asymmetry.The second example adds a pseudoscalar interaction to the Lagrangian of Eq. 3, which leads to an additional monopole-dipole and a dipole-dipole wake force.

Complex scalar
For simplicity, I follow the steps of the classical derivation in Sec.II A (the quantized treatment gives identical results) and consider a complex scalar field Φ with Lagrangian and corresponding equation of motion for Φ, with a source particle at the origin: As before, decompose the field into a G-independent background and a spherical perturbation linear in G from a source particle at the origin: Φ = Φ 0 e −i(ωt−k0•x) + δΦ(t, x) with ω 2 = m 2 + k 2 0 .The solution to the equation of motion to leading (first) order in G is: where I show more intermediate steps than for the realscalar derivation in Sec.II A. The wake potential is the obtained by retaining the term of the potential energy V = GmΦΦ † of ψ to second order in G, which after time averaging is: This has the same form as Eq. 8 for the real scalar, noting that the energy density in the nonrelativistic limit is ρ = 2m 2 Φ 0 Φ † 0 .The quantized treatment of the complex scalar yields consistent results, and is identical to Eqs. 16 & 17 with the replacement n → (n + n)-the wake forces of particles and antiparticles are additive.
The results above are a trivial extension from the case of a real scalar, since a complex scalar can be decomposed into two real scalars, e.g.Φ = (ϕ 1 + iϕ 2 )/ √ 2 so GmΦ † Φ ψψ = (Gm/2)(ϕ 2 1 + ϕ 2 2 ) ψψ, each of which produce an identical wake force.It is nevertheless instructive that wake forces do not depend on particle-antiparticle asymmetry, a result that also holds for fermion-mediated wake forces as shown in App.B, in contrast to fermion potentials linear in G (e.g.Eq. 45).

Pseudoscalar interactions
Consider the quadratic interactions of a real scalar or pseudoscalar ϕ with both the scalar and pseudoscalar currents of ψ with couplings G s and G p , respectively.In addition to the monopole-monopole wake force proportional to G 2 s derived in Sec.II, the addition of the pseudoscalar interaction also generates a monopole-dipole wake force proportional to G s G p and a dipole-dipole one proportional to G 2 p , which will be derived in the quantized treatment below.
With the inclusion of the pseudoscalar interaction, the amplitude corresponding to the diagrams in Fig. 2 Treating the ψ fermions nonrelativistically, the identities of Eqs.D8 & D9 can be used to simplify the amplitude: Setting the energy transfer between the fermions to zero (q µ = (q 0 , q) ≃ 0), selecting the iε prescription for the retarded propagator, and using the external fourmomentum to cancel the pole (k 2 = E 2 k − k 2 = m 2 ), the propagator denominators are of the form: (q ± k) 2 − m 2 ≃ −(q 2 ± 2q • k ∓ iε).The Fourier transform of the potential is related to the amplitude as 1.Dirac fermion, scalar current The simplest example of a fermion-mediated wake force is that of a Dirac fermion χ with a scalar current ψψ: Like for the scalar field in Eq. 11, quantize the Dirac field as: with annihilation (creation) operators a k (a † k ) for particles and b k (b † k ) for particles and antiparticles, respectively.The corresponding spinor solutions u s (k) and v s (k) satisfy the standard identities of Eqs.D4 & D5.The background is taken to be a mixed state with a number density n of χ particles and a number density n of χ antiparticles with unpolarized spins and momentum distribution f (k): The total energy density in the medium χ and χ particles is then just: The Feynman diagrams for a fermionic wake force are the same as for the spin-0 equivalent of Fig. 2, with the addition of a pair of diagrams for the antiparticles χ, all shown in Fig. 8.The combined matrix element evaluates to: in the same order: where χ1 χ 1 is the interaction of the current with ψ 1 , and χ2 χ 2 with ψ 2 .Only the signs are evaluated in Eqs.B7-B10.The top right diagram of Fig. 8 has no additional minus sign, because no anticommutators are needed to perform the contractions with the vacuum for the incoming and outgoing particles, nor for the internal propagator χ χ.For the top left diagram and the contractions in Eq.B7, one needs 2 anticommutators to move χ2 to the left, then 1 to move χ 1 to the right, and finally 1 to swap χ1 and χ 2 , for an even number of anticommutators.For the bottom right diagram and Eq.B10, 1 anticommutator to bring χ 1 leftward, 1 to bring χ2 rightward, and 1 to swap χ1 and χ 2 , for an odd number of anticommutators and thus an extra minus sign.The bottom left diagram and Eq.B9 are similarly negative.This explains the extra minus sign for the bottom (antiparticle) diagrams in Fig. 8, reflected in the last line of Eq.B6.
The the matrix element is similar to that of the scalar case, with the addition of spinor numerators such as s ūs (k) / q + / k + m u s (k) = −4q • k + 8m 2 for the first amplitude.The monopole-monopole wake potential is: with the form factor appropriate for fermions: In the nonrelativistic limit |k| ≪ m, this form factor is the same as the one for the monopole-monopole wake potential mediated by scalars in Eq. 17.Similarly, the overall fermionic monopole-monopole wake potential of Eq.B12 is the same as the scalar one of Eq. 16 in the coupling normalizations of this work (Eqs.3 & B1), with fermions and antifermions contributing additively.
For nonrelativistic particles at low occupation numbers, there is thus no material difference between monopolemonopole wake forces mediated by fermions and scalars.

Four-Fermi interactions of Dirac fermions
Consider the neutral-current four-Fermi interaction of a single Dirac neutrino in Eq. 1, first for a single mass eigenstate.The matrix element for the wake potential diagrams of Fig. 8 is similar to Eq. B6, except the spinor numerators are dressed with additional gamma matrices.Firstly, the linear mass term in the numerator of the fermion propagator does not contribute because the interactions preserve chirality, i.e. (1 − γ 5 )γ ν (1 − γ 5 ) = 0. Secondly, neutrinos and antineutrinos contribute additively as in all previous cases.The full wake potential between two fermions ψ 1 and ψ 2 with vector and axial couplings g and spins σ 1 , σ 2 is thus: where the I vector and Γ tensors are defined as: Monopole-monopole wake potential-Using the identity Γ 0ρ0 = 8 2E k η ρ0 − k ρ , the monopole-monopole wake potential is found to be quite similar to that of the quadratic interaction with the scalar current in Eq.B12: with the vector-vector form factor: In the nonrelativistic limit, this form factor reduces to that of Eqs. 17 & B13.
Monopole-dipole wake potential-Without loss of generality, consider the part of the wake monopole-dipole potential proportional to g where the form factor is now a vector function: In the nonrelativistic limit and r ≲ 1/|k 0 |, this form factor generically has a magnitude of order the typical velocity |k 0 |/m of the medium.

Majorana fermions
Suppose the fermions χ are their own antiparticles, and quantize the theory with creation and annihilation operators a † k and a k : Assume the modes are populated with number density n and momentum distribution f (k) according to Eq. B3.More contractions are possible for Majorana fields, namely χχ and χ χ in addition to the usual χ χ.
For quadratic interactions of the form χΓχ with Γ = {1, iγ 5 , γ µ γ 5 }, these additional contractions can all be "absorbed" by conventional factors of 1/2 in the vertices and propagators.For example, consider the Lagrangian of a Majorana fermion interacting with a scalar current: This is the analog of Eq.B1, and will yield the identical wake potential to Eq. B12 with this normalization of G and with the replacement (n + n) → n.Vertex structures with Γ = {γ µ , σ µν } yield zero. 2+1 dimensions-The average density profile n(x) of the 2 + 1D simulations is shown as a function of x = (x, y) and |x| (azimuthal average) in Figs. 9 and 10 respectively, and is normalized relative to the average number density ⟨n⟩ x over the 2-torus.The simulation parameters are σ k = 100 and Gmn ψ (x) = {−400, +200, +400, +4,000} for |x| < R = 0.1, and zero otherwise.The Gaussian random field initialization does not strictly correspond to a pure Maxwell-Boltzmann distribution in the vacuum outside the barrier, because a fraction πR 2 of the total density is initialized inside the barrier at a slightly higher energy.I conservatively estimate the systematic error of this effect on n(x)/⟨n⟩ x to be σ sys ∼ ϵπR 2 inside the barrier, and σ sys ∼ ϵ(πR 2 ) 2 outside.These systematic deviations decrease with smaller fractional volumes occupied by the barrier.
Eq. C1 can also be solved perturbatively, using the same methods as in the main text.Decompose the field into the perturbative series Ψ = Ψ (0) + Ψ (1) + . . .and similarly for the potential of Eq.C2: V ≃ V (0) + V Defining the (unperturbed) "number density" n ≡ α |Ψ α | 2 , the expected potentials to first-and secondorder in G are: The 2D form factor is: Perturbatively, the predicted number density can be written as n(x) ≃ 2(V (0) + V (1) )/G to leading order in G.These perturbative predictions are overlaid as red G V (0) + V (1) ±1σ th thermodynamic estimate FIG. 10.Number density n(|x|) relative to the 2-torus average ⟨n⟩x in 2 + 1D as in Fig. 9 but azimuthally averaged at a distance |x| from the center of the disk-shaped barrier/well with R = 0.1, and for a MB distribution with σ k = 100.The effective coupling inside the disk is varied over Gmn ψ = {−400, +200, +400, +4,000}, with the first three displayed with vertical shifts of 0.0025 in the left panel, and the last shown separately in the right panel.Black curves show the time-averaged means over approximately 5,000 simulations, and the green (yellow) bands indicate 1σstat (2σstat) statistical errors.Red curves depict the leading-order wake potential prediction; red bands include a 1σ th estimated theory error of σ th ≃ ϵ 2 ≡ (Gmn ψ /2σ 2 k ) 2 .The thermodynamic estimate from Eq. C9 is included in dashed light blue.The systematic error from partial field initialization inside the barrier is estimated to be σsys ∼ ϵπR 2 ∼ 0.006 for the right panel (and negligible in the left panel) but is not shown.solid lines on Fig. 10, where they are found to match the numerical results well when the perturbative expansion parameter ϵ ≡ |V (1) /V (0) | = |G|mn ψ /2σ 2 k , which is the same in 2 + 1D and 3 + 1D, is small in absolute value.The red bands in Fig. 10 are the corresponding fractional ±1σ th theory errors around the first-order prediction with σ th = ϵ 2 .The permille-level concordance between the numerical results and the perturbative prediction indicates that the perturbative expansion parameter is ϵ in two or more spatial dimensions, even for large sources.
Finally, a simple thermodynamic argument can be used to estimate the density modulation deep inside a medium-dependent potential.Suppose particles outside the potential have a number density n outside ∝ d ¯dk e −k 2 /2σ 2 k = d ¯dk e −E kin /T given by a Maxwell-Boltzmann distribution with effective temperature T = σ 2 k /m.Modes deep inside a constant potential should have an effective Boltzmann factor e −(E kin +Gn ψ /2)/T , thus yielding the nonperturbative prediction: assuming such equilibrium is indeed established.In the perturbative regime, the estimate agrees with that of the first-order wake potential: n inside /n outside − 1 = V (1) /V (0) = −Gmn ψ /2σ 2 k , which holds both in 2 + 1D and 3 + 1D.The nonperturbative thermodynamic estimate from Eq. C9 is indicated by the light blue dashed line in the right panel of Fig. 10.It yields a milder density decrease because the next-to-leading order in the perturbative expansion has opposite sign.
3+1 dimensions-The results of the 3 + 1D simulations are shown in Fig. 11 as black curves with green (yellow) statistical error bands at ±1σ stat (±2σ stat ).The chosen numerical resolution is lower for computational feasibility reasons, but the qualitative features are similar to the 2 + 1D case.The perturbative prediction for the first-order nonrelativistic wake potential was already derived in the main text (cfr.Eqs.16 and 19): The corresponding perturbative predictions for the density of ϕ particles, n(x) = (2/G)(Gn/2 + V (1) (x)) are shown as red curves with 1σ th = ϵ 2 theory error bands.The thermodynamic estimate from Eq. C9 is indicated by the light blue dashed line, and lines up precisely with the numerical results deep inside the barrier.The systematic uncertainty from partial field initialization inside the barrier is not shown, but should be of order σ sys ∼ ϵ(4π/3)R 3 for |x| < R, and is negligible for R = 0.1 but more appreciable for R = 0.2, with the hierarchy σ stat ≲ σ sys ≲ σ th in that case.Like in 2 + 1D, the nonperturbative numerical results in 3 + 1D are in quantitative agreement with perturbative wake potential calculations.10.The estimated systematic error from partial field initialization inside the barrier is σsys ∼ ϵ(4π/3)R 3 ∼ {0.0004, 0.003} for |x| < R but is not shown.

FIG. 2 .
FIG.2.Diagrams contributing to the wake force between distinguishable, nonrelativistic fermions ψ1 and ψ2 in a medium of real scalars ϕ.

FIG. 5 .
FIG. 5. Distance (r) dependence of the form factor magnitude |F(x)| for MB (top panel) and FD (bottom panel) momentum distributions.The black curves are the isotropic form factors of Eqs.19 & 21 (no boost).The red (gold) curves are for BMB and BFD distributions of Eqs.23 & 24 in the forward (backward) direction.The other colored curves (green, blue, purple) are for small mass-squared splittings ∆ 2 > 0, in the forward direction for the BMB distribution, and for the isotropic FD distribution.The boost/anisotropy alone lifts the form factor suppression to F ∝ 1/r 2 in all cases, while for small mass splittings, F ∼ ∆ 2 /|k0| 2 for r ≲ λosc ∼ |k0|/|∆ 2 | can be achieved, where |k0| is the typical momentum in the distribution.

− 2 ]
DM wake = gravity G s,N = G F s u p e r fl u i d H e σ e l = 1 0 − 4 0 c m 2 σ e l = 1 0 − 5 0 c m 2 s c r e e n e d w a k e w a k e f o r c e ( e x i s t i n g ) w a k e f o r c e ( T B , T = 1 0 m K ) w a k e f o r c e ( T B , S Q L ) r u c t u r e fo r m a t io n ↓ ↓ B B N ↓ DM Monopole-Monopole Wake Force k [m]
This work is supported by the National Science Foundation under Grant PHY-2210551.The Center for Computational Astrophysics at the Flatiron Institute is supported by the Simons Foundation.This work was performed in part at the Institute for Nuclear Theory at the University of Washington, supported in part by the U.S. Department of Energy grant No. DE-FG02-00ER41132; and the Aspen Center for Physics, supported by National Science Foundation grants PHY-1607611 and PHY-2210452.

FIG. 8 .
FIG.8.Diagrams contributing to the wake force in a medium of Dirac fermions χ and antifermions χ.

n = 3 ,FIG. 9 .
FIG. 9. Time-averaged number density n(x) = |Ψ(x)| 2 over a 2-torus in 2+1D simulations of nonrelativistic ϕ particles with isotropic Maxwell-Boltzmann distribution parametrized by a momentum spread σ k = 100.The particles are scattered by a uniform disk of radius R = 0.1 (dashed yellow outline) with an effective coupling strength Gmn ψ = +200 in its interior.The density profile is normalized to its 2-torus average ⟨n⟩x.Upon taking the mean over many coherence times and simulations, the fractional density profile has permille-level statistical density variations, revealing the fractional density modulation of order −Gmn ψ /2σ 2 k = −0.01due to the wake of the disk.The coherence length of the density variation/modulation is of O(σ −1 k ).
)J 0 (k|x|), (C8) with the last line for a 2D Maxwell-Boltzmann distribution.It scales logarithmically at short distances, and falls off exponentially for |x| ≫ σ −1 k , so that the 2 + 1D (like the 3 + 1D) wake force has a short range.
FIG.11.Number density n(|x|) relative to the 3-torus average ⟨n⟩x in the 3+1D simulations, averaged at a distance |x| from the center of a ball-shaped barrier with strength Gmn ψ = +80 and radius R = {0.1,0.2} (the former vertically offset by +0.04 for clarity), for an isotropic MB distribution with σ k = 20.Labels are as in Fig.10.The estimated systematic error from partial field initialization inside the barrier is σsys ∼ ϵ(4π/3)R 3 ∼ {0.0004, 0.003} for |x| < R but is not shown.