* [127] A. V. Chechkin, I. M. Zaid, M. Lomholt, I. M. Sokolov and R. Metzler, Bulk-mediated surface diffusion along a cylinder: Propagators and crossovers, Phys. Rev. E 79, 040105(R) (2009).
  * [128] O. Bénichou, D. S. Grebenkov, P. E. Levitz, C. Loverdo and R. Voituriez, Optimal Reaction Time for Surface-Mediated Diffusion, Phys. Rev. Lett. 105, 150606 (2010).
  * [129] O. Bénichou, D. S. Grebenkov, P. E. Levitz, C. Loverdo and R. Voituriez, Mean first-passage time of surface-mediated diffusion in spherical domains, J. Stat. Phys. 142, 657 (2011).
  * [130] F. Rojo and C. E. Budde, Enhanced diffusion through surface excursion: A master-equation approach to the narrow-escape-time problem, Phys. Rev. E 84, 021117 (2011).
  * [131] A. V. Chechkin, I. M. Zaid, M. A. Lomholt, I. M. Sokolov, and R. Metzler, Effective surface motion on a reactive cylinder of particles that perform intermittent bulk diffusion, J. Chem. Phys. 134, 204116 (2011).
  * [132] A. V. Chechkin, I. M. Zaid, M. Lomholt, I. M. Sokolov, and R. Metzler, Bulk-mediated diffusion on a planar surface: Full solution, Phys. Rev. E 86, 041101 (2012).
  * [133] J.-F. Rupprecht, O. Bénichou, D. S. Grebenkov and R. Voituriez, Exact mean exit time for surface-mediated diffusion, Phys. Rev. E 86, 041135 (2012).
  * [134] J.-F. Rupprecht, O. Bénichou, D. S. Grebenkov and R. Voituriez, Kinetics of active surface-mediated diffusion in spherically symmetric domains, J. Stat. Phys. 147, 891 (2012).
  * [135] A. M. Berezhkovskii, L. Dagdug, and S. M. Bezrukov, A new approach to the problem of bulk-mediated surface diffusion, J. Chem. Phys. 143, 084103 (2015).
  * [136] A. M. Berezhkovskii, L. Dagdug, and S. M. Bezrukov, Bulk-mediated surface transport in the presence of bias, J. Chem. Phys. 147, 014103 (2017).
  * [137] P. C. Bressloff, Renewal equation for single-particle diffusion through a semipermeable interface, Phys. Rev. E. 107, 014110 (2023).
  * [138] P. C. Bressloff, Renewal equations for single-particle diffusion in multi-layered media (Preprint, arXiv:2301.02895)
  * [139] D. S. Grebenkov, Probability distribution of the boundary local time of reflected Brownian motion in Euclidean domains, Phys. Rev. E 100, 062110 (2019).
  * [140] P. Levitz, M. Zinsmeister, P. Davidson, D. Constantin, and O. Poncelet, Intermittent Brownian dynamics over a rigid strand: Heavily tailed relocation statistics in a simple geometry, Phys. Rev. E 78, 030102(R) (2008).
  * [141] A. J. Bray, S. N. Majumdar and G. Schehr, Persistence and first-passage properties in non-equilibrium systems, Adv. Phys. 62, 225 (2013).
  * [142] N. Levernier, M. Dolgushev, O. Bénichou, R. Voituriez, and T. Guérin, Survival probability of stochastic processes beyond persistence exponents, Nat. Comm. 10, 2990 (2019).
  * [143] D. S. Grebenkov, Statistics of boundary encounters by a particle diffusing outside a compact planar domain, J. Phys. A.: Math. Theor. 54, 015003 (2021).
  * [144] A. Talbot, The accurate numerical inversion of Laplace transforms, IMA J. Appl. Math. 23, 97-120 (1979).
  * [145] L. L. Latour, K. Svoboda, P. P. Mitra, and C. H. Sotak, Time-dependent diffusion of water in a biological model system, Proc. Natl. Acad. Sci. USA 91, 1229-1233 (1994).
  * [146] P. N. Sen, Time-dependent diffusion coefficient as a probe of geometry, Conc. Magn. Reson. 23A, 1-21 (2004).
  * [147] J. Wu and K. M. Berland, Propagators and Time-Dependent Diffusion Coefficients for Anomalous Diffusion? Biophys. J. 95, 2049-2052 (2008).
  * [148] E. W. Montroll and G. H. Weiss, Random Walks on Lattices. II, J. Math. Phys. 6, 167 (1965).
  * [149] R. Metzler and J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach, Phys. Rep. 339, 1-77 (2000).
  * [150] J. Klafter and I. M. Sokolov, First Steps in Random Walks: From Tools to Applications (Oxford University Press, 2011).
  * [151] S. B. Yuste, E. Abad, and K. Lindenberg, Exploration and trapping of mortal random walkers, Phys. Rev. Lett. 110, 220603 (2013).
  * [152] B. Meerson and S. Redner, Mortality, Redundancy, and Diversity in Stochastic Search, Phys. Rev. Lett. 114, 198101 (2015).
  * [153] D. S. Grebenkov and J.-F. Rupprecht, The escape problem for mortal walkers, J. Chem. Phys. 146, 084106 (2017).
  * [154] B. Meerson, Mortal Brownian motion: Three short stories, Int. J. Mod. Phys. B 33, 1950172 (2019).
  * [155] M. V. Chubynsky and G. W. Slater, Diffusing diffusivity: a model for anomalous, yet Brownian, diffusion, Phys. Rev. Lett. 113, 098302 (2014).
  * [156] A. V. Chechkin, F. Seno, R Metzler, and I. M. Sokolov, Brownian yet non-Gaussian diffusion: from superstatistics to subordination of diffusing diffusivities, Phys. Rev. X 7, 021002 (2017).
  * [157] D. S. Grebenkov, An encounter-based approach for restricted diffusion with a gradient drift, J. Phys. A: Math. Theor. 55, 045203 (2022).
  * [158] M. R. Evans, S. N. Majumdar, and G. Schehr, Stochastic resetting and applications J. Phys. A: Math. Theor. 53, 193001 (2020).
  * [159] P. C. Bressloff, Diffusion-mediated surface reactions and stochastic resetting, J. Phys. A: Math. Theor. 55, 275002 (2022).
  * [160] Z. Benkhadaj and D. S. Grebenkov, Encounter-based approach to diffusion with resetting, Phys. Rev. E 106, 044121 (2022).
  * [161] D. S. Grebenkov, Encounter-based approach to the escape problem, Phys. Rev. E 107, 044105 (2023).
  * [162] D. S. Grebenkov, Depletion of Resources by a Population of Diffusing Species, Phys. Rev. E 105, 054402 (2022).
  * [163] P. C. Bressloff, Narrow capture problem: an encounter-based approach to partially reactive targets, Phys. Rev. E 105, 034141 (2022).
  * [164] D. S. Grebenkov, Statistics of diffusive encounters with a small target: Three complementary approaches, J. Stat. Mech. 083205 (2022).
  * [165] A. Szabo, Theory of diffusion-influenced fluorescence quenching, J. Phys. Chem. 93, 6929 (1989).
  * [166] D. S. Grebenkov, M. Filoche, and B. Sapoval, Spectral Properties of the Brownian Self-Transport Operator, Eur. Phys. J. B 36, 221-231 (2003).
  * [167] S. J. Chapman, R. Erban, and S. A. Isaacson, Reactive boundary conditions as limits of interaction potentials for Brownian and Langevin dynamics, SIAM J. Appl. Math 76, 368-390 (2016).
  * [168] D. S. Grebenkov, Imperfect Diffusion-Controlled Reactions, in ”Chemical Kinetics: Beyond the Textbook”, Eds. K. Lindenberg, R. Metzler, G. Oshanin (World Scientific, New Jersey, 2019), pp. 191-219.
  * [169] F. Piazza, The physics of boundary conditions in reaction-diffusion problems, J. Chem. Phys. 157, 234110 (2022).


# Uncovering a hidden black hole binary from secular eccentricity variations
of a tertiary star

Bin Liu<EMAIL_ADDRESS>Niels Bohr International Academy, Niels Bohr
Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark  Daniel J. D’Orazio  Niels
Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100
Copenhagen, Denmark  Alejandro Vigna-Gómez  Niels Bohr International Academy,
Niels Bohr Institute, Blegdamsvej 17, 2100 Copenhagen, Denmark  Johan Samsing
Niels Bohr International Academy, Niels Bohr Institute, Blegdamsvej 17, 2100
Copenhagen, Denmark

###### Abstract

We study the dynamics of a solar-type star orbiting around a black hole binary
(BHB) in a nearly coplanar system. We present a novel effect that can prompt a
growth and significant oscillations of the eccentricity of the stellar orbit
when the system encounters an “apsidal precession resonance”, where the
apsidal precession rate of the outer stellar orbit matches that of the inner
BHB. The eccentricity excitation requires the inner binary to have a non-zero
eccentricity and unequal masses, and can be created even in non-coplanar
triples. We show that the secular variability of the stellar orbit’s
apocenter, induced by the changing eccentricity, could be potentially
detectable by Gaia. Detection is favorable for BHBs emitting gravitational
waves in the frequency band of the Laser Interferometer Space Antenna (LISA),
hence providing a distinctive, multi-messenger probe on the existence of
stellar-mass BHBs in the Milky Way.

##  I Introduction

About 90 double compact object merger events have been detected by the LIGO-
Virgo-KAGRA Collaboration in the first three observing runs LIGO-2021 . The
avenues to produce stellar-mass black-hole binary (BHB) mergers include
different formation channels and environments. Some of them are isolated
binary evolution Lipunov-1997 ; Lipunov-2007 ; Podsiadlowski-2003 ;
Belczynski-2010 ; Belczynski-2016 ; Dominik-2012 ; Dominik-2013 ; Dominik-2015
; Alejandro-2017 , chemically homogeneous evolution Mandel-2016 ;
Marchant-2016 ; duBuisson ; Riley , and multiple-body evolution in the gas
disks of active galactic nuclei Baruteau-2011 ; McKernan-2012 ; McKernan-2018
; Bartos-2017 ; Stone-2017 ; Leigh-2018 ; Secunda-2019 ; Yang-2019 ;
Grobner-2020 ; Ishibashi-2020 ; Tagawa-2020 ; Liyaping-2021 ; Ford-2021 ;
Samsing-Nature ; Lirixin-2022 ; Lijiaru-2022 . Additionally,there are various
flavors of dynamical channels that involve either strong gravitational
scatterings in dense clusters Zwart(2000) ; OLeary(2006) ; Miller(2009) ;
Banerjee(2010) ; Downing(2010) ; Ziosi(2014) ; Rodriguez(2015) ; Samsing(2017)
; Samsing(2018) ; Rodriguez(2018) ; Gondan(2018) , tertiary-induced mergers
via von Zeipel-Lidov-Kozai (ZLK) oscillations vonZeipel ; Lidov ; Kozai ;
Smadar ; Miller-2002 ; Wen-2003 ; Antonini-2012 ; Antonini(2017) ;
Silsbee(2017) ; Petrovich-2017 ; Liu-ApJ ; Xianyu-2018 ; Hoang-2018 ; Liu-
Quadruple ; Fragione-Quadruple ; Fragione-nulearcluster ; Zevin-2019 ; Liu-
HierarchicalMerger  or flyby-induced mergers Michaely-2019 ; Michaely-2020 .
However, the relative contribution of each channel and the astrophysical
origin of the detected mergers is still unclear.

BHB progenitors are expected to be numerous but remain undetected as an
abundant population in our Universe. Searching for these inspiraling BHBs is
of great importance to understand the origin of gravitational-wave (GW)
sources. If the BHBs are not accreting, these quiescent sources can be
detected via GWs. Since BHBs in the inspiral phase are still far from merger,
the associated GWs are in the low-frequency band that can be explored by
future spaceborne GW observatories, such as LISA LISA , TianQin TianQin ,
Taiji TaiJi , B-DECIGO DECIGO , Decihertz Observatories DeciHZ , and TianGO
TianGo . Also, since a significant fraction of compact BHBs may be members of
hierarchical systems Tokovinin ; Raghavan ; Fuhrmann , the motion of a nearby
visible object (such as a star or a pulsar) can be used to search for BHBs. In
this scenario, the inner BHB can perturb the outer orbit, either inducing
short-term orbital oscillations (tertiary orbit becomes quasi-Keplerian), or
causing long-term oscillations of the eccentricity and the orientation of the
angular momentum when the tertiary orbit is highly inclined Naoz-2017 ;
Chiang-2018 ; Naoz-2019 . If the tertiary object is bright enough, radial
velocity measurements could be used to determine the short and long term
deviations from a Keplerian orbit Suto-1 ; Suto-2 ; Suto-3 .

Currently, observations show that most triple-star systems are less inclined
or nearly coplanar Ransom ; Thompson ; Eisner . Here, we consider a solar-type
star orbiting around a BHB and study the secular evolution of the stellar
orbit when the triple system is coplanar. We show that the outer binary may
experience an eccentricity growth driven by the “apsidal precession resonance”
Liu-Yuan ; Liu-2020-PRD . Compared to the complex evolution of the orientation
of the angular momentum (i.e., precession or nutation), the secular change of
the eccentricity of the outer stellar orbit could provide distinctive evidence
to reveal the presence of a BHB.

##  II Apsidal precession resonance

We consider an inner BHB with masses $m_{1}$, $m_{2}$, and a solar-type star
($m_{\star}=1\ \rm{M_{\odot}}$) that moves around the center of mass of the
inner bodies. The reduced mass for the inner binary is
$\mu_{\mathrm{in}}\equiv m_{1}m_{2}/m_{12}$, with $m_{12}\equiv m_{1}+m_{2}$.
Similarly, the outer binary has
$\mu_{\star}\equiv(m_{12}m_{\star})/(m_{12}+m_{\star})$. The semi-major axes
and eccentricities are denoted by $a_{\mathrm{in}}$, $a_{\star}$ and
$e_{\mathrm{in}}$, $e_{\star}$, respectively. The orbital angular momenta of
two orbits are thus given by
$\textbf{L}_{\mathrm{in}}=\mathrm{L}_{\mathrm{in}}\hat{\textbf{L}}_{\mathrm{in}}=\mu_{\mathrm{in}}\sqrt{Gm_{12}a_{\mathrm{in}}(1-e_{\mathrm{in}}^{2})}\,\hat{\textbf{L}}_{\mathrm{in}}$
and
$\textbf{L}_{\star}=\mathrm{L}_{\star}\hat{\textbf{L}}_{\star}=\mu_{\star}\sqrt{G(m_{12}+m_{\star})a_{\star}(1-e_{\star}^{2})}\,\hat{\textbf{L}}_{\star}$.

When the triple system is less inclined or nearly coplanar, the ZLK
oscillations are not allowed to occur, but a significant eccentricity
excitation of the inner binary may still be induced Ford ; Naoz-PN . A
secular, “apsidal precession resonance” plays a dominant role if the total
apsidal precession of the inner binary matches the precession rate of the
outer binary Liu-Yuan ; Liu-2020-PRD . Precession of both inner and outer
binaries is driven by Newtonian and general relativistic (GR) effects . Such
resonance allows efficient angular momentum exchange between the inner and
outer binaries.

Figure 1: Apsidal precession resonance in coplanar triple systems where the
outer orbital period ($P_{\star}$) is set to be 15 days. All the results are
obtained by integrating the SA secular equations including GR effects (but
without GW emission). Left panel: Evolution examples of the orbital
eccentricities in the inner (black) and outer (red) binaries, with
$m_{12}=50\rm{M}_{\odot}$ and initial $e_{\mathrm{in}}^{0}=0.9$ and
$e_{\star}^{0}=0$. Middle panel: The ratio of apsidal precession rates
(Equations 3-4; including the dependence of finite eccentricities, i.e.,
$e_{\mathrm{in}}^{0}$, $e_{\star}>0$) and the maximum change of eccentricity
$\Delta e_{\star}^{\mathrm{max}}$ of the outer stellar orbit as a function of
the semimajor axis $a_{\mathrm{in}}$. The cross-hatched region corresponds to
dynamically unstable triple systems Kiseleva-1996 . The parameters are the
same as in the left panel, except $a_{\mathrm{in}}$ is relaxed to a range of
values. The labeled GW frequency is the peak frequency at pericenter Wen-2003
. The numerical results are obtained by integrating the SA secular equations
over several timescales (as labeled), and the analytical result is given by
the energy and angular momentum conservation laws (dashed blue line). Right
panel: $\Delta e_{\star}^{\mathrm{max}}$ induced by the apsidal precession
resonance in the ($m_{2}/m_{1}-a_{\mathrm{in}}$) plane. We use the same
parameters as in the middle panel, taking into account the full range of mass
ratios of the BHB. The three black contours (solid, dashed and dot-dashed)
specify $\Delta e_{\star}^{\mathrm{max}}=0.1,0.2$ and $0.5$, respectively. The
gray crosses indicate the significant change of $e_{\star}$ that leads to
instability.

Here, we extend our previous studies to an “inverse” secular problem, and
address the question of how the apsidal precession resonance modifies the
eccentricity evolution of the tertiary for the first time. Since we are
interested in the long-term orbital evolution, we adopt the single-averaged
(SA; only averaging over the inner orbital period) secular equations of
motion, taking into account the contributions from the Newtonian effect up to
the octupole level of approximation and the leading order GR effect in both
inner and outer orbits. The explicit SA equations are provided in Liu-ApJ ;
Blanchet .

In Figure 1, the left panels show the evolution of eccentricities of both
inner and outer orbits. Starting with the eccentric inner binary, we see that
$e_{\star}$ can be excited from the circular orbit and undergoes oscillations.
In the middle panel, we find that a large fraction of systems can develop
eccentricities ($0.04\lesssim a_{\mathrm{in}}/\mathrm{AU}\lesssim 0.08$), and
an evident peak, $\Delta e_{\star}^{\mathrm{max}}\simeq 0.8$, can be resolved
when the evolution is sufficiently long ($\gtrsim 10$yrs). The right panel
illustrates the level of $e_{\star}-$excitation in the
($m_{2}/m_{1}-a_{\mathrm{in}}$) plane. The eccentricity of the stellar orbits
can be excited for $a_{\mathrm{in}}\gtrsim 0.04\ \mathrm{AU}$, and the systems
with smaller mass ratio tend to have larger $\Delta e_{\star}^{\mathrm{max}}$.
This is because the evolution of $e_{\star}$ is determined by the octupole-
order secular interactions, which can be quantified by terms proportional to
Smadar :

$\varepsilon_{\mathrm{oct}}\equiv\frac{m_{1}-m_{2}}{m_{12}}\frac{a_{\mathrm{in}}}{a_{\star}}\frac{e_{\star}}{1-e_{\star}^{2}}.$
(1)

We see that the eccentricities of some outer orbits reach significantly large
$\Delta e_{\star}^{\mathrm{max}}$ close to unity, leading to unbound orbits.

Figure 2:  The maximum change of eccentricity as a function of the inner
binary semimajor axis. We consider the inner BHB with the total mass
$m_{12}=50\ \rm{M}_{\odot}$ and mass ratio $m_{2}/m_{1}=0.2$ (same as the
middle panel of Figure 1), and the stellar orbits with different orbital
period (as labeled). The results (solid lines) are obtained by the numerical
integrations of SA equations for 25 yrs (left panel), 150 yrs (middle panel)
and 300 yrs (right panel), and the results (dashed lines) are all from 10 yrs.
To compare with the fiducial example (black lines), we fix the initial
$e_{\star}^{0}$ ($e_{\mathrm{in}}^{0}$) in the upper (lower) panels (as
labeled).

For coplanar ($\hat{\mathbf{L}}=\hat{\mathbf{L}}_{\star}$), non-dissipative
(no gravitational radiation) systems, the secular dynamics can be understood
analytically. When $e_{\mathrm{in}}$, $e_{\star}\ll 1$, the evolution of
$\mathbf{e}_{\mathrm{in}}$ and $\mathbf{e}_{\star}$ is governed by the linear
Laplace-Lagrange equations MD ; YanqingWu . If we define the complex
eccentricity variables as $\mathcal{E}_{\mathrm{in}}\equiv
e_{\mathrm{in}}\mathrm{exp}(i\varpi_{\mathrm{in}})$ and
$\mathcal{E}_{\star}\equiv e_{\star}\mathrm{exp}(i\varpi_{\star})$, where
$\varpi_{\mathrm{in}}$, $\varpi_{\star}$ are the longitude of pericenter of
the inner and outer orbits, then the evolution equations are reduced to

$\displaystyle\frac{d}{dt}\begin{pmatrix}\mathcal{E}_{\mathrm{in}}\\\
\mathcal{E}_{\star}\end{pmatrix}$
$\displaystyle=i\begin{pmatrix}\omega_{\mathrm{in}}&\nu_{\mathrm{in},\star}\\\
\nu_{\star,\mathrm{in}}&\omega_{\star}\end{pmatrix}\begin{pmatrix}\mathcal{E}_{\mathrm{in}}\\\
\mathcal{E}_{\star}\end{pmatrix},$ (2)

with

$\displaystyle\omega_{\mathrm{in}}=\frac{3}{4}n_{\mathrm{in}}\frac{m_{\star}}{m_{12}}\bigg{(}\frac{a_{\mathrm{in}}}{a_{\star}}\bigg{)}^{3}+\omega_{\mathrm{GR},\mathrm{in}},$
(3)
$\displaystyle\omega_{\star}=\frac{3}{4}n_{\star}\frac{m_{1}m_{2}}{m_{12}^{2}}\bigg{(}\frac{a_{\mathrm{in}}}{a_{\star}}\bigg{)}^{2}+\omega_{\mathrm{GR},\star},$
(4)
$\displaystyle\nu_{\mathrm{in},\star}=-\frac{15}{16}n_{\mathrm{in}}\bigg{(}\frac{a_{\mathrm{in}}}{a_{\star}}\bigg{)}^{4}\frac{m_{\star}(m_{1}-m_{2})}{m_{12}^{2}},$
(5)
$\displaystyle\nu_{\star,\mathrm{in}}=-\frac{15}{16}n_{\star}\bigg{(}\frac{a_{\mathrm{in}}}{a_{\star}}\bigg{)}^{3}\frac{m_{1}m_{2}(m_{1}-m_{2})}{m_{12}^{3}},$
(6)

where $n_{\mathrm{in}(\star)}=(Gm_{12}/a_{\mathrm{in}(\star)}^{3})^{1/2}$ and
$\omega_{\mathrm{GR},{\mathrm{in}(\star)}}=3G^{3/2}m_{12}^{3/2}/[c^{2}a_{\mathrm{in}(\star)}^{5/2}(1-e_{\mathrm{in}(\star)}^{2})]$
are the mean motion and the GR-induced pericenter-precession frequency of the
inner (outer) binary for $m_{12}\gg m_{\star}$, respectively.

Starting with $e_{\mathrm{in}}=e_{\mathrm{in}}^{0}$, $e_{\star}=0$ at $t=0$,
Equation (2) can be solved to determine the time evolution of $e_{\star}(t)$
(see Liu-Yuan, ), which oscillates between $0$ and $e_{\star}^{\mathrm{max}}$,
where

$e_{\star}^{\mathrm{max}}=2e_{\mathrm{in}}^{0}\frac{|\nu_{\star,\mathrm{in}}|}{\sqrt{(\omega_{\mathrm{in}}-\omega_{\star})^{2}+4\nu_{\mathrm{in},\star}\nu_{\star,\mathrm{in}}}}.$
(7)

Clearly, $e_{\star}^{\mathrm{max}}$ attains its peak value when
$\omega_{\mathrm{in}}=\omega_{\star}$ occurs, and then we have

$\displaystyle
e^{\mathrm{peak}}_{\star}=e_{\mathrm{in}}^{0}\frac{\mu^{1/2}_{\mathrm{in}}}{\mu^{1/4}_{\star}m^{1/4}_{\star}}\bigg{(}\frac{a_{\mathrm{in}}}{a_{\star}}\bigg{)}^{1/4}.$
(8)

Figure 3:  Similar to Figure 2, but we set the initial eccentricities as
$e_{\mathrm{in}}^{0}=0.9$ and $e_{\star}^{0}=0$, and consider different
initial inclinations (as labeled; $I_{\star}$ is the inclination angle between
$\textbf{L}_{\mathrm{in}}$ and $\textbf{L}_{\star}$). The changes of
$e_{\star}$ and $I_{\star}$ are shown in the upper and lower panels,
respectively.

Note that the linear theory is valid for the low$-e$ systems. For the example
in Figure 1, non-zero $e_{\mathrm{in}}^{0}$ may lead to unphysically large
$e^{\mathrm{peak}}_{\star}$; however, Equations (7) and (8) are useful in the
sense that we can expect: i) $e_{\star}-$excitation appears when
$\omega_{\mathrm{in}}\simeq\omega_{\star}$ (see the middle panel of Figure 1;
the resonance occurs with finite $e_{\star}$); ii) $e_{\star}-$excitation
becomes stronger with increasing $e_{\mathrm{in}}$ (see also Figure 2 below).

For finite eccentricities, Equation (2) breaks down. However, in the case of
exact coplanarity, the maximum eccentricity of a triple can still be
calculated algebraically, using energy and angular momentum conservations Lee
. This method works well for two orbits with arbitrary eccentricities, but it
cannot show the time evolution and resonance features. A full derivation can
be found in Liu-2020-PRD  and the solution is presented in Figure 1.

Figure 2 shows $e_{\star}-$excitation with the arbitrary initial
eccentricities for coplanar triples. We again consider the fiducial example,
in which the inner BHB has the total mass $m_{12}=50\ \rm{M}_{\odot}$ with
mass ratio $m_{2}/m_{1}=0.2$. We choose three values of the semimajor axis of
the outer orbit ($a_{\star}$) and consider a range of $a_{\mathrm{in}}$ that
satisfies the stability criterion. To illustrate the role of the eccentricity,
we start with the same initial configurations, namely, the argument of
periapse, the longitude of ascending nodes and the true anomaly of the outer
orbit are set to be the same at $t=0$. Each system is evolved for a long
timescale (to achieve the highest value of $e_{\star}$) and a short timescale
(10 yrs), respectively. The maximum change of $e_{\star}$ is picked only for
the system remains stable.

In the upper left panel of Figure 2, we see that all the initial circular
outer orbits can become eccentric when the resonance occurs, and the maximum
change $\Delta e_{\star}^{\mathrm{max}}$ grows as $e_{\mathrm{in}}^{0}$
increases. In the lower left panel, $e_{\star}-$excitation can still occur for
the initially eccentric outer orbits. However, due to the stability, only a
fraction of systems (with a narrow range of $a_{\mathrm{in}}$; when
$e_{\star}^{0}\gtrsim 0.6$) may undergo significant eccentricity oscillations.
Note that for the wider stellar orbits ($P_{\star}=50,100$ days), as shown in
the middle and right panels, the system has to be evolved for a sufficiently
long time. This is because the timescale of $e_{\star}-$excitation is of the
order of Liuetal-2015

$T_{e_{\star}}\big{|}_{e_{\mathrm{in}},e_{\star}\ll
1}\simeq\frac{m_{12}}{n_{\mathrm{in}}\varepsilon_{\mathrm{oct}}}\bigg{(}\frac{a_{\star}}{a_{\mathrm{in}}}\bigg{)}^{3}\frac{L_{\star}}{L_{\mathrm{in}}}.$
(9)

When the inner and outer orbits are mutually inclined, no simple analytical
result can be derived, and the long-term evolution of the outer orbits can
only be studied numerically. Previous studies Naoz-2017 ; Chiang-2018 ;
Naoz-2019  showed that in inclined triple systems, the eccentricity of the
outer orbit can oscillate moderately and the angular momentum
($\textbf{L}_{\star}$) undergoes nodal precession/nutation around the inner
one ($\textbf{L}_{\mathrm{in}}$). Moreover, $\textbf{L}_{\star}$ might
experience a flip if the tertiary is a test particle.

Figure 3 presents the results of the triples with a series of initial
inclinations. We see that in the upper panels, regardless of the values of the
initial inclinations, $e_{\star}-$excitation due to resonance always occurs,
and the resonance location shifts when $I_{\star}^{0}$ changes
($I_{\star}^{0}$ is the initial inclination angle between
$\textbf{L}_{\mathrm{in}}$ and $\textbf{L}_{\star}$). In the lower panels, we
find that the inclination varies for a wider range of $a_{\mathrm{in}}$
compared to the change of $e_{\star}$. In particular, $\Delta I_{\star}$
always undergoes an additional excitation when $\Delta
e_{\star}^{\mathrm{max}}$ approaches the peak value for the inclined systems.

##  III Resonance in stellar orbits and detectability

We now focus on dependence of the change in eccentricity on parameters of the
outer stellar orbit, considering the inner BHB is a LISA source. To explore
the observability of this effect, we consider the detectability of a maximum
eccentricity change, $\Delta e_{\star}^{\mathrm{peak}}$, within a certain
timescale.

We initialize systems with $e_{\mathrm{in}}^{0}=0.9$ and $e_{\star}^{0}=0$ in
a coplanar configuration. For the inner BHB, we choose three values of the
total mass ($m_{12}/\rm{M_{\odot}}=20,50,100$), and allow the mass-ratio range
to take on all values such that both masses are consistent with being a BH,
$m_{1},m_{2}\gtrsim 5\rm{M}_{\odot}$. To develop the joint detection with GW
detectors in the coming future, we focus on the BHBs radiating GWs in the LISA
frequency band. Thus, the semimajor axis $a_{\mathrm{in}}$ is chosen from a
uniform distribution that satisfies $f_{\mathrm{GW}}\geq 10^{-4}$ Hz and with
merger time (due to GW emission) larger than $10^{3}$ yrs Peters-1964 . Then,
for the (outer) stellar orbits, we sample the semimajor axis $a_{\star}$ from
a range of $P_{\star}\sim[1,180]$ days, considering only systems that are
dynamically stable. Each system evolves to 10, 30 or 100 yrs, using the SA
equations of motion. The maximum change of the $e_{\star}$ is recorded if the
star remains gravitationally bound and stable during the evolution. Finally,
since the orbital evolution relies on the initial geometry of the triples, to
cover all possibilities, we randomly sample the longitude of pericenter of the
inner orbit and the true anomaly of the outer orbit. Each triple system is
evolved with 100 different initial geometries.

Figure 4: Top panel: Parameter space in [$a_{\mathrm{in}}-a_{\star}$ or
$(P_{\star})$] plane, where the apsidal precession resonance occurs. We
consider the BHB with $m_{12}=50\ \rm{M}_{\odot}$, $m_{2}/m_{1}=0.2$ and the
initial $e_{\mathrm{in}}^{0}=0.9$, $e_{\star}^{0}=0$. The blue region is given
by Equations 3-4 including the finite $e_{\mathrm{in}}^{0}$ and $e_{\star}$
(as labeled). The dashed lines characterize the frequency of GWs emitted by
the inner BHB. Lower panels: Maximum of $\Delta e_{\star}^{\mathrm{max}}$
(i.e., $\Delta e_{\star}^{\mathrm{peak}}$) as a function of $a_{\star}$ for
different $m_{12}$. The solid and dashed lines are obtained by the numerical
integrations for 10, 30 and 100 yrs, respectively. The vertical blue line
refers to $a_{\star}=a_{\star}^{\mathrm{Res}}$ at $e_{\star}=0$ via
$f_{\mathrm{GW}}=10^{-4}$ Hz. The cross-hatched region in the bottom panel
indicates that the systems eventually become unbound when evolving for a
longer timescale.

For a given set of parameters, the criterion of apsidal precession resonance
($\omega_{\mathrm{in}}=\omega_{\star}$) provides a good estimate for the
resonance radius: resonance occurs at the location
$a_{\star}=a_{\star}^{\mathrm{Res}}$ for a given $a_{\mathrm{in}}$. In Figure
4, the upper panel clarifies the resonance locations when
$m_{12}=50\rm{M}_{\odot}$. The region of interest where the outer orbit
potentially undergoing $e_{\star}-$excitation due to the apsidal precession
resonance is located within a wide range of $a_{\star}$ (or $P_{\star}$).

The lower panels of Figure 4 present the results of the averaged maximum
changes of the outer eccentricity ($\Delta e_{\star}^{\mathrm{peak}}$), over
all 100 runs, as a function of $a_{\star}$ for different masses of BHBs. A
remarkable eccentricity excitation is achieved for all tested values of the
BHB mass. For the $m_{12}=20\rm{M}_{\odot}$ case, the maximum eccentricity
growth is found for the smallest $a_{\star}$ and decreases for larger
$a_{\star}$. For the BHBs ($m_{12}/\rm{M_{\odot}}=50,100$), the perturbation
becomes stronger, and the induced peak eccentricity can be so high (because of
$e^{\mathrm{peak}}_{\star}\varpropto\mu_{\mathrm{in}}$) that some of the
systems with small $a_{\star}$ can become unbound. The peak of $\Delta
e_{\star}^{\mathrm{peak}}$ appears to shift to larger $a_{\star}$ for the
larger BHB masses. The decrease of $\Delta e^{\mathrm{peak}}$ with increasing
$a_{\star}$ is the result of the secular timescale of generating the
eccentricity growth (See Figure 3). Note that, in principle, the resonance can
occur as $a_{\star}\leq a_{\star}^{\mathrm{Res}}$ (the vertical blue line;
$P_{\star}\lesssim 160$ days). As shown, more systems can have larger $\Delta
e_{\star}^{\mathrm{peak}}$ when the evolving time is longer (See Figure 2). We
find that $e_{\star}-$excitation is mainly contributed by the inner BHBs with
relatively small mass ratios ($m_{2}/m_{1}\lesssim 0.5$).

The secular variability of the orbital eccentricity induces a change in the
projected orbit that may be probed via astrometric monitoring with surveys
such as Gaia Gaia . To determine detectability with Gaia, we compute a signal-
to-noise ratio (SNR) for astrometric detection,
$\rho=\theta_{\mathrm{signal}}/\theta_{\mathrm{Gaia}}$. Assuming the signal
can be well-approximated by the change in apocenter, the maximum such signal
for a system at distance $D$ is approximated by

$\theta_{\mathrm{signal}}(D)\simeq\frac{a_{\star}\Delta
e_{\star}^{\mathrm{peak}}}{D}.$ (10)

To compute noise, we follow Dan-Gaia , which draws on Casertano-1995 ;
Bernstein ; Casertano-2008 ; Perryman ; Ranalli  to evaluate the Gaia
astrometric precision for a single-scan. Here, we assume that our source is a
solar-type star with absolute V-band magnitude of $-26.8$, and $V-I_{c}=0.688$
Holmberg . Following (Perryman, ), the single-scan precision is computed from
the end-of-mission, sky-averaged parallax uncertainty
$\theta_{\mathrm{eom}}(D)$, as,

$\theta_{\mathrm{Gaia}}(D)\equiv\frac{\sqrt{140}}{1.1\times
2.15}\theta_{\mathrm{eom}}(D).$ (11)

The pre-factors account for an average 140 Gaia visits over 10 years, a
geometrical averaging factor of $2.15$, and a contingency margin of $1.1$
Perryman . Note that the contingency margin is 1.1 instead of the value of 1.2
chosen in (Perryman, ). This is based on newest information from Gaia EDR3
uncertainties as expressed in Section 1 of (GaiaFactSheet, ). To compute the
end-of-mission astrometric precision, we use the most up-to-date fitting
formula from the Gaia expected science performance document (GaiaFactSheet, ),

$\displaystyle\theta_{\mathrm{eom}}(D)$ $\displaystyle=$ $\displaystyle
0.527\left[40+800Z+30Z^{2}\right]^{1/2}\mu\mathrm{as}$ (12) $\displaystyle Z$
$\displaystyle\equiv$
$\displaystyle\max\left\\{10^{0.4\left[13.0-15.0\right]},10^{0.4\left[G-15.0\right]}\right\\}$
$\displaystyle G$ $\displaystyle=$ $\displaystyle
m_{V}(D)-0.01746+0.008092(V-I_{C})$
$\displaystyle-0.2810(V-I_{C})^{2}+0.03655(V-I_{C})^{3}.$

Here the $0.527$ prefactor is for a 10 year mission (referred to as Gaia DR5
in GaiaFactSheet ), the conversion between Gaia-G- and V-magnitudes, in the
last line above, is given in Table A2 of (Evans+2018, ), and
$m_{V}(D)=-26.8+2.5\log_{10}[(D/\mathrm{AU})^{2}]$ is the apparent magnitude
of a sun-like star at distance $D$ in Astronomical Units (AU).

Figure 5: Detectability of the change of the outer stellar eccentricity in the
($P_{\star}-D$) plane with Gaia, which is evaluated by Equations (10)-(11) and
$\Delta e_{\star}^{\mathrm{peak}}$ is given by the data from Figure 4 for 10
yrs. The subfigure in the top panel illustrates an example of detectability
where the source is at $D=900$ pc, and the shaded region corresponds to the
systems with $\theta_{\mathrm{signal}}\geq\theta_{\mathrm{Gaia}}$. The white
dashed lines specify the region with $\rho>2$.

Figure 5 summarizes the optimal-characteristic SNR as a function of distance
to the source and outer orbital period, for three different inner binary
masses. Here, we use the data from Figure 4 that gives the maximum change of
eccentricity for a range of $a_{\star}$ over 10 yrs (for $m_{12}=100\
\rm{M}_{\odot}$, we assume $\Delta e_{\star}^{\mathrm{peak}}\lesssim 0.3$ due
to the long term instability). We see that the overall SNR improves as the
distance $D$ decreases, and the boundary of detectability is marked at $\sim
10^{3}$ pc given by $\rho=2$ (Ranalli, ). Two peaks of high SNR are the
results of the pure $e_{\star}-$excitation (Figure 4) and the enhancement of
$\theta_{\mathrm{signal}}$ from wide binaries. Although the detectability here
is evaluated based only on $\Delta e_{\star}^{\mathrm{peak}}$, the outer
pericenter argument can vary in the actual detection Suto-2 . When the triple
system is inclined, a combination of predicted changes in orbital inclination
and line-of-sight orbital projection can also increase or decrease our
estimate here. Our results are characteristic of the optimal astrometric
signatures over the course of the Gaia mission.

Note that the average time between visits is $\sim 26$ days, hence, the
shorter period systems may be sampled at sub-orbital frequencies. However,
because this is a secular effect, even a longer than orbital sampling rate
could probe the longer-timescale secular change of the projected orbit.
Follow-up analysis must simulate mock systems to determine the required
orbital sampling rate and SNR needed to reliably detect this effect. Further
work is also required to determine in which cases the current Gaia pipeline
would flag such secularly evolving systems as binaries with an unseen
companion, or miss identify them due to difficulty in fitting to a Keplerian
orbit.

##  IV Summary and Discussion

We have studied a novel secular dynamical effect of a solar-type star around a
compact BHB in a nearly coplanar triple configuration. We point out that the
stellar orbit may undergo significant oscillations of the orbital eccentricity
if the system satisfies an “apsidal precession resonance”. Starting with an
eccentric inner BHB, the outer eccentricity can be excited due to the
resonance and the enhancement becomes stronger as the inner eccentricity
increases. This effect, which can drive the eccentricity of the outer orbits
close to unity (even unbound the binary), is overlooked in all previous
studies and may have applications for other types of systems, such as the
planet around a binary star and star/compact object around a supermassive BH
binary with/without a gaseous disc.

Figure 6: Similar to the top panel of Figure 4, but we include two
combinations of initial eccentricities (as labeled). Here, the mass of BHB is
set to be $m_{12}=50\ \rm{M}_{\odot}$ and the mass ratio is $m_{2}/m_{1}=0.2$.

Figure 7: Similar to the middle panel of Figure 1, but we include the change
of eccentricities of both inner and outer binaries. The system parameters are
$m_{12}=50\ \rm{M}_{\odot}$, $m_{2}/m_{1}=0.2$, $e_{\mathrm{in}}^{0}=0.9$,
$e_{\star}^{0}=0$, $P_{\star}=15$ days (left panel), $e_{\mathrm{in}}^{0}=0$,
$e_{\star}^{0}=0.9$, $P_{\star}=800$ days (right panel). Upper panels: the
ratio of apsidal precession rate. Lower panels: the results (black and red
lines) are obtained by the numerical integrations of SA equations for 25 yrs
(left panel) and $6\times 10^{3}$ yrs (right panel).

Note that the apsidal precession resonance allows angular momentum exchange to
occur efficiently between the inner and outer orbits, leading to a transfer of
eccentricity. One case that we do not include in our study is the triple
systems with $e_{\mathrm{in}}^{0}=0$ and $e_{\star}^{0}=0.9$. In this
situation, the inner (outer) binary is expected to become (less) eccentric as
the resonance occurs.

Figure 6 shows the parameter space where the apsidal precession resonance can
play a role, on the basis of
$\omega_{\mathrm{in}}(e_{\mathrm{in}})=\omega_{\star}(e_{\star})$. Compared to
the fiducial cases (left panel), the resonance region for the initially stable
systems in the right panel is shifted to larger $a_{\star}$ (or $P_{\star}$).
Therefore, the corresponding semimajor axis of the inner BHB
($a_{\mathrm{in}}$) becomes larger, leading to a much lower GW frequency range
(i.e., out of the LISA band).

Figure 7 presents the change of eccentricities due to resonance for both inner
and outer orbits, taking into account the examples identified in the Figure 6.
In the left panel (fiducial case), we see that the inner eccentricity
decreases (i.e., $\Delta e_{\mathrm{in}}^{\mathrm{max}}\lesssim 0.1$) when the
outer eccentricity is excited. In the right panel, we find that the
eccentricities evolve in an opposite way. The inner binary can efficiently
gain some eccentricity from the outer binary during the resonance (the peak
value is about $\Delta e_{\mathrm{in}}^{\mathrm{max}}\sim 0.5$), while the the
outer eccentricity can decrease by a factor of $\lesssim 0.2$.

Since we are interested in the LISA source and the evidence change of
$e_{\star}$ within $\sim 10$ years, we do not include this type of system
(with $e_{\mathrm{in}}^{0}=0$ and $e_{\star}^{0}=0.9$) in our main study.

The formation of compact massive BHB$+$star systems may be challenging. The
progenitor stars of these BHBs usually expand hundreds or thousand of solar
radii throughout their evolution, likely dynamically interacting with the
stellar tertiary. However, some low-metallicity massive stellar binaries might
remain compact throughout their evolution Marchant-2016 , allowing for
dynamically-stable compact triples Alejandro2021 . Alternatively, a BHB can be
formed first and eventually capture a long-lived low-mass star. Note that our
analysis is not restricted to any specific formation scenarios, and can be
adapted to other types of systems by applying scaling relations.

We find that the secular variability of the stellar orbit’s apocenter induced
by the changing eccentricity is detectable by Gaia for inner BHBs emitting GWs
in the LISA frequency range. Assuming that the formation and merger rates of
BHBs are in equilibrium, we expect to have hundreds of BHBs in the LISA band
in the Galaxy based on the LIGO detection rate XiaoFang ; Christian ; Wagg ;
Kremer ; Xianyu-LISA ; Hoang . Since our proposed secular variability can only
be resolved by Gaia within several kpc, the expected number of sources that
LISA and Gaia could see becomes $\sim$a few. Although the actual number of
LISA sources accompanied by a tertiary star in our Galaxy is quite uncertain,
identifying the secular motion of stellar orbits in current (Gaia Data Release
3; Gaia ) and future Gaia data is timely.

Our proof-of-concept calculations demonstrate that the long-term evolution of
eccentricity of a nearby stellar orbit can serve as a distinctive imprint of
such an unseen binary companion. Precise measurements of secular variability
are therefore an independent approach to reveal hidden BHBs, in addition to GW
detection. We emphasize that the inner binary systems which generate Gaia-
detectable variations in the orbits of their stellar tertiary also emit GWs in
the LISA band. Therefore, Gaia may provide candidates LISA sources before LISA
flies (planned for the 2030’s). In this sense, a joint detection with Gaia and
LISA Silvia-LISA ; Dan-LISA ; Breivik  would be a unique multi-messenger tool
to understand the evolution, fate and configurations of compact BHBs.

Acknowledgement.— We thank Dong Lai and Bence Kocsis for useful discussions.
We also thank Hagai Perets and Jihad Touma for careful reading the manuscript
and helpful comments given. B.L. gratefully acknowledges support from the
European Union’s Horizon 2021 research and innovation programme under the
Marie Sklodowska-Curie grant agreement No. 101065374. D.J.D. received funding
from the European Union’s Horizon 2020 research and innovation programme under
Marie Sklodowska-Curie grant agreement No. 101029157, and from the Danish
Independent Research Fund through Sapere Aude Starting Grant No. 121587.
A.V.G. received support through Villum Fonden grant No. 29466. J.S. is
supported by the Villum Fonden grant No. 29466.

## References

  * (1) R. Abbott et al. (LIGO Scientific, VIRGO), arXiv:2108.01045 [gr-qc] (2021).
  * (2) V. M. Lipunov, K. A. Postnov, and M. E. Prokhorov, Astron. Lett. 23, 492 (1997).
  * (3) V. M. Lipunov, V. Kornilov, E. Gorbovskoy, D. A. H. Buckley, N. Tiurina, P. Balanutsa, A. Kuznetsov, J. Greiner, V. Vladimirov, D. Vlasenko et al., Mon. Not. R. Astron. Soc. 465, 3656 (2017).
  * (4) P. Podsiadlowski, S. Rappaport, and Z. Han, Mon. Not. R. Astron. Soc. 341, 385 (2003).
  * (5) K. Belczynski, M. Dominik, T. Bulik, R. O’Shaughnessy, C. Fryer, and D. E. Holz, Astrophys. J. Lett. 715, L138 (2010).
  * (6) K. Belczynski, D. E. Holz, T. Bulik, and R. O’Shaughnessy, Nature (London) 534, 512 (2016).
  * (7) M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel, and R. O’Shaughnessy, Astrophys. J. 759, 52 (2012).
  * (8) M. Dominik, K. Belczynski, C. Fryer, D. E. Holz, E. Berti, T. Bulik, I. Mandel, and R. O’Shaughnessy, Astrophys. J. 779, 72 (2013).
  * (9) M. Dominik, E. Berti, R. O’Shaughnessy, I. Mandel, K. Belczynski, C. Fryer, D. E. Holz, T. Bulik, and F. Pannarale, Astrophys. J. 806, 263 (2015).
  * (10) S. Stevenson, A. Vigna-Gómez, I. Mandel, J. W. Barrett, C. J. Neijssel, D. Perkins, and S. E. de Mink, Nature Commun. 8, 14906 (2017).
  * (11) I. Mandel and S. E. De Mink, Mon. Not. R. Astron. Soc. 458, 2634 (2016).
  * (12) P. Marchant, N. Langer, P. Podsiadlowski, T. M. Tauris, and T. J. Moriya, Astron. Astrophys. 588, A50 (2016).
  * (13) L. du Buisson, P. Marchant, P. Podsiadlowski, C. Kobayashi, F. B. Abdalla, P. Taylor, I. Mandel, S. E. de Mink, T. J. Moriya, and N. Langer, Mon. Not. R. Astron. Soc. 499, 5941 (2020).
  * (14) J. Riley, I. Mandel, P. Marchant, E. Butler, K. Nathaniel, C. Neijssel, Shortt S., and A. Vigna-Gómez, Mon. Not. R. Astron. Soc. 505, 663 (2021).
  * (15) C. Baruteau, J. Cuadra, and D. N. C. Lin, Astrophys. J. 726, 28 (2011).
  * (16) B. McKernan, K. E. S. Ford, W. Lyra, and H. B. Perets, Mon. Not. Roy. Astron. Soc. 425, 460 (2012).
  * (17) B. McKernan, K. E. S. Ford, J. Bellovary, N. W. C. Leigh, Z. Haiman, B. Kocsis, W. Lyra, M. M. MacLow, B. Metzger, M. O’Dowd, S. Endlich, and D. J. Rosen, Astrophys. J. 866, 66 (2018).
  * (18) I. Bartos, B. Kocsis, Z. Haiman, and S. Márka, Astrophys. J. 835, 165 (2017).
  * (19) N. C. Stone, B. D. Metzger, and Z. Haiman, Mon. Not. Roy. Astron. Soc. 464, 946 (2017).
  * (20) N. W. C. Leigh, A. M. Geller, B. McKernan, K. E. S. Ford, M. M. Mac Low, J. Bellovary, Z. Haiman, W. Lyra, J. Samsing, M. O’Dowd, B. Kocsis, and S. Endlich, Mon. Not. Roy. Astron. Soc. 474, 5672 (2018).
  * (21) A. Secunda, J. Bellovary, M.-M. Mac Low, K. E. S. Ford, B. McKernan, N. W. C. Leigh, W. Lyra, and Z. Sándor, Astrophys. J. 878, 85 (2019).
  * (22) Y. Yang, I. Bartos, V. Gayathri, K. E. S. Ford, Z. Haiman, S. Klimenko, B. Kocsis, S. Márka, Z. Márka, B. McKernan, and R. O’Shaughnessy, Phys. Rev. Lett. 123, 181101 (2019).
  * (23) M. Gröbner, W. Ishibashi, S. Tiwari, M. Haney, and P. Jetzer, Astron. Astrophys. 638, A119 (2020).
  * (24) W. Ishibashi and M. Gröbner, Astron. Astrophys. 639, A108 (2020).
  * (25) H. Tagawa, Z. Haiman, and B. Kocsis, Astrophys. J. 898, 25 (2020).
  * (26) Y. P. Li, A. M. Dempsey, S. Li, H. Li, and J. Li, Astrophys. J. 911, 124 (2021).
  * (27) K. E. S. Ford and B. McKernan, arXiv:2109.03212 [astro-ph.HE] (2021).
  * (28) J. Samsing, I. Bartos, D. J. D’Orazio, Z. Haiman, B. Kocsis, N. W. C. Leigh, B. Liu, M. E. Pessah, and H. Tagawa, Nature, 603, 237 (2022).
  * (29) R. Li and D. Lai, arXiv:2202.07633 [astro-ph.HE] (2022).
  * (30) J. Li and D. Lai, arXiv:2203.05584 [astro-ph.HE] (2022).
  * (31) S. F. P. Zwart and S. L. W. McMillan, Astrophys. J. Lett. 528, L17 (2000).
  * (32) R. M. O’Leary, F. A. Rasio, J. M. Fregeau, N. Ivanova, and R. OShaughnessy, Astrophys. J. 637, 937 (2006).
  * (33) M. C. Miller and V. Lauburg, Astrophys. J. 692, 917 (2009).
  * (34) S. Banerjee, H. Baumgardt, and P. Kroupa, Mon. Not. R. Astron. Soc. 402, 371 (2010).
  * (35) J. M. B. Downing, M. J. Benacquista, M. Giersz, and R. Spurzem, Mon. Not. R. Astron. Soc. 407, 1946 (2010).
  * (36) B. M. Ziosi, M. Mapelli, M. Branchesi, and G. Tormen, Mon. Not. R. Astron. Soc. 441, 3703 (2014).
  * (37) C. L. Rodriguez, M. Morscher, B. Pattabiraman, S. Chatterjee, C.-J. Haster, and F. A. Rasio, Phys. Rev. Lett. 115, 051101 (2015).
  * (38) J. Samsing and E. Ramirez-Ruiz, Astrophys. J. Lett. 840, L14 (2017).
  * (39) J. Samsing and D. J. D’Orazio, Mon. Not. R. Astron. Soc. 481, 5445 (2018).
  * (40) C. L. Rodriguez, P. Amaro-Seoane, S. Chatterjee, and F. A. Rasio, Phys. Rev. Lett. 120, 151101 (2018).
  * (41) L. Gondán, B. Kocsis, P. Raffai, and Z. Frei, Astrophys. J. 860, 5 (2018).
  * (42) H. von Zeipel, Astronomische Nachrichten 183, 345 (1910).
  * (43) M. L. Lidov, Planetary and Space Science 9, 719 (1962).
  * (44) Y. Kozai, Astron. J. 67, 591 (1962).
  * (45) S. Naoz, Annu. Rev. Astron. Astrophys. 54, 441 (2016).
  * (46) M. C. Miller and D. P. Hamilton, Astrophys. J. 576, 894 (2002).
  * (47) L. Wen, Astrophys. J. 598, 419 (2003).
  * (48) F. Antonini and H. B. Perets, Astrophys. J. 757, 27 (2012).
  * (49) F. Antonini, S. Toonen, and A. S. Hamers, Astrophys. J. 841, 77 (2017).
  * (50) K. Silsbee and S. Tremaine, Astrophys. J. 836, 39 (2017).
  * (51) C. Petrovich and F. Antonini, Astrophys. J. 846, 146 (2017).
  * (52) B. Liu and D. Lai, Astrophys. J. 863, 68 (2018).
  * (53) L. Randall and Z.-Z. Xianyu, Astrophys. J. 853, 93 (2018).
  * (54) B.-M. Hoang, S. Naoz, B. Kocsis, F. A. Rasio, and F. Dosopoulou, Astrophys. J. 856, 140 (2018).
  * (55) B. Liu and D. Lai, Mon. Not. R. Astron. Soc. 483, 4060 (2019).
  * (56) G. Fragione and B. Kocsis, Mon. Not. R. Astron. Soc. 486, 4781 (2019).
  * (57) G. Fragione, E. Grishin, N. W. C. Leigh, H. B. Perets, and R. Perna, Mon. Not. R. Astron. Soc. 488, 47 (2019).
  * (58) M. Zevin, J. Samsing, C. Rodriguez, C.-J. Haster, and E. Ramirez-Ruiz, Astrophys. J. 871, 91 (2019).
  * (59) B. Liu and D. Lai, Mon. Not. R. Astron. Soc. 502, 2049 (2021).
  * (60) E. Michaely and H. B. Perets, Astrophys. J. Lett. 887, L36 (2019).
  * (61) E. Michaely and H. B. Perets, Mon. Not. R. Astron. Soc. 498, 4924 (2020).
  * (62) P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy et al., arXiv:1702.00786 [astro-ph.IM] (2017).
  * (63) J. Luo, L.-S. Chen, H.-Z. Duan, Y.-G. Gong, S. Hu, J. Ji, Q. Liu, J. Mei, V. Milyukov, M. Sazhin, C.-G. Shao, V. T. Toth, H.-B. Tu, Y. Wang, Y. Wang, H.-C. Yeh, M.-S. Zhan, Y. Zhang, V. Zharov, and Z.-B. Zhou, Classical and Quantum Gravity 33, 035010 (2016).
  * (64) W.-R. Hu and Y.-L. Wu, National Science Review 4, 685 (2017).
  * (65) T. Nakamura, M. Ando, T. Kinugawa, H. Nakano, K. Eda, S. Sato, M. Musha, T. Akutsu, T. Tanaka, N. Seto, N. Kanda, and Y. Itoh, Progress of Theoretical and Experimental Physics 2016, 093E01 (2016).
  * (66) M. Arca Sedda, C. P. L. Berry, K. Jani, P. AmaroSeoane, P. Auclair, J. Baird, T. Baker, E. Berti et al., arXiv:1908.11375 [gr-qc] (2019).
  * (67) K. A. Kuns, H. Yu, Y. Chen, and R. X. Adhikari, Phys. Rev. D 102, 043001 (2020).
  * (68) A. Tokovinin, S. Thomas, M. Sterzik, and S. Udry, Astron. Astrophys. 450, 681 (2006).
  * (69) D. Raghavan, H. A. McAlister, T. J. Henry, D. W. Latham, G. W. Marcy, B. D. Mason, D. R. Gies, R. J. White, and T. A. ten Brummelaar, Astrophys. J. Suppl. Ser. 190, 1 (2010).
  * (70) K. Fuhrmann, R. Chini, L. Kaderhandt, and Z. Chen, Astrophys. J. 836, 139 (2017).
  * (71) S. Naoz, G. Li, M. Zanardi, G. C. de Elía, and R. P. Di Sisto, Astron. J. 154, 18 (2017).
  * (72) B. R. Vinson and E. Chiang, Mon. Not. R. Astron. Soc. 474, 4855 (2018).
  * (73) G. C. de Elía, M. Zanardi, A. Dugaro, S. Naoz, Astron. Astrophys. 627, A17 (2019).
  * (74) T. Hayashi, S. Wang, and Y. Suto, Astrophys. J. 890, 112 (2020).
  * (75) T. Hayashi and Y. Suto ,Astrophys. J. 897, 29 (2020).
  * (76) T. Hayashi and Y. Suto, Astrophys. J. 907, 48 (2021).
  * (77) S. M. Ransom, I. H. Stairs, A. M. Archibald, J. W. T. Hessels, D. L. Kaplan, M. H. van Kerkwijk, J. Boyles, A. T. Deller, S. Chatterjee, and A. Schechtman-Rook et al., Nature 505, 520 (2014).
  * (78) T. A. Thompson, C. S. Kochanek, K. Z. Stanek, C. Badenes, R. S. Post, T. Jayasinghe, D. W. Latham, A. Bieryla, G. A. Esquerdo, P. Berlind, M. L. Calkins, J. Tayar, L. Lindegren, J. A. Johnson, T. W. S. Holoien, K. Auchettl, and K. Covey, Science 366, 637 (2019)
  * (79) N. L. Eisner, C. Johnston, S. Toonen, A. J. Frost, S. Janssens, C. J. Lintott, S. Aigrain, et al., Mon. Not. R. Astron. Soc. 511, 4710 (2022).
  * (80) B. Liu, D. Lai, and Y.-F. Yuan, Phys. Rev. D 92, 124048 (2015).
  * (81) B. Liu, and D. Lai, Phys. Rev. D 102, 023020 (2020).
  * (82) E. B. Ford, B. Kozinsky, and F. A. Rasio, Astrophys. J. 535, 385 (2000).
  * (83) S. Naoz, B. Kocsis, A. Loeb, and N. Yunes, Astrophys. J. 773, 187 (2013).
  * (84) L. Blanchet, Living Rev. Rel. 17, 2 (2014).
  * (85) L. G. Kiseleva, S. J. Aarseth, P. P. Eggleton, and R. de La Fuente Marcos, in ASP Conf. Ser. 90, The Origins, Evolution, and Destinies of Binary Stars in Clusters, ed. E. F. Milone and J.-C. Mermilliod (San Francisco, CA: ASP), 433 (1996).
  * (86) C. D. Murray, and S. F. Dermott, Solar System Dynamics, Cambridge U. Press, NY (1999).
  * (87) Y. Wu, and P. Goldreich, Astrophys. J. 564, 1024 (2002).
  * (88) M. H.Lee, S. J. Peale, Astrophys. J. 592, 1201 (2003).
  * (89) B. Liu, D. J. Muñoz, and D. Lai, Mon. Not. Roy. Astron. Soc. 447, 747 (2015).
  * (90) P. C. Peters, Phys. Rev. 136, B1224 (1964).
  * (91) Gaia Collaboration, F. Arenou, C. Babusiaux, M. A. Barstow, S. Faigler, A. Jorissen, P. Kervella, et al., arXiv:2206.05595[astro-ph.SR] (2022).
  * (92) D. J. D’Orazio and A. Loeb, Phys. Rev. D, 100, 103016 (2019).
  * (93) S. Casertano, M. G. Lattanzi, and M. A. C. Perryman, in Future Possibilities for Astrometry in Space, ESA Special Publication, Vol. 379, edited by M. A. C. Perryman and F. van Leeuwen (1995) pp. 47–54.
  * (94) H.-H. Bernstein and U. Bastian, in Future Possibilities for Astrometry in Space, ESA Special Publication, Vol. 379, edited by M. A. C. Perryman and F. van Leeuwen (1995) p. 55.
  * (95) S. Casertano, M. G. Lattanzi, A. Sozzetti, A. Spagna, S. Jancart, R. Morbidelli, R. Pannunzio, D. Pourbaix, and D. Queloz, Astron. Astrophys. 482, 699 (2008).
  * (96) M. Perryman, J. Hartman, G. Á. Bakos, and L. Lindegren, Astrophys. J. 797, 14 (2014).
  * (97) P. Ranalli, D. Hobbs, and L. Lindegren, Astron. Astrophys. 614, A30 (2018).
  * (98) J. Holmberg, C. Flynn, and L. Portinari, Mon. Not. Roy. Astron. Soc. 367, 449 (2006).
  * (99) Gaia collaboration, Expected Science Performance for the nominal and the extended mission, https://www.cosmos.esa.int/web/gaia/science-performance.
  * (100) Gaia Collaboration, A. G. A. Brown, A. Vallenari, T. Prusti, J. H. J. de Bruijne, C. Babusiaux, C. A. L. Bailer-Jones, M. Biermann, D. W. Evans, L. Eyer et al., Astron. Astrophys. 616, A1 (2018).
  * (101) A. Vigna-Gómez, S. Toonen, E. Ramirez-Ruiz, N. W. C. Leigh, J. Riley, and C.-J. Haster, Astrophys. J. Lett. 907, L19 (2021).
  * (102) X. Fang, T. A. Thompson, and C. M. Hirata, Astrophys. J. 875, 75 (2019).
  * (103) P. Christian and A. Loeb, Mon. Not. Roy. Astron. Soc. 469, 930 (2017).
  * (104) T. Wagg, F. S. Broekgaarden, S. E. de Mink, L. A. C. van Son, N. Frankel, S. Justham, arXiv:2111.13704[astro-ph.HE] (2021).
  * (105) K. Kremer, S. Chatterjee, K. Breivik, C. L. Rodriguez, S. L. Larson, and F. A. Rasio, Phys. Rev. Lett. 120, 191103 (2018).
  * (106) L. Randall and Z.-Z. Xianyu, arXiv:1902.08604[astroph.HE] (2019).
  * (107) B.-M. Hoang, S. Naoz, B. Kocsis, M. F. Will, and J. Mclver, Astrophys. J. Lett. 875, L31 (2019).
  * (108) T. Robson, N. J. Cornish, N. Tamanini, and S. Toonen, Phys. Rev. D 98,064012 (2018)
  * (109) D. J. D’Orazio and J. Samsing, Mon. Not. R. Astron. Soc. 481, 4775 (2018).
  * (110) K. Breivik, K. Kremer, M. Bueno, S. L. Larson, S. Coughlin and V. Kalogera, Astrophys. J. Lett. 854, L1 (2018).


# Long-range Velocity Correlations from Active Dopants

Leila Abbaspour  Institute for the Dynamics of Complex Systems, University of
Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany.  Max Planck
School Matter to Life, University of Göttingen, Friedrich-Hund-Platz 1, 37077
Göttingen, Germany.  Max Planck Institute for Dynamics and Self-Organization,
Am Fassberg 17, 37077 Göttingen, Germany.  Rituparno Mandal  Institute for
Theoretical Physics, University of Göttingen, Friedrich-Hund-Platz 1, 37077
Göttingen, Germany  Peter Sollich  Institute for Theoretical Physics,
University of Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany
Department of Mathematics, King’s College London, London WC2R 2LS, UK  Stefan
Klumpp  Institute for the Dynamics of Complex Systems, University of
Göttingen, Friedrich-Hund-Platz 1, 37077 Göttingen, Germany.  Max Planck
School Matter to Life, University of Göttingen, Friedrich-Hund-Platz 1, 37077
Göttingen, Germany.

###### Abstract

One of the most remarkable observations in dense active matter systems is the
appearance of long-range velocity correlations without any explicit aligning
interaction (of e.g. Vicsek type). Here we show that this kind of long range
velocity correlation can also be generated in a dense athermal passive system
by the inclusion of a very small fraction of active Brownian particles. We
develop a continuum theory to explain the emergence of velocity correlations
generated via such active dopants. We validate the predictions for the effects
of magnitude and persistence time of the active force and the area fractions
of active or passive particles using extensive Brownian dynamics simulation of
a canonical active-passive mixture. Our work decouples the roles that density
and activity play in generating long range velocity correlations in such
exotic non-equilibrium steady states.

### Introduction:

Active matter systems are one of the best-known examples of non-equilibrium
systems and are famous for their fascinating collective behaviour across a
diverse range of length and time scales Marchetti  _et al._ (2013); Bechinger
_et al._ (2016), from the cytoskeleton to bacterial colonies, tissues, flocks
of birds to animal herds.

Systems of active Brownian particles (ABP), i.e. particles exhibiting self-
propulsion, are a canonical example of active matter Fily and Marchetti
(2012); Takatori and Brady (2015); Cates and Tailleur (2015); Levis  _et al._
(2017); Solon  _et al._ (2018). These systems exhibit two types of non-
equilibrium pattern formation: in the presence of aligning interactions
between the directions of the self-propelled motion of the particles, they
show flocking, i.e. collective directed motion of groups of such particles
Mora  _et al._ (2016); Ballerini  _et al._ (2008). In the absence of such
aligning interactions, they exhibit motility induced phase separation (MIPS),
crucially without the need for attractive interactions Fily and Marchetti
(2012); Cates and Tailleur (2015); Takatori and Brady (2015). It is important
to note that this phase separation is not associated with any macroscopic
order in the orientation of the self-propulsion directions of the particles.
With external driving, on the other hand, such order can appear, see e.g. Ref.
Mandal and Sollich (2021).

Figure 1: Schematic illustration of (a) a dense assembly of active particles,
and (b) a dense mixture of active (blue) and passive (red) particles.
Snapshots of the system with each particle coloured according to the
orientation of its instantaneous velocity vector, for (c) $\phi_{a}=0.7$
($\phi_{p}=0$), (d) $\phi_{a}=0.01$, $\phi_{p}=0.69$; both systems show strong
emergent velocity correlations.

Very recently it has been discovered that while systems without aligning
interactions show no macroscopic orientational ordering, there are
spectacularly large spatial structures in the instantaneous velocity field,
especially within the dense phase created by motility-induced phase separation
(MIPS) Caprini  _et al._ (2020a); Caprini and Marini Bettolo Marconi (2021);
Caprini  _et al._ (2020b); Henkes  _et al._ (2020); Szamel and Flenner (2021).
It has been shown using both analytical calculations and numerical models (and
confirmed in experiments in dense tissues Henkes  _et al._ (2020)) that in
such a scenario, a dense assembly of active particles generates long range
velocity correlations in the large persistence time limit; the corresponding
correlation length grows as a power law $\sim\sqrt{\tau_{p}}$ with increasing
persistence time $\tau_{p}$ Caprini  _et al._ (2020a); Caprini and Marini
Bettolo Marconi (2021); Caprini  _et al._ (2020b); Henkes  _et al._ (2020);
Szamel and Flenner (2021). The emergence of such non-equilibrium velocity
correlation has always been attributed to a (highly persistent) dense active
matter system and thus taken to require a high density of active particles.

These results raise the question of whether the above two conditions, of high
density and high activity, can be decoupled. In particular, could long-range
velocity correlations be generated in a dense system of passive particles, by
introducing activity only through a small fraction (e.g. much lower than the
percolation density) of active particles? A related question of interest is
how ordered non-equilibrium states are affected by the inclusion of defect
particles (e.g. static defects or motile non-aligning agents known as
dissenters) that do not participate in the processes that induce the order.
The role of both quenched Yllanes  _et al._ (2017); Das  _et al._ (2018);
Martinez  _et al._ (2018) and annealed Yllanes  _et al._ (2017); Martinez  _et
al._ (2018); Bera and Sood (2020) disorder in the context of Vicsek-like
models has been investigated very recently and it has been shown that presence
of both types of disorder tends to destroy the ordered flocking state Yllanes
_et al._ (2017); Das  _et al._ (2018); Martinez  _et al._ (2018); Bera and
Sood (2020). Our work explores a similar line of questions, but in a system
where long-range velocity correlations appear without any explicit alignment
interaction. We ask in particular whether the long-range order in the
instantaneous velocity field seen in such systems is stabilised or
destabilized by the inclusion of a large fraction of passive particles.

Motivated by the above questions, in this paper we study mixtures of active
and passive particles to explore whether inclusion of passive particles
enhances or suppresses local orientational order and to see whether we can
decouple the roles of activity and density in generating long-range velocity
correlations. Using extensive particle-based simulation of an active-passive
mixture, we demonstrate that velocity ordering is enhanced by an increasing
density of passive particles, and show that long-range velocity correlations
can be generated in an athermal passive medium by a tiny fraction of active
insertions (dopants) as long as the medium is dense enough. We also construct
an analytical theory to explain the physics of velocity correlations in a
dense passive medium with active dopants. Our hydrodynamic theory predicts
that the amplitude of the velocity correlations is proportional to $f^{2}$
where $f$ is the magnitude of the propulsion force acting on each active
particle, proportional to $\tau_{p}^{-3/4}$ where $\tau_{p}$ is the
persistence time of the active particles, and proportional to the density of
active particles $\phi_{a}$. The hydrodynamic theory also predicts that the
correlation length $\xi_{L}$ only depends on $\tau_{p}$, as $\sqrt{\tau_{p}}$
Caprini  _et al._ (2020a); Caprini and Marini Bettolo Marconi (2021); Caprini
_et al._ (2020b); Henkes  _et al._ (2020); Szamel and Flenner (2021), not on
the active forcing magnitude $f$. We verify these theoretical predictions by
performing further targeted simulations. The explicit form of the correlation
function that we derive theoretically decouples the roles that density and
activity play in generating long range velocity correlations in a non-
equilibrium steady state. This insight will be useful in understanding long
range ordering in e.g. dense passive colloidal systems driven by a few self-
propelled Janus colloids, or assembly of dead bacteria churned up by few
living ones.

### Particle-based Model:

We consider a binary mixture Abbaspour and Klumpp (2021) of passive and active
Brownian particles (ABP) Fily and Marchetti (2012); Takatori and Brady (2015);
Levis  _et al._ (2017); Solon  _et al._ (2018); Cates and Tailleur (2015)
moving in two dimensions and occupying area fractions of $\phi_{p}$ and
$\phi_{a}$, respectively. The dynamical evolution of the particle positions is
described by the overdamped equations of motion:

$\gamma\dot{\mathbf{r}}_{i}=\mathbf{F}_{i}+f\mathbf{n}_{i}\,\Delta_{i}(\mathcal{A})$
(1)

where $\mathbf{r}_{i}$ is the position vector of the $i$-th particle and
$\gamma$ is the constant drag coefficient governing the friction force acting
on each particle. The factor $\Delta_{i}(\mathcal{A})$ ($\Delta_{i}=1$ for the
active particles and $\Delta_{i}=0$ for the passive particles) restricts the
active forces $f\mathbf{n}_{i}$ to the particles in the subset $\mathbcal{A}$
of active particles. The orientation vectors of the active forces are
$\mathbf{n}_{i}=(\cos{\theta_{i}},\sin{\theta_{i}})$, and the orientation
angles $\theta_{i}$ of the active forcing follow the dynamics

$\mathbf{\dot{\theta}}_{i}=\sqrt{\frac{2}{\tau_{p}}}\zeta_{i}\quad\text{for}\quad
i\in\mathbcal{A}$ (2)

where the noise $\zeta_{i}$ has zero mean and time correlations
$\langle\zeta_{i}(t)\zeta_{j}(t^{\prime})\rangle=\delta_{ij}\delta(t-t^{\prime}).$
In Eq. 2, $\tau_{p}$ is the persistence time. All particles in the system
interact only through steric interactions described by the forces
$\mathbf{F}_{i}=-\nabla_{i}U$ where $U$ is a repulsive WCA (Weeks-Chandler-
Anderson) interaction potential Weeks  _et al._ (1971) (see supplementary
information for the details of the potential).

Figure 2: Snapshot of the system (binary mixture of active and passive
particles) for various combinations of values of the area fractions of active
and passive particles, $\phi_{a}$ and $\phi_{p}$, respectively. The snapshots
at the bottom right corner shows that a small fraction of active inclusions
can cause long-range velocity correlations in a dense athermal and almost
entirely passive system.

### Velocity correlations in active–passive mixture:

First we reproduce the long range velocity correlations in a completely active
system (see Fig.1(a) for a schematic and Fig.1(b) for a snapshot of the system
showing the long range ordering of the instantaneous velocities). This effect
has been reported before in the context of different _dense active matter_
systems, both in simulations and in experiments Caprini  _et al._ (2020a);
Caprini and Marini Bettolo Marconi (2021); Caprini  _et al._ (2020b); Henkes
_et al._ (2020); Szamel and Flenner (2021). Here we want to explore a
different scenario (see Fig.1(c) for a schematic) and test whether a dense
_passive_ system that is driven by just a few active particles can show
similar ordering. Indeed, Fig.1(d) shows remarkably similar order in such a
dense passive system that is driven by a very small fraction ($\phi_{a}=0.01$)
of active Brownian particles. To explore this further, we also ran simulations
for combinations of different fractions of active and passive particles (see
Fig. 2 for snapshots). Strong local velocity correlations are seen to emerge
as long as the total density $\phi_{\rm tot}=\phi_{a}+\phi_{p}$ is high enough
and there exists a non-zero fraction of active particles. To get a better
understanding of this velocity ordering in an active-passive mixture we
developed a hydrodynamic theory that we describe in the next paragraph.

Figure 3: (a) Scaled velocity-velocity spatial auto correlation as a function
of spatial distance $r$, for different active force magnitudes $f$ as
indicated in the legend, and $\tau_{p}=3$ (b) The correlation length $\xi_{L}$
as a function of the active force for different persistence times $\tau_{p}$
as shown; the plot demonstrates that, at fixed persistence time, the
correlation length is independent of the magnitude of the active force.

### Theory:

To derive the correlation of the hydrodynamic velocity field of the system, we
extend the approach used by Henkes et al. for a dense active system Henkes
_et al._ (2020) to the case of a passive system with active dopants. We
consider the entire collection of passive particles as a dense medium having a
smooth velocity field, with the active particles as random point-like defects
that are the sources of a force density of the form

$\mathbf{F}_{i}(\mathbf{r},t)=a^{2}f\mathbf{n}_{i}(t)\delta(\mathbf{r}-\mathbf{r}_{i}(t))$
(3)

for the $i$-th active particle. Here $a$ is a microscopic length scale given
by the particle size, $f$ is the magnitude of the propulsion force acting on
each active particle as before, and $\textbf{n}_{i}$ is the unit vector
associated with the orientation of this active force. Assuming now that in the
large persistence time limit and at sufficiently high density, the active
particles deform the elastic solid-like medium by pushing other particles
without changing their positions significantly, we can evaluate the
correlation between the forces from active particles $i$ and $j$ in Fourier
space ($\mathbf{q},\omega$) as

$\langle\mathbf{F}_{i}(\mathbf{q},\omega)\cdot\mathbf{F}_{j}(\mathbf{q^{\prime}},\omega^{\prime})\rangle=\frac{4\pi\tau_{p}a^{4}{f}^{2}}{1+(\tau_{p}\omega)^{2}}e^{i\mathbf{q}\cdot\mathbf{r}_{i}}e^{i\mathbf{q}^{\prime}\cdot\mathbf{r}_{j}}\delta_{ij}\delta(\omega+\omega^{\prime})$
(4)

where the frequency dependence arises from the dynamics of the active force
orientation $\mathbf{n}_{i}(t)$ Henkes  _et al._ (2020). We then sum over all
the particles (over the indices $i$ and $j$) and also over the steady state
probability distribution of the active particles’ positions $\mathbf{r}_{i}$,
which we take as uniform, to arrive at the total force correlator

$\langle\mathbf{F}(\mathbf{q},\omega)\cdot\mathbf{F}(\mathbf{q}^{\prime},\omega^{\prime})\rangle=\frac{N_{a}a^{2}f^{2}(2\pi)^{3}}{N}\delta(\mathbf{q}+\mathbf{q^{\prime}})\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}$
(5)

where $N_{a}$ is the number of active particles and $N$ the total number of
particles in the system. After taking an angular average, the form of the
velocity correlation function can be derived Henkes  _et al._ (2020) (see
supplementary information for the full calculation and the final form of the
correlation function). Its behaviour can be approximated at large distances as

$C(\mathbf{r})\approx\left(\frac{\phi_{a}}{\phi_{a}+\phi_{p}}\right)\left(\frac{a^{2}{f}^{2}}{4\pi\zeta^{2}}\sqrt{\frac{\pi}{2r}}\right)\frac{1}{{\xi_{L}}^{3/2}}e^{-{r/\xi_{L}}}$
(6)

where $\xi_{L}=\left(\frac{(B+\mu)\tau_{p}}{\zeta}\right)^{1/2}$, $r$ is the
spatial distance, $\zeta$ is the friction coefficient as before, and $B$ and
$\mu$ are the bulk and shear moduli of the overall medium. These are expected
to be dependent only on the total area fraction $\phi_{\rm
tot}=\phi_{a}+\phi_{p}$ of active and passive particles. Therefore we now have
a testable prediction about the equal time velocity auto correlation function
from our hydrodynamic theory in terms of the control parameters
$f,\tau_{p},\phi_{a},\phi_{p}$ etc.

Figure 4: (a) The scaled equal-time velocity-velocity correlation as a
function of spatial distance $r$ for different persistence times $\tau_{p}$ as
indicated, for active force magnitude $f=0.25$. (b) The correlation length
$\xi_{L}$ as a function of the persistence time $\tau_{p}$ for different
values of the active force magnitude $f$ shows that the correlation length
grows as a power law (with exponent $\frac{1}{2}$) with the persistence time
$\tau_{p}$ and this behaviour is independent of the magnitude of the active
force $f$, which is shown in the legend.

### Comparison with Theory:

To validate the predictions of our theory we ran further simulations to test
explicitly the effects of varying active force magnitude $f$, persistence time
$\tau_{p}$, and area fractions of active $\phi_{a}$ and passive particles
$\phi_{p}$, and compared the correlation functions from those simulations with
those calculated using the hydrodynamic theory.

Eq. (6) predicts that the prefactor of the velocity-velocity spatial auto-
correlation function will scale as $f^{2}$ when we vary the active force
magnitude $f$, without any associated variation in the correlation length.
Therefore the values of the scaled auto-correlation function $C(r)/f^{2}$ are
expected to collapse into a single curve for different values of the active
forcing $f$ as long as the other parameters are kept constant. Fig. 3(a)
clearly shows that the scaled correlation functions for different magnitudes
of the active force do indeed collapse on top of each other. In Fig. 3(b) we
show the correlation length ($\xi_{L}$) as a function of the active force $f$
for different persistence times ($\tau_{p}$), which provides further evidence
that the correlation length is independent of the magnitude of the active
force when the persistence time is kept constant.

To shed light on the effects of the persistence time scale we vary $\tau_{p}$
next, keeping the both area fractions ($\phi_{a},\phi_{p}$) and the active
forcing magnitude $f$ constant. Our hydrodynamic theory (see Eq. 6) suggests
that the prefactor of the velocity-velocity spatial auto-correlation function
scales as $\tau_{p}^{-3/4}$ (due to the dependence of $\xi_{L}$ on
$\tau_{p}$), apart from the standard exponential dependence on $r/\xi_{L}$.
Indeed as Fig. 4(a) shows the simulation data points for different persistence
times $\tau_{p}$ for a given active force $f$ fall nicely on the same curve
once we scale $C(r)$ appropriately, i.e. $C$ by $\tau_{p}^{-3/4}$ and $r$ by
$\xi_{L}$. Fig. 4(b) further indicates that the correlation length $\xi_{L}$
grows as the square root of persistence time $\tau_{p}$ regardless of the
magnitude of the active force. This is also consistent with our theory and the
earlier studies Caprini  _et al._ (2020a); Caprini and Marini Bettolo Marconi
(2021); Caprini  _et al._ (2020b); Henkes  _et al._ (2020); Szamel and Flenner
(2021).

We finally explore the dependence on the area fractions of passive and active
particles $\phi_{p}$ and $\phi_{a}$, respectively, while keeping the active
forcing parameters ($f,\tau_{p}$) constant. The theory suggests that apart
from the linear dependence on the fraction of active particles $\phi_{a}$
there is no separate dependence on $\phi_{a}$, i.e. all other density
dependences of the correlation function $C(r)$ appear only via the total area
fraction $\phi_{\rm tot}=\phi_{a}+\phi_{p}$ of the binary mixture. Our
simulations, which involve different mixture compositions, confirm this
prediction in Fig. 5 (a) where we scale the correlation function by the area
fraction of active particles $\phi_{a}$ for a fixed value of the total density
$\phi_{\rm tot}$. Fig. 5 (b) shows that the system has practically the same
correlation length regardless of the fraction $\phi_{a}$ of active particles
(or the fraction $\phi_{p}$ of passive particles) when the total area fraction
$\phi_{\rm tot}$ is sufficiently high and kept constant.

Figure 5: (a) The scaled velocity-velocity spatial auto correlation function
as a function of spatial distance $r$ for different combinations of $\phi_{a}$
and $\phi_{p}$ for a fixed value of $\phi_{\rm tot}=0.95$. (b) Data for the
correlation length are consistent with $\xi_{L}$ being constant across the
same combinations of $\phi_{a}$ and $\phi_{p}$; the color scheme is the same
as in (a).

### Conclusion:

In this article we have demonstrated that long-range velocity correlations
(which have only been observed in dense active mater system until now Caprini
_et al._ (2020a); Caprini and Marini Bettolo Marconi (2021); Caprini  _et al._
(2020b); Henkes  _et al._ (2020); Szamel and Flenner (2021)) can be generated
in a dense athermal passive system by including a very small fraction of
persistent active Brownian particles. This observation conceptually decouples
the roles played by density and activity parameters in generating such non-
equilibrium ordering effects. Also, our results extend the discussion on
whether inclusion of disorder can increase order in a system or whether
conversely it tends to destroy order Yllanes  _et al._ (2017); Das  _et al._
(2018); Martinez  _et al._ (2018); Bera and Sood (2020).

We started by providing evidence that with a very small amount of active
inclusions or dopants, an otherwise passive, dense athermal system can exhibit
long range velocity correlations similar to a pure dense assembly of active
particles. We explored the degree of velocity correlation for different
numbers of active and passive particles particles in such a mixture. We then
derived the hydrodynamic theory to calculate the equal time velocity auto-
correlation function in terms of the microscopic system parameters such as
$f$, $\tau_{p}$ etc. This theory made testable predictions that we confirmed
via further molecular dynamics simulations. We examined the impact of
different parameters on the velocity correlations and found good agreement
between the simulation results and the hydrodynamic theory. Specifically, we
found that the correlation length depends only on the overall area fraction of
active and passive particles and grows as $\sqrt{\tau_{p}}$ with the
persistence time of self-propulsion. The latter result is in agreement with
previous findings on purely active systems Caprini  _et al._ (2020a); Caprini
and Marini Bettolo Marconi (2021); Caprini  _et al._ (2020b); Henkes  _et al._
(2020); Szamel and Flenner (2021).

Our predictions and results can be further tested both in simulations Stürmer
_et al._ (2019); Abbaspour and Klumpp (2021); Banerjee  _et al._ (2022) and in
experiments, e.g. on mixtures of microbes and passive colloidal particles
Williams  _et al._ (2022), assemblies of active and passive colloids Singh
_et al._ (2017); Mu  _et al._ (2022), mixture of mobile and immobile bacteria
Patteson  _et al._ (2018), or active granular mixtures Kumar  _et al._ (2014);
Gupta  _et al._ (2022). Understanding the decoupling of density and activity
makes it possible not only to reproduce a long-range velocity correlation in
different active-passive mixtures or assemblies but also paves the way for
designing and controlling active matter for practical purposes, e.g. in the
context of transport and mixing.

## Acknowledgement

This research was conducted within the Max Planck School Matter to Life,
supported by the German Federal Ministry of Education and Research (BMBF) in
collaboration with the Max Planck Society. The simulations were run on the
GoeGrid cluster at the University of Göttingen, which is supported by the DFG
(grant INST 186/1353-1 FUGG) and MWK Niedersachsen (grant no. 45-10-19-F-02).
This project has received funding from the European Union’s Horizon 2020
research and innovation programme under Marie Skłodowska-Curie grant agreement
no. 893128.

## References

  * Marchetti  _et al._ (2013) M. C. Marchetti, J. F. Joanny, S. Ramaswamy, T. B. Liverpool, J. Prost, M. Rao, and R. A. Simha, Rev. Mod. Phys. 85, 1143 (2013).
  * Bechinger  _et al._ (2016) C. Bechinger, R. Di Leonardo, H. Löwen, C. Reichhardt, G. Volpe, and G. Volpe, Rev. Mod. Phys. 88, 045006 (2016).
  * Fily and Marchetti (2012) Y. Fily and M. C. Marchetti, Phys. Rev. Lett. 108, 235702 (2012).
  * Takatori and Brady (2015) S. C. Takatori and J. F. Brady, Phys. Rev. E 91, 032117 (2015).
  * Cates and Tailleur (2015) M. E. Cates and J. Tailleur, Annu. Rev. Condens. Matter Phys. 6, 219 (2015).
  * Levis  _et al._ (2017) D. Levis, J. Codina, and I. Pagonabarraga, Soft Matter 13, 8113 (2017).
  * Solon  _et al._ (2018) A. P. Solon, J. Stenhammar, M. E. Cates, Y. Kafri, and J. Tailleur, New J. Phys. 20, 075001 (2018).
  * Mora  _et al._ (2016) T. Mora, A. M. Walczak, L. Del Castello, F. Ginelli, S. Melillo, L. Parisi, M. Viale, A. Cavagna, and I. Giardina, Nature physics 12, 1153 (2016).
  * Ballerini  _et al._ (2008) M. Ballerini, N. Cabibbo, R. Candelier, A. Cavagna, E. Cisbani, I. Giardina, V. Lecomte, A. Orlandi, G. Parisi, A. Procaccini, M. Viale, and V. Zdravkovic, Proceedings of the National Academy of Sciences 105, 1232 (2008).
  * Mandal and Sollich (2021) R. Mandal and P. Sollich, Proc. Natl. Acad. Sci. 118, e2101964118 (2021).
  * Caprini  _et al._ (2020a) L. Caprini, U. Marini Bettolo Marconi, and A. Puglisi, Phys. Rev. Lett. 124, 078001 (2020a).
  * Caprini and Marini Bettolo Marconi (2021) L. Caprini and U. Marini Bettolo Marconi, Soft Matter 17, 4109 (2021).
  * Caprini  _et al._ (2020b) L. Caprini, U. M. B. Marconi, C. Maggi, M. Paoluzzi, and A. Puglisi, Phys. Rev. Research 2, 023321 (2020b).
  * Henkes  _et al._ (2020) S. Henkes, K. Kostanjevec, J. M. Collinson, R. Sknepnek, and E. Bertin, Nature Communications 11, 1 (2020).
  * Szamel and Flenner (2021) G. Szamel and E. Flenner, EPL (Europhysics Letters) 133, 60002 (2021).
  * Yllanes  _et al._ (2017) D. Yllanes, M. Leoni, and M. C. Marchetti, New Journal of Physics 19, 103026 (2017).
  * Das  _et al._ (2018) R. Das, M. Kumar, and S. Mishra, Phys. Rev. E 98, 060602 (2018).
  * Martinez  _et al._ (2018) R. Martinez, F. Alarcon, D. R. Rodriguez, J. L. Aragones, and C. Valeriani, The European Physical Journal E 41, 1 (2018).
  * Bera and Sood (2020) P. K. Bera and A. K. Sood, Phys. Rev. E 101, 052615 (2020).
  * Abbaspour and Klumpp (2021) L. Abbaspour and S. Klumpp, Phys. Rev. E 103, 052601 (2021).
  * Weeks  _et al._ (1971) J. D. Weeks, D. Chandler, and H. C. Andersen, The Journal of chemical physics 54, 5237 (1971).
  * Stürmer  _et al._ (2019) J. Stürmer, M. Seyrich, and H. Stark, The Journal of chemical physics 150, 214901 (2019).
  * Banerjee  _et al._ (2022) J. P. Banerjee, R. Mandal, D. S. Banerjee, S. Thutupalli, and M. Rao, Nature Communications 13, 4533 (2022).
  * Williams  _et al._ (2022) S. Williams, R. Jeanneret, I. Tuval, and M. Polin, nature communications 13, 4776 (2022).
  * Singh  _et al._ (2017) D. P. Singh, U. Choudhury, P. Fischer, and A. G. Mark, Advanced Materials 29, 1701328 (2017).
  * Mu  _et al._ (2022) Y. Mu, L. Lei, J. Zheng, W. Duan, Z. Wang, J. Tang, Y. Gao, and Y. Wang, ACS nano 16, 6801 (2022).
  * Patteson  _et al._ (2018) A. E. Patteson, A. Gopinath, and P. E. Arratia, Nature communications 9, 1 (2018).
  * Kumar  _et al._ (2014) N. Kumar, H. Soni, S. Ramaswamy, and A. Sood, Nature communications 5, 1 (2014).
  * Gupta  _et al._ (2022) R. K. Gupta, R. Kant, H. Soni, A. K. Sood, and S. Ramaswamy, Phys. Rev. E 105, 064602 (2022).
  * Higham (2001) D. J. Higham, SIAM review 43, 525 (2001).

Supplemental Materials:
Long-range Velocity Correlation from Active Dopants

##  I Hydrodynamic Theory

We follow the continuum approach used by Henkes et al. Henkes  _et al._ (2020)
but with a crucial difference in terms of describing the medium, which will be
explained in detail below. We start by writing the equation of motion of the
displacement field,

$\zeta\dot{\mathbf{u}}(\mathbf{r})=\mathbf{F}(\mathbf{r},t)+\nabla\cdot\hat{\sigma}$
(S1)

where $\mathbf{u}(\mathbf{r})$ is the displacement field of the whole medium
(where we do not distinguish whether the medium is made of active or passive
particles), $\zeta$ is the friction coefficient and $\hat{\sigma}$ is the
stress tensor. We also define
$\mathbf{F}(\mathbf{r},t)=\sum_{i}\mathbf{F}_{i}(\mathbf{r},t)$ where
$\mathbf{F}_{i}(\mathbf{r},t)$ is the force due to the $i$-th active particle
which is of the form,

$\mathbf{F}_{i}(\textbf{r},t)=a^{2}f\>\mathbf{n}_{i}(t)\>\delta(\mathbf{r}-\mathbf{r}_{i}(t))$
(S2)

where $a$ is a microscopic length scale typically given by the particle
parameter, $f$ is the magnitude of propulsion force acting on each active
particle, $\textbf{n}_{i}$ is the unit vector associated with the orientation
of the active force of the $i$-th active particle and $\textbf{r}_{i}$ is the
position of the $i$-th active particle.

We then write down the Fourier transform (both in space and time) of the
active forces as,

$\mathbf{F}_{i}(\mathbf{q},\omega)=a^{2}f\int\int\mathbf{n}_{i}(t)\>\delta(\mathbf{r}-\mathbf{r}_{i}(t))\>e^{i\omega
t}e^{i\mathbf{q}\cdot\mathbf{r}}dt\>\>d^{2}\mathbf{r}.$ (S3)

Next, considering a spatially static active particle configuration (described
by $\mathbf{r}_{i}(t)\equiv\mathbf{r}_{i}$) where the orientation of the
active particles can still evolve, we get

$\mathbf{F}_{i}(\mathbf{q},\omega)=a^{2}f\>e^{i\mathbf{q}\cdot\mathbf{r}_{i}}\int\mathbf{n}_{i}(t)\>e^{i\omega
t}\>dt\ .$ (S4)

We can now calculate the force correlation as,

$\mathbf{F}_{i}(\mathbf{q},\omega)\cdot\mathbf{F}_{j}(\mathbf{q^{\prime}},\omega^{\prime})=a^{4}{f}^{2}\>e^{i\mathbf{q}\cdot\mathbf{r}_{i}}\>e^{i\mathbf{q^{\prime}}\cdot\mathbf{r}_{j}}\int\int\mathbf{n}_{i}(t)\cdot\mathbf{n}_{j}(t^{\prime})e^{i\omega
t}\>e^{i\omega^{\prime}t^{\prime}}dt\>dt^{\prime}$ (S5)

and with an appropriate average over an ensemble of configurations of force
orientations we get,

$\displaystyle\langle\mathbf{F}_{i}(\mathbf{q},\omega)\cdot\mathbf{F}_{j}(\mathbf{q^{\prime}},\omega^{\prime})\rangle$
$\displaystyle=a^{4}{f}^{2}\>e^{i\mathbf{q}\cdot\mathbf{r}_{i}}\>e^{i\mathbf{q^{\prime}}\cdot\mathbf{r}_{j}}\int\int\langle\mathbf{n}_{i}(t)\cdot\mathbf{n}_{j}(t^{\prime})\rangle
e^{i\omega t}\>e^{i\omega^{\prime}t^{\prime}}dt\>dt^{\prime}$
$\displaystyle=a^{4}{f}^{2}\>e^{i\mathbf{q}\cdot\mathbf{r}_{i}}\>e^{i\mathbf{q^{\prime}}\cdot\mathbf{r}_{j}}\int\int
e^{-\frac{|t-t^{\prime}|}{\tau_{p}}}\>\delta_{ij}e^{i\omega
t}\>e^{i\omega^{\prime}t^{\prime}}dt\>dt^{\prime}$
$\displaystyle=a^{4}{f}^{2}\>\delta_{ij}\>e^{i\mathbf{q}\cdot\mathbf{r}_{i}}\>e^{i\mathbf{q^{\prime}}\cdot\mathbf{r}_{j}}\>2\pi\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime})$

Note that this correlation still depends on a specific configuration in terms
of the positions of the active particles $\\{\mathbf{r}_{i}\\}$. Summing over
the particles

$\displaystyle\sum_{ij}\langle\mathbf{F}_{i}(\mathbf{q},\omega)\cdot\mathbf{F}_{j}(\mathbf{q^{\prime}},\omega^{\prime})\rangle$
$\displaystyle=\sum_{ij}a^{4}{f}^{2}\>\delta_{ij}\>e^{i\mathbf{q}\cdot\mathbf{r}_{i}}\>e^{i\mathbf{q^{\prime}}\cdot\mathbf{r}_{j}}\>2\pi\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime})$
$\displaystyle=\sum_{i}a^{4}{f}^{2}\>e^{i(\mathbf{q}+\mathbf{q^{\prime}})\cdot\mathbf{r}_{i}}\>2\pi\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime}),$

and averaging over their positions gives the force correlation

$\displaystyle\langle\mathbf{F}(\mathbf{q},\omega)\cdot\mathbf{F}(\mathbf{q^{\prime}},\omega^{\prime})\rangle$
$\displaystyle=a^{4}\sum_{\mathbf{r}_{i}}P(\mathbf{r}_{i}){f}^{2}\>e^{i(\mathbf{q}+\mathbf{q^{\prime}})\cdot\mathbf{r}_{i}}\>2\pi\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime})$
$\displaystyle=N_{a}a^{4}\sum_{\mathbf{r}}P(\mathbf{r}){f}^{2}\>e^{i(\mathbf{q}+\mathbf{q^{\prime}})\cdot\mathbf{r}}\>2\pi\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime})$
$\displaystyle=\frac{N_{a}{f}^{2}a^{4}}{L^{2}}\int
d^{2}\mathbf{r}\>e^{i(\mathbf{q}+\mathbf{q^{\prime}})\cdot\mathbf{r}}\>2\pi\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime})$
$\displaystyle=\frac{N_{a}{f}^{2}a^{4}(2\pi)^{3}}{L^{2}}\delta(\mathbf{q}+\mathbf{q^{\prime}})\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime})$
$\displaystyle=\frac{N_{a}{f}^{2}(2\pi)^{3}}{N}\frac{a^{4}N}{L^{2}}\delta(\mathbf{q}+\mathbf{q^{\prime}})\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime})$
$\displaystyle=\frac{N_{a}}{N}a^{2}{f}^{2}(2\pi)^{3}\delta(\mathbf{q}+\mathbf{q^{\prime}})\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime})$

where $N_{a}$ is the number of active particles, $N$ is the total number of
particles and $L$ is the linear dimension of the system. This finally gives us
the force correlation in the following form

$\langle\mathbf{F}(\mathbf{q},\omega)\cdot\mathbf{F}(\mathbf{q^{\prime}},\omega^{\prime})\rangle=\frac{\phi_{a}}{\phi_{a}+\phi_{p}}a^{2}{f}^{2}(2\pi)^{3}\delta(\mathbf{q}+\mathbf{q^{\prime}})\>\frac{2\tau_{p}}{1+(\tau_{p}\omega)^{2}}\>\delta(\omega+\omega^{\prime})$
(S6)

where $\phi_{a},\phi_{p}$ are area fractions of active particles and passive
particles respectively. Apart from the factor
$\frac{\phi_{a}}{\phi_{a}+\phi_{p}}$, this equation has the same form as
derived in Ref. Henkes  _et al._ (2020) for a purely active system.

The force correlation being known, we can follow Ref. Henkes  _et al._ (2020)
to obtain the velocity correlations: taking the Fourier transform of the
equation of motion S1 and decomposing it into transversal and longitudinal
modes, one can write the velocity correlation (in Fourier space) in terms of
the force correlation (also in Fourier space) that we derived in Eq. S6.
Taking the inverse Fourier transform of the velocity correlation, we obtain

$\begin{split}C_{vv}(\mathbf{r})\approx\left(\frac{\phi_{a}}{\phi_{a}+\phi_{p}}\right)\left(\frac{a^{2}{f}^{2}}{4\pi\zeta^{2}}\sqrt{\frac{\pi}{2r}}\right)\left(\frac{1}{{\xi_{L}}^{3/2}}e^{-{r/\xi_{L}}}+\frac{1}{{\xi_{T}}^{3/2}}e^{-{r/\xi_{T}}}\right)\end{split}$
(S7)

where $r$ is the spatial distance. Here the correlation lengths are defined as
$\xi_{T}=\left(\frac{\mu\tau_{p}}{\zeta}\right)^{1/2}$ and
$\xi_{L}=\left(\frac{(B+\mu)\tau_{p}}{\zeta}\right)^{1/2}$ where $B$ and $\mu$
are the bulk and the shear modulus of the system Henkes  _et al._ (2020),
which parameterize the stress tensor in Eq. S1. If we make the further
assumption that $B\gg\mu$ (which is most accurate for area fractions close to
the jamming transition) this leads to the final expression for the velocity
correlation function

$C(\mathbf{r})\approx\left(\frac{\phi_{a}}{\phi_{a}+\phi_{p}}\right)\left(\frac{a^{2}{f}^{2}}{4\pi\zeta^{2}}\sqrt{\frac{\pi}{2r}}\right)\frac{1}{{\xi_{L}}^{3/2}}e^{-{r/\xi_{L}}}.$
(S8)

Eq. S8 shows the dependence on density ($\phi_{a}$, $\phi_{p}$) and activity
parameters ($f$, $\tau_{p}$) explicitly in a closed-form expression.

##  II Details of the Particle Based Model

The overdamped Brownian dynamics simulations were performed in two spatial
dimensions. The equations of motion in Eq. 1 were integrated with an explicit
Euler-Maruyama scheme Higham (2001) with a time step $\Delta t=7\times
10^{-4}$ in a square box with periodic boundary conditions. All the particles
in the system interact only through the steric interactions where the force is
defined as, $\mathbf{F}_{i}=-\nabla_{i}U_{tot}$. $U_{tot}$, the purely
repulsive interaction potential of WCA type Weeks  _et al._ (1971) is given by

$U_{tot}(r_{ij})=\begin{cases}4\epsilon\left[\left(\frac{\sigma}{r_{ij}}\right)^{12}-\left(\frac{\sigma}{r_{ij}}\right)^{6}\right]+\epsilon,&r_{ij}<2^{1/6}\sigma\\\
0,&\text{otherwise}\end{cases}$ (S9)

where $r_{ij}=|\mathbf{r}_{i}-\mathbf{r}_{j}|$ is the magnitude of separation
vector of the pair of particles $i$ and $j$ and $\mathbf{r}_{i}$,
$\mathbf{r}_{j}$ are the position vectors of the same pair. Here $\sigma$
stands for the parameter representing particle diameter and $\epsilon$ is the
energy scale of the interactions. By using $\sigma=1.0$ and $\epsilon=1.0$ we
set the scales of length and energy in our description to $\sigma$ and
$\epsilon$, respectively. We choose a fixed box size of $L=107\sigma$ and
various values of the area fraction $\phi_{\rm tot}$. The number of particles
($N$) in the system is then determined by $N=\frac{\phi_{\rm
tot}4L^{2}}{\pi\sigma^{2}}$. The number of active and passive particles
($N_{a}$ and $N_{b}$ respectively) are chosen using the relations
$N_{a}=\frac{\phi_{a}}{\phi_{\rm tot}}N$ and $N_{p}=\frac{\phi_{p}}{\phi_{\rm
tot}}N$, where $\phi_{a}$ and $\phi_{p}$ are the fraction of active and
passive particles that we choose for a particular simulation. To obtain a
reliable estimate of the correlation length, we performed 50 independent
simulations, each with a randomly chosen initial configuration. After each
simulation reached a steady state in terms of both energy and correlation
length, we obtained 50 snapshots of the system at regular intervals. In total,
we obtained 2500 snapshots by averaging over the 50 simulations and the 50
snapshots from each simulation. By averaging over a large number of snapshots
and simulations, we were able to obtain accurate estimates of the correlation
length.

##  III System size effect

We conducted a systematic investigation of finite-size effects for a system
with a fixed total density of $\phi_{\rm tot}=0.95$
($\phi_{a}=0.01,\phi_{p}=0.94$), while using the same values for persistence
time and self-propulsion force as in the plots in the main text. Specifically,
we explored three different system sizes ($L=53\sigma$, $106\sigma$,
$177\sigma$) and found that the results exhibit a considerable degree of
overlap across these sizes (Figure S1).

Figure S1: (a) Correlation length $\xi_{L}$ as a function of the magnitude of
the active force for two different persistence times ($\tau_{p}=10$ and $100$)
and three different system sizes (indicated by different colors and markers,
see legend). The results indicate that the correlation length is unaffected by
changes in the active force, for fixed values of persistence time and system
size. (b) Correlation length $\xi_{L}$ as a function of persistence time
$\tau_{p}$ for different values of the active force magnitude ($f=0.1$,
$0.25$, $0.5$, $1$) and three different system sizes (indicated by different
colors and markers). The results demonstrate that the correlation length
exhibits a power-law growth (with exponent $\frac{1}{2}$) with increasing
persistence time, and that this behavior is independent of the magnitude of
the active force and system size.

##  IV Movies

We have supplemented the information provided in the paper with a set of
movies.

  * •

Movie M1.mov displays a binary mixture of active and passive particles across
different ($\phi_{a}$, $\phi_{p}$) values over time. As above, $\phi_{a}$ and
$\phi_{p}$ denote the area fractions of active and passive particles,
respectively. Both $\phi_{a}$ and $\phi_{p}$ values for each snapshot can be
found in the corresponding row and column. The video showcases a significant
correlation in velocity orientation when the overall density of the active and
passive mixture is high enough. The self-propulsion force $f$ and persistence
time $\tau_{p}$ are kept constant across the different densities, with values
of 0.5 and 1000, respectively.

In the subsequent movies, we maintained a fixed area fraction of active
particles at 0.01 while adding passive particles with varying area fractions
to the system. The corresponding area fraction of passive particles for each
movie is given below. In the visualization on the left hand side, blue
particles represent active particles, while red particles represent passive
particles. On the right hand side, the same system is shown with a color code
that encodes the velocity orientation of the particles. In all of the movies,
the self-propulsion force $f$ is set to 0.5, and the persistence time
$\tau_{p}$ is 1000.

  * •

Movie M2.mov depicts a system with a passive particle area fraction of 0.15
and an active particle area fraction of 0.01 . In this case, the system
remains in a homogeneous state, with active particles freely navigating the
system. There is no apparent correlation in the velocity orientation of the
particles in the system.

  * •

Movie M3.mov displays a system with a passive particle area fraction of 0.5
and an active particle area fraction of 0.01. In this case, the excluded
volume interaction between the particles becomes significant. Small domains
with aligned velocity emerge in the system (on the right-hand side of the
movie, small correlated domains with the same color can be observed). Notably,
the active particles carve their own path through the passive medium by
plowing through the particles (on the left-hand side of the movie).

  * •

Movie M4.mov shows a system with a passive particle area fraction of 0.79 and
an active particle area fraction of 0.01. It is evident that the average size
of domains with aligned velocity increases, indicating that the active
particles start to push the passive particles to move together coherently.
This emphasizes the importance of a sufficiently large total area fraction of
particles, so that strong velocity correlations can be observed.

  * •

Movie M5.mov shows a system with a passive particle area fraction of 0.99 and
an active particle area fraction of 0.01. Remarkably, it shows that the entire
system becomes correlated, emphasizing the strong influence of active
particles even at very low area fractions, particularly in systems with high
total area fractions. Note that in this case, a fraction of only $1\%$ of
active particles can correlate the entire system. The crystalline structure of
the system is not essential here: similar results can be obtained in an
amorphous systems consisting of a mixture of particle sizes.