How do scalar-field dark matter haloes react to orbiting bodies?

TL;DR

This paper investigates how scalar-field dark matter halos respond to orbiting masses, revealing the formation of rotating structures and quantifying dissipation effects, which enhances understanding of dark matter interactions.

Contribution

It introduces a detailed analysis of the dynamical response of scalar-field dark matter halos to orbiting bodies, including the formation of non-spherical structures and dissipation mechanisms.

Findings

Orbiting mass induces rotating non-spherical clumps.

Dissipation force scales with the square of the mass.

Scalar field structure reacts dynamically to perturbations.

Abstract

Low-energy, self-gravitating solutions of a scalar field coupled to gravity, described by the Schrodinger-Poisson system, are good candidates for realistic astrophysical structures, being particularly suited to describe dark matter halos. In this work we study the scenario in which one of these structures is gravitationally perturbed by a point-like mass. We analyse the effects that the body has on the distribution of the scalar field and how it backreacts on the body's motion. We show that an initially static, spherical structure can develop rotating non-spherical clumps, the amplitude and the velocity of which are directly related to the mass of the orbiting particle. We also study the dissipation mechanisms involved in the transit of the point-like particle across the scalar field structure and we observe that the force responsible for the dissipation scales as the square of the mass of the particle.

Figures

Equations

$$S=\int d^{4}x\sqrt{-g}\left(\frac{R}{16\pi }-\frac{1}{2}{g}^{\alpha \beta }{\partial }_{\alpha }{\Phi }^{\ast }{\partial }_{\beta }\Phi -\frac{1}{2}{m}_s^2\Phi {\Phi }^{\ast }\right),$$

$$R_{\alpha \beta }-\frac{R}{2}{g}_{\alpha \beta }=8\pi T_{\alpha \beta },$$

$$T_{\alpha \beta }=\frac{1}{2}\left({\Phi }_{,\alpha }{\Phi }_{,\beta }^{\ast }+{\Phi }_{,\alpha }^{\ast }{\Phi }_{,\beta }\right)-\frac{g_{\alpha \beta }}{2}\left({\Phi }^{,\gamma }{\Phi }_{,\sigma }^{\ast }+m_s^2\Phi {\Phi }^{\ast }\right),$$

$$\frac{1}{\sqrt{-g}}{\partial }_{\alpha }\left(\sqrt{-g}{g}^{\alpha \beta }{\partial }_{\beta }\Phi \right)=m_s^2\Phi ,$$

$$g_{00}=-1+2U(r),g_{jk}=\left[1+2U(r)\right]{\delta }_{jk},$$

$$\Phi \sim \exp (-\mathrm{i}{m}_{s}t)\chi (t,r),$$

$$\mathrm{i}{\partial }_t\chi$$

$${\nabla }^{2}U$$

$$t\to m_{s}t,r\to m_{s}r,\chi \to \frac{1}{\sqrt{4\pi }}\psi ,$$

$$\mathrm{i}{\partial }_t\psi$$

$${\nabla }^{2}U$$

$$(t,r,U,\psi )\to ({\lambda }^{-2}\hat{t},{\lambda }^{-1}\hat{r},{\lambda }^2\hat{U},{\lambda }^2\hat{\psi }),$$

$${\lambda }^2<10^{-3}.$$

$$\psi =\exp (-\mathrm{i}\gamma t)f(r),$$

$$f^{\prime \prime }(r)+\frac{2}{r}{f}^{\prime }(r)+2(\gamma -U(r))f(r)=0,$$

$$U^{\prime \prime }(r)+\frac{2}{r}{U}^{\prime }(r)=f(r)f^{\ast }(r).$$

$$M_{\mathrm{DM}}=\int \rho (t,r)dV,$$

$$M_f=\int {\rho }_f(t,r)dV=\frac{1}{4\pi }\int \psi {\psi }^{\ast }dV,$$

$$M_f=\frac{1}{4\pi }{\int }_0^{\infty }{\psi }^{\ast }\psi dV={\int }_0^{\infty }f(r)f^{\ast }(r)r^{2}dr.$$

$$\underset{r\to \infty }{\lim }U(r)=-\frac{M_f}{r}.$$

$$\gamma$$

$$M_f$$

$$T_P^{\mu \nu }=\frac{1}{\sqrt{-g}}{M}_p\int \frac{d{x}_P^{\mu }}{d\tau }\frac{d{x}_P^{\nu }}{d\tau }{\delta }^{(4)}(x^{\alpha }-x_P^{\alpha }(\tau ))\mathrm{d}\tau ,$$

$$\mathrm{i}{\partial }_t\psi$$

$${\nabla }^{2}U$$

$$P(\mathbf{x},{\mathbf{x}}_P)=\frac{4\pi M_P}{r^2}\delta (r-r_P)\delta (\cos \theta -\cos {\theta }_P)\delta (\varphi -{\varphi }_P),$$

$$\delta (\cos \theta -\cos {\theta }_P)\delta (\varphi -{\varphi }_P)=\sum _{\ell m}{Y}_{\ell m}^{\ast }({\theta }_P,{\varphi }_P)Y_{\ell m}(\theta ,\varphi ),$$

$$P(\mathbf{x},{\mathbf{x}}_P)=p_A(r,{\mathbf{x}}_P)+\sum _{\ell ,m}{p}_{\ell ,m}(r,{\mathbf{x}}_P)Y_{\ell m}(\theta ,\varphi ),$$

$$p_A(r,{\mathbf{x}}_P)$$

$$p_{\ell ,m}(r,{\mathbf{x}}_P)$$

$$\psi (r,\theta ,\varphi )$$

$$U(r,\theta ,\varphi )$$

$$\{\begin{array}{ll}{\partial }_t{\phi }_A& =\frac{\mathrm{i}}{2}{\nabla }_r^2{\phi }_A-\mathrm{i}\int \mathrm{d}\Omega [U\psi ]\\ {\nabla }_r^2{V}_A& =p_A+\int \mathrm{d}\Omega [\psi {\psi }^{\ast }]\end{array},$$

$$\{\begin{array}{l}{\partial }_t{\phi }_{\ell ,m}=\frac{\mathrm{i}}{2}\left({\nabla }_r^2-\frac{\ell (\ell +1)}{r^2}\right){\phi }_{\ell ,m}-\mathrm{i}\int \mathrm{d}\Omega [U\psi ]Y_{\ell ,m}^{\ast }\\ \left({\nabla }_r^2-\frac{\ell (\ell +1)}{r^2}\right)V_{\ell ,m}=p_{\ell ,m}+\int \mathrm{d}\Omega [\psi {\psi }^{\ast }]Y_{\ell ,m}^{\ast }\end{array}.$$

$${\phi }_A(t=0)$$

Full text

1. I Introduction

2. II Framework

3. II.1 Low energy scalar field structures

4. II.2 The point-like particle

5. II.3 The evolution equations

6. III Running the simulation

7. III.1 The initial conditions

8. III.2 The boundary conditions of the SP-system

9. III.3 Time and space discretization of the system

10. IV Results

11. IV.1 General evolution of the field

12. IV.2 Effects on the orbiting particle

13. IV.2.1 Friction force

14. IV.2.2 Loss of angular momentum

15. IV.3 Changes in the density distribution of the field

16. IV.3.1 Time dependence of the non-spherical density

17. IV.3.2 Magnitude of the non-spherical density

18. V Conclusion

19. A Testing the code

20. A.1 Evolving a stationary scalar field solution

21. A.2 Testing the code evolving the orbiting particle

22. B Discretizing the Dirac delta

23. C Derivation of the Schrodinger-Poisson system with orbiting particle

24. C.1 Poisson equation

25. C.2 Klein-Gordon equation

I Introduction

Scalar fields are the simplest objects one can have in a field theory. Their simplicity is the key for their versatility which puts them in the center of the most important debates of modern physics. Scalar fields are present in many areas of study, either as by-products in low energy limits of string theory Svrcek and Witten (2006); Arvanitaki et al. (2010) or as postulated ingredients of models that try to explain poorly-understood phenomena, such as inflation Nanopoulos et al. (1983); Baer et al. (2015), quintessence Ratra and Peebles (1988); Caldwell et al. (1998) or the QCD’s strong CP problem Peccei and Quinn (1977).

The plethora of scalar fields that appear in these scenarios are predicted in such a way that they interact very weakly with baryonic matter, making them hard to reveal in particle detectors. In any case, all of these fields have to interact gravitationally and it is at astrophysical scales that their gravitational influence may become big enough for their presence to be detected Kesden et al. (2005); Macedo et al. (2013a); Khmelnitsky and Rubakov (2014); Diemer et al. (2013); Macedo et al. (2013b); Blas et al. (2017); Ferreira et al. (2017); Bošković et al. (2018); Boskovic et al. (2018).

One case of interest happens when scalar fields develop structures – dubbed “scalar clouds” – around black holes (BHs) as a result of the activation of a superradiance mechanism Brito et al. (2015a). These clouds can have a strong impact on the gravitational wave-driven dynamics of BHs, creating a fertile ground for astrophysical tests Brito et al. (2015b); Arvanitaki et al. (2015); Herdeiro and Radu (2015); Cunha et al. (2015); Ferreira et al. (2017); Brito et al. (2017a, b); Ikeda et al. (2018). Another astrophysically relevant scenario involving scalar fields is the formation of self-gravitating structures known as boson stars or oscillatons, for complex and real scalars, respectively Kaup (1968); Ruffini and Bonazzola (1969); Brito et al. (2015c, 2016); Seidel and Suen (1991); Liebling and Palenzuela (2012). These scalar structures are solutions of the Einstein-Klein-Gordon system of equations and their characteristics can vary substantially depending on the scalar field potential that is considered. Indeed, different potentials give rise to structures that can either resemble BHs Guzman and Rueda-Becerril (2009); Palenzuela et al. (2017) or big galactic haloes Urena-Lopez and Bernal (2010); Lee and Koh (1996). The latter case has motivated the use of scalar fields as a viable model to explain the phenomena atributed to dark matter Mielke and Schunck (2002); Hui et al. (2017); Hu et al. (2000). The formation of boson stars and oscillatons has been studied Seidel and Suen (1994); Madsen and Liddle (1990); Khlopov et al. (1985); Guzman and Urena-Lopez (2006); Okawa et al. (2014); Hidalgo et al. (2017); Veltmaat et al. (2018)) as well as their stability Jetzer (1989); Gleiser (1988); Lee and Pang (1989); Seidel and Suen (1990); Balakrishna et al. (1998); Guzman and Urena-Lopez (2004); Grandclement et al. (2011); Alcubierre et al. (2003); Brito et al. (2017a, b). They can exist in a ground state (a 0-node solution) or they can be formed in an excited state (a $n$-node solution, $n\neq 0$) which either decays to the ground state, collapses to a BH or dissipates away Guzman and Urena-Lopez (2006, 2004); Balakrishna et al. (2008); Guzman (2004).

Infering the existence of scalar fields from its self-gravitating structures requires that one knows what observational imprints those structures may leave. This line of research is a very active one, some exemples of which are the analysis of geodesics around these structures Diemer et al. (2013); Macedo et al. (2013b); Kesden et al. (2005); Macedo et al. (2013a); Brihaye et al. (2014); Grandclement et al. (2014); Becerril et al. (2006); Bošković et al. (2018), the study of their gravitational lensing effects Schunck et al. (2006); Dabrowski and Schunck (2000), their emission of gravitational radiation Balakrishna et al. (2006); Ferrell and Gleiser (1989); Kesden et al. (2005), its behavior when involved in collisions with other structures Eby et al. (2017); Palenzuela et al. (2007, 2017); Clough et al. (2018) and their effects on the pulsar timing signal measurements Khmelnitsky and Rubakov (2014); Blas et al. (2017).

In the context of scalar field dark matter haloes Mielke and Schunck (2002); Hui et al. (2017); Hu et al. (2000), only the self-gravitating structures which are characterized by low compacteness and large radii are of interest. Boson stars and oscillatons which comply with these requirements can be studied within a Newtonian setup. In fact, when the magnitude of the scalar field is very small, the solutions of the Einstein-Klein-Gordon system have extremely large spatial extent and the underlying spacetime geometry is very close to Minkowski. It was shown Ferrell and Gleiser (1989); Guzman and Urena-Lopez (2004); Bošković et al. (2018) that to study this regime, it is sufficient to work with the simpler Schrodinger-Poisson (SP) system instead of the Einstein-Klein-Gordon one. Moreover, the SP system describes the dynamics of both complex and real scalar fields Guzman and Urena-Lopez (2004) in the low-energy limit.

Since the low-energy scalar configurations in the ground state are stable, one can assume that they may be astrophysically relevant as components of bounded gravitational systems. We are particularly interested in a two-body system in which one of the bodies is described by a stable scalar field configuration and the other is pointlike. In this scenario, the scalar field configuration is tidally deformed due to the presence of the point-particle and that deformation affects the evolution of the whole system. We developed a numerical code to study this two-body scenario and we are able to describe the tidal deformations of the scalar field density as well as the gravitational friction force that the point-particle feels as it traverses the scalar configuration.

We use Planck units – $c=\hbar=G=1$ – unless otherwise stated.

II Framework

We consider the scenario in which a low-energy, stable, scalar field configuration of the Einstein-Klein Gordon system is perturbed by a pointlike mass. This two-component system will be evolved separately and the only interaction between the components is gravitational. The time evolution of the low-energy scalar struture will be described by the SP equations, which will contain the effect of the presence of the point-particle, whereas the movement of the point-particle will be rendered from the gravitational potential of the scalar field configuration by the laws of Newtonian mechanics.

II.1 Low energy scalar field structures

A complex scalar field of mass $m_{s}$, minimally coupled to gravity, is described by the action

[TABLE]

Minimizing this action is equivalent to solving the Einstein-Klein-Gordon system

[TABLE]

where $R_{\alpha\beta}$ and $R$ represent the Ricci tensor and Ricci scalar, respectively, $g_{\alpha\beta}$ is the metric tensor and $g$ is its determinant. We are interested in the low-energy limit, i.e., in situations in which the scalar field energy is given almost entirely by its rest-energy and the space-time metric deviates slightly from the Minkowski metric. In the harmonic gauge, $\partial_{\nu}\left(\sqrt{-g}g^{\mu\nu}\right)=0$, the metric tensor can be written as

[TABLE]

and $\quad g_{0j}=0$. The scalar field has an harmonic time dependence

[TABLE]

with $\chi(t,r)$ a slowly varying function of time. Notice that the spatial dependence is solely radial, since we’re assuming spherical symmetry. Substituting Eqs. (5) and (6) in the Einstein-Klein-Gordon system and focusing only on the leading order terms, one obtains the Schrodinger-Poisson (SP) system

[TABLE]

Making the following re-scaling

[TABLE]

the previous equations can be written as

[TABLE]

It is known Guzman and Urena-Lopez (2004) that a transformation of the form

[TABLE]

leaves the SP system of equations unchanged. Using this property, one can normalize the system by working only with the “hat-variables”, meaning that the order of magnitude of all the quantities involved in this problem is hidden in the parameter $\lambda$. Notice that Ref. Urena-Lopez et al. (2002) fixes the value of $\lambda$ for which the behavior described by the SP sytem coincides with the behavior of the Einstein-Klein-Gordon system (i.e. the scale of the low-energy limit)

[TABLE]

This is the limit of validity of all the statements that stem from the analysis of the SP system: these calculations only cover scalar fields whose magnitude is compatible with the previous limit. So, having established the limits of validity of our working system, we will, from now on and unless otherwise states, work in terms of the “hat-quantities” of Eq. (10).

Stationary solutions for the SP system have been found Guzman and Urena-Lopez (2004) and their stability was studied. These solutions are spherically symmetric and are classified by the number of nodes they present, i.e., by the number of times that the scalar field function changes sign. It was found that only the fundamental (0-node) solutions are stable; all the other stationary solutions decay into nodeless configurations by losing mass through scalar radiation Guzman and Urena-Lopez (2004). To obtain these stationary solutions, we consider a scalar field of the form

[TABLE]

i.e. $\gamma$ is the difference between the total energy of the field and its rest energy (notice that $\psi$ corresponds to the slow-varying part of the scalar field, the total scalar field function being $\Phi\sim\exp(-\mathrm{i}m_{s}t)\exp(-\mathrm{i}\gamma t)f(r)$; thus, its energy is $E=m_{s}+\gamma$)111Since we are in a low-energy limit, this difference is very small, however, since we are working with the “hat-quantities” (see Eq. (10)), the value of $\gamma$ that we present in the equations is related with that normalized system of coordinates; to convert it back to Planck units, one has to multiply it by $\lambda^{2}$, which will consistently make it small and since stationary configurations are bound states of the system, one expects222Using the non-scaled variables of Eq. (7), one can see that since $\frac{\gamma}{m_{s}}=\frac{\omega}{m_{s}}-1$, and because the bound state condition is $\frac{\omega}{m_{s}}<1$, a bound state must have $\frac{\gamma}{m_{s}}<0$. that $\gamma<0$. By substituting the previous ansatz on the SP system, one obtains

[TABLE]

These equations are used to find the profiles $f(r)$ of the stationary configurations of the scalar field. To do it, one has to impose boundary conditions that come from two reasonable physical requirements: the profile has to be regular and finite. Regularity is enforced by demanding that $f^{\prime}(r)$ and $U^{\prime}(r)$ are zero in the origin; finiteness is then guaranteed by insisting that $\lim_{r\rightarrow\infty}f(r)=0$. Moreover, one must demand that the resulting gravitational potential $U(r)$, when measured at infinity, describes, as it should, the effect of the total mass of the scalar field configuration. This quantity is calculated with the volume integral

[TABLE]

where $\rho=m_{s}^{2}\Phi\Phi^{*}$ is the leading order term of the the weak-field limit of the $00$ component of the scalar field energy-momentum tensor. Writing the previous integral in terms of the hat-quantities of Eq. (10), we obtain $M_{\rm DM}=M_{f}/m_{s}$ where

[TABLE]

which in the case of the stationary configuration, can be simplified to

[TABLE]

Having calculated the mass, one can write the boundary behavior for the gravitational potential in the hat-quantities as

[TABLE]

Imposing these conditions along with the value for the scalar field at the origin (since we are working with the normalized “hat quantities” it is enough to consider $f(0)=1$), one obtains profiles $f(r)$ which can be characterized by the number of nodes. To each of the profiles corresponds a unique value of the quantities $\gamma$ and $U(r=0)$. In Fig. 1 we show two stationary configurations and their respective potential. Since we are only interested in the 0-node, ground state solutions, we quote here only its characteristic values,

[TABLE]

where $R_{99}$ is the radial position up to which $99\%$ of the mass of the scalar configuration is contained.

II.2 The point-like particle

We now want to understand how an orbiting mass $M_{p}$ disturbs, dynamically, the previous self-gravitating massive scalar structure. We model the orbiting mass as pointlike 333Due to the nonlinear nature of Einstein’s field equations, problems arise in the definition of point particles in such context. However, we will always be working in the Newtonian limit where such idealization is acceptable. and use its energy-momentum tensor in the Einstein-Klein-Gordon system of Eqs. (4) (see section 6.5 of Ref. Straumann (2012))

[TABLE]

where $x^{\alpha}_{P}$ are the spacetime coordinates of the point-particle, $\tau$ is its proper-time and $\delta^{(4)}$ is the Dirac-delta. A low-energy analysis of the Einstein-Klein-Gordon system with the pointlike mass included yields (see Appendix C)

[TABLE]

with

[TABLE]

where $r_{P}$, $\theta_{P}$ and $\phi_{P}$ correspond to the spherical coordinates indicating the position of the point-like particle. Notice that since we are using the “hat quantities” of Eq. (10), the value of the mass of the orbiting particle, $M_{p}$, when converted to Planck units, must also be multiplied by $m_{s}$. Using the spherical harmonics closure relation Arfken et al. (2011)

[TABLE]

Eq. (21) can be re-written as a sum of spherical harmonics

[TABLE]

where

[TABLE]

This expansion motivates us to write the other quantities involved in the problem in a similar way, i.e.

[TABLE]

II.3 The evolution equations

When the pointlike particle is put in orbit, both the scalar field and the gravitational potential will, in general, develop a multipolar strucuture which will translate in the development of non-trivial profiles for the $(\ell,m)$ components of the expansions of Eqs. (24) and (25). This development will be coordinated by the Schrodinger-Poisson system of equations, where we introduce the aforementioned expressions for the scalar field and gravitational potential, which are then projected in each $(\ell,m)$ component. By doing that, we obtain two equations for each mode, one for $\varphi_{\ell,m}$, coming from the projection of the Schrodinger equation, and one for $V_{\ell,m}$ coming from the projection of the Poisson equation. Particularly, in order to harvest the zeroth order component (which we indicate with “A”) we integrate both sides of it over the whole sphere and we obtain

[TABLE]

where $d\Omega=\sin\theta d\phi d\theta$ and $\nabla^{2}_{r}\equiv\frac{\partial^{2}}{\partial r^{2}}+\frac{2}{r}\frac{\partial}{\partial r}$. To obtain the equations corresponding to a general $(\ell,m)$ mode we integrate each side of the equations multiplied by the corresponding spherical harmonic function and we obtain

[TABLE]

III Running the simulation

In order to study the evolution of our two body system we evolve the SP-equations for each $(\ell,m)$ mode to account for the scalar field structure, along with the equations of motion of the point particle, which are given by Newton’s laws of motion. This is accomplished with the following algorithm:

Given the initial conditions, calculate the gravitational potential of the “scalar field + particle” system using a sparse matrix solver included in SciPy Jones et al. (01); 2. 2.

From $t=0$ to $t=\Delta t$, use the iterated Crank-Nicolson method Teukolsky (2000) to advance the scalar field and the Euler’s method Hairer et al. (2008) to advance the position and velocity of the particle; 3. 3.

Calculate the gravitational potential of the “scalar field + particle” system using a sparse matrix solver; 4. 4.

From $t=\Delta t$ to $t=2\Delta t$, use the iterated Crank-Nicolson method to advance the scalar field and the two-step Adams-Bashforth Hairer et al. (2008) method to advance the position and velocity of the particle; 5. 5.

Repeat the former and the latter steps until $t=700$.

III.1 The initial conditions

At $t=0$, the scalar field is given by the groundstate, 0-node, stable configuration that was calculated in Section II.1. The point particle’s degrees of freedom are its initial position and velocity. We performed several simulations to study the influence of varying these parameters on the evolution of the system. In each simulation, the initial condition for the field and the gravitational potential are

[TABLE]

where $f_{E}$ and $U_{E}$ are given by the 0-node solutions of Eq. (13).

The point particle will be thrown at the scalar self-gravitating structure, a setup characterized by (see Fig. 2),

the impact parameter $y_{i}$, 2. 2.

the mass $M_{P}$ of the particle being thrown, 3. 3.

the velocity $v_{i}$ with which the particle is thrown.

To write the initial conditions of this motion, we consider the plane that contains the position and velocity vectors of the particle (for all purposes it can be the $\theta=\pi/2$ plane). In this plane, we put the center of coordinates in the center of the scalar field structure, we use $(x_{i},y_{i})$ to indicate the initial position of the perturbing particle and $(-v_{i},0)$ to indicate its initial velocity. Then, we can obtain the initial conditions in polar coordinates of the plane (see Fig. 2):

[TABLE]

We run 27 ($3\times 3\times 3$) simulations, spanning the following set of initial conditions

[TABLE]

To make it easier to refer to each simulation, we assign a unique code to each of them. The first character of the code refers to the mass of the particle – “L”, “M” or “S”, i.e., large, medium or small – for, respectively, $M_{p}=0.1,10^{-3}$ or $10^{-5}$; the remaining characters will indicate explicitly the impact parameter $y_{i}$ and the initial velocity $v_{i}$. As an example, the simulation that has as initial conditions $y_{i}=1.0,v_{i}=0.3,M_{P}=0.1$ is called “simulation LY1V03”.

We also run 9 ($3\times 3$) simulations in which we evolve the equations of motion of the particle and the SP-system side-by-side without considering the effect of the particle on the scalar field. We will call these the “control tests” and we will refer to them with a similar code indicating the impact parameter and the initial velocity. So, the “control test” that has as initial conditions $y_{i}=1.0,v_{i}=0.3$ is referred to by “control test Y1V03”.

III.2 The boundary conditions of the SP-system

Regarding the scalar field components, and in line to the description in Section II.1, we demand regularity at the boundaries of all quantities. At the origin, regularity is guaranteed by fixing $\varphi^{\prime}_{A}(0)=0$, $\varphi_{\ell,m}(0)=\varphi^{\prime}_{\ell,m}(0)=0$ and $U^{\prime}(0)=0$. To treat the boundary condition at infinity, we benefit from the careful analysis made in Ref. Guzman and Urena-Lopez (2004). The authors show that in order to treat spatial infinity in this system, one must either put it far enough from the active zone or add a sponge to the simulation so that no reflections at the infinity boundary occur. We choose the first option. All the runs of our code were made with a spatial grid that extends up to $r=1000$ (in the agreed units). We consider that at the infinity boundary all components of the scalar field and of the potential are zero except the spherical component of the potential, $V_{A}$. In fact, as the mass of the scalar field structure has to be conserved in the grid, we impose that $V_{A}(r=1000)=-M_{f}/1000$, where $M_{f}$ is given in Eq. (18).

III.3 Time and space discretization of the system

We tested the code – with results that can be seen in Appendix A – and based on that we decided to conduct all the simulations using a spatial grid of $\Delta r=0.1$ and time step of $\Delta t=10^{-3}$.

IV Results

In the simulations of the “scalar field + particle” system, the evolution of both components encodes information about the whole system. We will show details regarding the movement of the particle and the struture of the scalar field configuration as a function of time. Particularly, we will analyse how the backreactions of the field affect the movement of the particle and how the non-spherical components of the field evolve.

IV.1 General evolution of the field

We verified that in the simulations we ran, the description of all the quantities involved – the scalar field, the gravitational potential, and the trajectory of the point-particle – were dominated by the $\ell=0$ and $\ell=1$ terms of the expansions in Eqs. (24) and (25). Particularly, we verify that for the cases with particle mass given by $M_{p}=10^{-3}$ and $M_{p}=10^{-5}$ the terms $\ell\geq 2$ are completely negligible, whereas for the case with $M_{p}=0.1$, we verify that $\mathrm{max}[|\varphi_{2,m}|/|\varphi_{1,m}|]\sim\mathrm{max}[|V_{2,m}|/|V_{1,m}|]\sim\mathcal{O}(10^{-2})$, which is the upper bound for any other ratio of the form $\mathrm{max}[|\varphi_{\ell+1,m}|/|\varphi_{\ell,m}|]$ or $\mathrm{max}[|V_{\ell+1,m}|/|V_{\ell,m}|]$, for $\ell>1$ in the $M_{p}=0.1$ case. The latter fact is translated in a slight change in the numerical values of some of the quantities that are calculated in what follows. However, since our focus will be on orders of magnitude and not in exact numerical values, we will, for the sake of simplicity and economy in the length of the expressions, use throughout this section a truncated series to describe the meaningful quantities, i.e.,

[TABLE]

Notice that the term $m=0$ isn’t considered in this expansion. This is due to the fact that the orbital plane is taken to be $\theta=\pi/2$ (see Section III.1), and so the $m=0$ components are identically zero.

IV.2 Effects on the orbiting particle

The set of initial conditions of the orbiting particle (see Section III.1), gives rise to bounded and unbounded orbits of the equilibrium scalar field structure. Technically, an unbounded orbit has energy per unit mass,

[TABLE]

larger than zero. In such case the expression

[TABLE]

is always non-negative, where $U_{\text{eff}}$ is given by

[TABLE]

for initial angular momentum per unit mass444This classification is valid for the equilibrium configuration, in which angular momentum is conserved. $\mathcal{L}$. In our case, however, whenever we say that an orbit is unbounded, we simply mean that the apocenter of the orbit wasn’t observed in the grid during the whole simulation time. In this sense, we find that all simulations result in bounded orbits except MY1V07, SY1V07, Y1V07, MY3V07, SY3V07, Y3V07, LY5V07, MY5V07, SY5V07 and Y5V07. Notice that all the unbounded orbits correspond to initial velocity $v_{i}=0.7$ and that for the cases $y_{i}=1.0,3.0$ the control test is unbounded whereas the corresponding simulation with $M_{P}=0.1$ is bounded. These cases show that the reaction of the scalar field to the presence of the massive particle alters significantly its trajectory.

IV.2.1 Friction force

The backreaction of the scalar configuration on the motion of the particle can be computed through the calculation of the effective force that appears in the movement of the particle as it travels through the scalar cloud. We call this a “friction force”. To calculate it, we are going to compare the acceleration vector in the orbital plane of the simulations with backreactions (indicated by the subscript “sim”) with the respective “control tests”555Remember the classification introduced in subsection III.1 (indicated by the subscript “control”), i.e., our friction force is written as

[TABLE]

with $\bm{F}_{\text{control}}=M_{p}\bm{a}_{\text{control}}$ (same thing for “sim” component); the acceleration vector is written as

[TABLE]

So, we will write

[TABLE]

Another way of quantifying this force is through

[TABLE]

where $U$ represents the gravitational potential. The difference between the two should reflect the extra force that appears as a finite mass effect. The gravitational potential in Eq. (34) can be written as666We verify that the source terms of the Poisson equation for the $\ell=1$ components of the gravitational potential – see Eq. (25) –, given by $s_{1,\pm 1}=\int\int(\psi\psi^{*}+p_{1,\pm 1}Y_{1,\pm 1})Y_{1,\pm 1}^{*}\sin\theta\mathrm{d}\theta\mathrm{d}\phi$ where $p_{1,\pm 1}$ comes from the point-particle contribution (see Eqs. (20) and (22)), satisfy the following identity $(s_{1,-1}+s_{1,-1}^{*})$ = $-(s_{1,1}+s_{1,1}^{*})$ and $(s_{1,-1}-s_{1,-1}^{*})$ = $(s_{1,1}-s_{1,1}^{*})$, which means that $\mathrm{Re}[s_{1,1}]=-\mathrm{Re}[s_{1,-1}]$ and $\mathrm{Im}[s_{1,1}]=\mathrm{Im}[s_{1,-1}]$. This relation between the real and imaginary parts of the source terms allows us to write $V_{1,-1}=R(r)+\mathrm{i}I(r)$ and $V_{1,1}=-V_{1,-1}^{*}$.

[TABLE]

where $R\equiv\mathrm{Re}[V_{1,-1}]$ and $I\equiv\mathrm{Im}[V_{1,-1}]$. In the plane of motion, which we consider to be $\theta=\pi/2$, we can write $\bm{f}_{f}$ as

[TABLE]

We will calculate the quantity

[TABLE]

where $t_{\text{in}},t_{\text{out}}$ are, respectively, the time in which the particle penetrated and left the scalar field structure for the first time. In order to understand the influence of each of the parameters on the movement of the particle, we run tests in groups of simulations where only one of the parameters is changing:

same mass and initial velocity: MV1 = {LY1V07, LY3V07, LY5V07}, MV2 = {MY1V05, MY3V05, MY5V05}; 2. 2.

same mass and impact parameter: MI1 = {LY1V03, LY1V05, LY1V07}, MI2 = {MY1V03, MY1V05, MY1V07}; 3. 3.

same velocity and impact parameter: VI1 = {LY1V07, MY1V07, SY1V07}, VI2 = {LY1V05, MY1V05, SY1V05}.

One can make two comments on the results of this exercise. The first is that the mass is the factor that influences the most the magnitude of the friction force777This is not surprising in our context given that the range of variation of the value of the mass is much bigger than the one of the other two parameters. We can say that we are conditioned by the size of the scalar field structure. If the impact parameter and/or the velocity are too big or too small the simulations won’t work properly either because the particle wouldn’t spend enough time close to the scalar field structure or because it would pass too far from it in such a way that the interaction wouldn’t produce measurable effects.. As can be seen in Fig. 3, the variation of the initial velocity $v_{i}$ and the impact parameter $y_{i}$ does not affect significantly the order of magnitude of the friction force, whereas one can see a systematic and clear variation of this value as one changes the value of the mass of the particle. Specifically, one verifies that

[TABLE]

where $\alpha_{r}$ depends on the initial velocity and the impact parameter, while $\alpha_{\phi}$ almost doesn’t depend on those parameters, presenting a value of the order $10^{-1}$ in all instances. The dependence of $\alpha_{r}$ in the $y_{i}$ and $v_{i}$ parameters is asymmetrical: the variation of the initial velocity doesn’t affect significantly its value – the order of magnitude does not change – whereas the bigger the impact parameter the smaller is the order of magnitude of the coefficient. In fact, while for cases in which $y_{i}=1,3$ we verify $\alpha_{r}\sim\mathcal{O}(1)$ for $y_{i}=5$ we observe $\alpha_{r}\sim\mathcal{O}(10^{-1})$. This comes as no surprise since the bigger the impact parameter the farther from the center of the scalar structure the particle will pass which means that the particle crosses regions where the scalar field is more and more diluted, decreasing the value of the friction force. Moreover, since the value we are calculating corresponds to a force per unit mass, which by the relations of Eq. (45) is proportional to the mass of the incoming particle, we conclude that the total force $\bm{F}$ scales as the square of the incoming particle $M_{p}$. This is not a new result (see, for instance Ostriker (1999)), but provides a connection between our study and the study of the drag force in self-interacting media.

The second comment has to do with the direction of the average force. We obtain that the average friction force has negative components in both planar directions.

IV.2.2 Loss of angular momentum

Another way to describe the effect of the interaction between the orbiting particle and the scalar field is to study the angular momentum of the former. We don’t use the energy because the energy depends on the gravitational potential which in our scenario is dynamical and so it doesn’t provide a good measure to characterize the movement of the orbiting particle. The angular momentum, however, is a good measure in the sense that it depends only on kinematic variables. Having said that, we are going to study the quantity

[TABLE]

where $\mathcal{L}_{\text{out}}$, $\mathcal{L}_{\text{in}}$ represent the angular momentum per unit mass of the orbiting particle when it leaves and when it enters the scalar field structure, respectively. We show in Fig. 4, the angular momentum per unit mass a function of time. From the picture, we can see that, as expected, the loss of angular momentum is mainly affected by the mass of the orbiting particle. Particularly, the following relation can be found

[TABLE]

with the proportionality factor, $\sigma$, varying slightly with the initial velocity and impact factor being, however, always of order $10^{-1}$.

IV.3 Changes in the density distribution of the field

The appearance of the friction force and the loss of angular momentum are related to the dynamical reaction of the scalar field to the presence of the incoming particle. In this respect, it is verified that the scalar field structure develops non-spherical over-densities that are time dependent. In order to appreciate this behavior, we will isolate the different components of the scalar field density. Using Eqs. (33) and (34) we can write the quantity $\rho_{f}=(4\pi)^{-1}\psi\psi^{*}$ (see Eq. (15)) as

[TABLE]

where

[TABLE]

To simplify the expression for $\rho_{f}$, we observe that (see footnote 6)

[TABLE]

which allows us to write, without loss of generality,

[TABLE]

and with that we can rewrite $\rho_{f}$ as

[TABLE]

where we fixed the value $\theta=\pi/2$ for the orbital plane.

IV.3.1 Time dependence of the non-spherical density

The time dependence of the functions $A(t,r)$ and $B(t,r)$ (see Eq. (53)) is illustrated in Fig. 5. In the figure, we see that the profile of the non-spherical components of the density evolves with time, a behavior that, combined with the angular dependence conveyed by the sinusoidal functions (see Eq. (54)), will result in rotating and oscillating non-spherical component of the density of the scalar field. A dramatic example of such behavior can be appreciated in Fig. 6 in which it is displayed the value of the density function in the plane $\theta=\pi/2$ for the simulation LY5V07. In this particular simulation, the incoming particle is scattered by the scalar field structure and moves past it. However, the short time in which the particle is close to the center of the scalar field structure is enough to give rise to a rotating over-density, as can be seen in the last contour plot presented in the respective figure. This example is simple to represent, however, when the particle stays in a bounded orbit, the non-spherical over-densities have a less organized behavior. In this case, there are two competing effects: on the one hand the scalar field structure is dictating the evolution of the non-spherical densities through the SP system and on the other hand the movement of the bounded particle introduces an oscillating “forcing term” on top of that.

In order to illustrate this behavior we will consider the position of the center of mass of the scalar field configuration, which, as shown below, is directly related to the functions $A$ and $B$ (see Eq. (53)). The position of the center of mass of the scalar field configuration $\bm{r}_{CM}=(x_{CM},y_{CM},z_{CM})$ is given by

[TABLE]

where $\rho_{f}(\bm{r})$ is given by Eq. (54). Then, the denominator is written as

[TABLE]

and since $W_{1}=\iint Y_{1-1}Y_{1,1}\sin\theta\mathrm{d}\theta\mathrm{d}\phi=-1$ we obtain that

[TABLE]

In all our simulations the movement is planar, so it suffices to calculate the $(x,y)$ coordinates of the center of mass (the origin of the coordinates is at the center of the initial configuration of the scalar field). We obtain that

[TABLE]

and

[TABLE]

We verify that the center of mass of the scalar configuration oscillates around its initial position – $x_{\mathrm{CM}}=y_{\mathrm{CM}}=0$ – and the magnitude of the oscillation depends mainly on the mass of the particle. As we can see in Figs. 7 and 8, bounded orbits will produce less organized oscillations of the coordinates of the center of mass whereas unbounded orbits create a more organized, regular pattern. Moreover, independently of the other parameters, one verifies that $\mathcal{O}(x_{\mathrm{CM}})\sim 10^{-1}$ for $M_{P}\sim 10^{-1}$, $\mathcal{O}(x_{\mathrm{CM}})\sim 10^{-2}$ for $M_{P}\sim 10^{-3}$ and $\mathcal{O}(x_{\mathrm{CM}})\sim 10^{-4}$ for $M_{P}\sim 10^{-5}$; the same relations hold for $y_{\mathrm{CM}}$.

As became clear in the study of the friction force and the loss of angular momentum, the mass of the particle is the most important factor determining the change in its dynamics. Taking that into account, we will focus on the characteristics of the movement of the center of mass of the scalar field configuration using the simulations SY5V07, MY5V07 and LY5V07. Using this set of simulations is appropriate for two reasons: 1) all the orbits are unbounded, which allows a clearer and simpler analysis of the dynamical aspects of the center of mass; 2) the simulations differ from each other by the value of the mass of the incoming particle, which is exactly the parameter that influences the most all the details of the dynamics of the system. We calculate the radial and angular velocity of the center of mass by using the values of $x_{\mathrm{CM}}$ and $y_{\mathrm{CM}}$ that we calculate directly from the simulation files (see Eqs. (58) and (59)). To do it, we use the following expressions

[TABLE]

in which $\dot{x}\equiv dx/dt$. We present the results of these calculations in Fig. 9 and from there two things are evident: 1) the behavior of the simulations with the smaller values of the mass of the particle are very similar, except for the frequency of oscillation and the magnitude; 2) the magnitude of the velocity components scales with the mass of the incoming particle, and we verify that

[TABLE]

IV.3.2 Magnitude of the non-spherical density

In Fig. 5 it is plotted the profile of the non-spherical components of the scalar field density, namely $A(t,r)$ and $B(t,r)$ for two different moments in time. The magnitude of these non-spherical components depends on the mass of the neighboring particle: the bigger the mass $M_{P}$, the bigger the magnitude of these componentes, as can be seen in Fig. 10. There, we represent the maximum value of the magnitude of functions $A$ and $B$ (see Eq. (53)) in the different sets of simulations MV1, MV2, MI1, MI2, VI1 and VI2. We see that, again, of the three initial parameters of the incoming particle, the mass has the strongest influence on the magnitude of the non-spherical components of the scalar field density. In fact, similarly to the case of the friction force, one can write

[TABLE]

i.e., the maximum magnitude of both $A$ and $B$ is directly proportional to the mass of the incoming particle, with the proportionality factors $\beta_{A,B}$ presenting values between 3 and 5 without any correlation with the initial velocity and impact parameter of the particle.

V Conclusion

We studied what happens to a stable, low-energy scalar field configuration when a particle-like body passes in its neighborhood. We observe that the mass of such body is the most important factor parametrizing the interaction between the two involved parties. Specifically, we verify that once the incoming particle dives in the scalar field structure, it is affected by a friction-like force that scales as $M_{p}^{2}$. This result bridges our calculations with other studies of friction force in self-interacting media. We also studied the effect of this friction force on the loss of angular momentum by the orbiting particle, with results that show the same scaling , i.e., the loss of angular momentum per unit mass scales with $M_{p}$. We verify the development of non-spherical components of the initially spherical scalar field structure due to the presence of the neighboring particle, their magnitude scaling with the particle’s mass $M_{p}$. Furthermore, we observe that even when the incoming body passes by the scalar field structure, i.e. describing an unbounded orbit, a non-spherical component develops and stays rotating around the center of the scalar strucure. The velocity of rotation of these non-spherical components scales with the mass of the orbiting particle.

The apparent strong connection between the mass of the incoming particle and the effects it leaves on the scalar field structure, may become a good lead in investigations about the history of such structures. In the scenarios in which DM is described by such scalar field structures, the detection of rotating clumps of DM in galaxies may be explained by a primordial encounter between the DM aggregate and a passing massive particle-like body. In such spirit, future work should include a more detailed analysis of the numerical approximations that were considered in the present work, as well as a simulation of a more rich scenario, in which the scalar field structure is orbited by stars as the encounter with the particle-like body happens.

Acknowledgements.

We thank Miguel Zilhão, Caio Macedo, David Hilditch and, specially, Vítor Cardoso for their insightful comments and feedback. We thankfully acknowledge the computer resources, technical expertise and assistance provided by Sérgio Almeida and Manuel Torrinha at CENTRA/IST. Computations were performed at the cluster “Baltasar-Sete-Sóis”, and supported by the MaGRaTh–646597 ERC Consolidator Grant. The author acknowledges financial support provided by Fundação para a Ciência e a Tecnologia Grant number PD/BD/113481/2015 awarded in the framework of the Doctoral Programme IDPASC - Portugal.

Appendix A Testing the code

We evolve the two components of our two-body system using different techniques. To solve SP system of equations, of the form of Eqs. (26) and (27), we use a centered finite difference stencil to write the derivatives. Particularly, at a generic point $u_{j}=j\Delta u$, we discretize the first derivatives as

[TABLE]

and the second derivatives as

[TABLE]

for a general function $H(u)$, indicating $H(u_{j})=H_{j}$. Having discretized the equations, we apply the iterated Crank Nicolson method with two iterations, following the conclusions of Ref. Teukolsky (2000). To solve the equations of motion of the point-particle, which can be cast in the generic form $dv/dt=G(t,v)$, we use Euler’s method, with the evolution step given by

[TABLE]

and the two-step Adams-Bashforth method given by

[TABLE]

A.1 Evolving a stationary scalar field solution

Using a timestep $\Delta t=10^{-3}$, we run a test with three different grid spacings - $\Delta r=0.2,0.1,0.05$. To quantify the effect of the grid spacing in the evolution of the field, define

[TABLE]

where $\rho_{E}(r)=f_{E}f_{E}^{*}$ is the equilibrium density of the scalar field (see Eq. (13)) and $\rho(t,r)$ is the density of the field that is evolved in time using our code. The results of this evolution are shown in Fig. 11.

The test allows us to conclude that with decreasing resolution, the magnitude of the deviations from the initial stationary configuration decreases. Moreover, we obtain that $\mathrm{max}\left[\Delta\rho\right]\sim\left(\Delta r\right)^{2}$.

A.2 Testing the code evolving the orbiting particle

The evolution of the particle will be made with the same time step as the one used for the SP-equations. To correctly describe this evolution, the code has to guarantee the conservation of the energy and angular momentum per unit mass for the control tests. To visualize that conservation, we calculate the following quantities

[TABLE]

where $\varepsilon(0)$ and $\mathcal{L}(0)$ represent the initial energy and angular momentum per unit mass, respectively.

In Fig. 12 we show the evolution of these quantities for the control test Y3V03. We observe that both the energy and the angular momentum are conserved, in the worst case, up to the percent level.

Appendix B Discretizing the Dirac delta

In order to describe the perturbing mass orbiting the scalar configuration as a point particle, it is necessary to use the Dirac delta. To describe it in a numerical grid, we follow an approach used in previous works (see Ref. Sundararajan et al. (2007) and references therein) in which the construction of the discretized version of the Dirac delta is made considering that its defining feature is integrability. This means that we obtain the discretized Dirac delta by studying the expression

[TABLE]

for some “well-behaved” (continuous with continuous derivatives) function $f(x)$. The finite difference version of the previous expression is given by

[TABLE]

where $f_{i}$ and $\delta_{i}$ represent the values of function $f$ and the Dirac delta, respectively, at the grid point $i$. In our case it suffices to consider that the point $X_{0}$ is always a grid point such that $X_{0}/dx=i*$. With this setup, we can say that the only point of the grid in which the Dirac delta takes a non-zero value is precisely the grid point corresponding to $X_{0}$. This implies that the finite difference formula of Eq. (73) can be written as

[TABLE]

from where we can read that the Dirac delta has the following finite difference representation

[TABLE]

This definition agrees, in the respective limits, with the simplest definition for the Dirac delta of Sundararajan et al. (2007). We decided to use a one-point-only discretized Dirac delta for two reasons: a) it works well; b) given the scales of the problem at hand, we want to reinforce as much as possible the localized nature of the perturbing mass.

Appendix C Derivation of the Schrodinger-Poisson system with orbiting particle

In this appendix, we reintroduce the constants $G$, $\hbar$ and $c$ and we will be using the following post-Newtonian metric

[TABLE]

C.1 Poisson equation

The Einstein tensor, with the metric of Eq. (76), is given by

[TABLE]

The energy momentum tensor of the scalar field and the point-particle is written as

[TABLE]

where

[TABLE]

corresponds to the energy-momentum tensor of the scalar field, and888This result is obtained after integrating in time the expression of Eq. (19).

[TABLE]

where $g$ is the determinant of the metric, is the energy-momentum tensor due to a point particle. Using these definitions, one can write

[TABLE]

where we separate the term of order $c^{-4}$ from the others because this will be important to define the low-energy regime and where $\partial_{0}\equiv\partial/\partial(ct)$ with $t$ indicating the time coordinate. For the particle, the energy momentum is given by

[TABLE]

Considering that the velocity of the particle is much smaller than the velocity of light, its energy momentum tensor reduces to

[TABLE]

which is a reflection of the fact that the behavior of the particle is dominated by its rest mass. For the scalar field, the low-energy regime also means that its behavior is dominated by its rest mass. Since a scalar field of mass $m_{s}$ and momentum $p$ oscillates with a frequency $\omega=E/\hbar$, i.e.

[TABLE]

where $E=\sqrt{p^{2}c^{2}+m_{s}^{2}c^{4}}$ is the energy of the field999This expression is valid only for free fields, but the rest of the calculations will be valid as long as, up to first order, $E\sim m_{s}c^{2}$, which means that the calculation is also valid for bound states of the scalar field whose energy is dominated by the rest-mass term., in a low energy regime we rightly consider that

[TABLE]

where $\psi(t,\bm{x})$ constains not only the spatial dependence of the field but also a slowly varying dependence in time. Because we know that in the low-energy limit the rest-mass dominates, we apply the condition $p^{2}\ll m_{s}c^{2}$ such that

[TABLE]

which allows one to infer that

[TABLE]

This form of the scalar field reprensents the low-energy limit and the derivative terms in the energy-momentum tensor can thus be simplified as (remember that $\partial_{0}\equiv\partial/\partial(ct)$)

[TABLE]

Plugging this result in the expression for the energy-momentum tensor, we obtain

[TABLE]

Notice that $\partial_{t}\psi\sim p^{2}/2m_{s}$ so that

[TABLE]

because, again, we are in the low energy limit. Finally, the energy momentum tensor of the scalar field in the low energy limit is given by

[TABLE]

Having defined the low-energy limit, we can now equate both sides of the Einstein equations

[TABLE]

which reduce to its 00-component in the low-energy regime:

[TABLE]

C.2 Klein-Gordon equation

Plugging in the Klein-Gordon equation

[TABLE]

the post-Newtonian metric of Eq. (76) and the scalar field

[TABLE]

we obtain

[TABLE]

From here, we see that the dynamics described by the Klein-Gordon equation is dominated by the Schrodinger equation in the low-energy limit.