# A Dark Matter Probe in Accreting Pulsar-Black Hole Binaries

Ali Akil<EMAIL_ADDRESS>Department of Physics, Southern University of
Science and Technology, Shenzhen 518055, China  Department of Physics, The
Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Jockey Club Institute for Advanced
Study, The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Qianhang Ding
<EMAIL_ADDRESS>Department of Physics, The Hong Kong University of
Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China  Jockey Club Institute for Advanced
Study, The Hong Kong University of Science and Technology,
Clear Water Bay, Kowloon, Hong Kong, China

###### Abstract

The accretion of dark matter (DM) into astrophysical black holes slowly
increases their mass. The rate of this mass accretion depends on the DM model
and the model parameters. If this mass accretion effect can be measured
accurately enough, it is possible to rule out some DM models, and, with the
sufficient technology and the help of other DM constraints, possibly confirm
one model. We propose a DM probe based on accreting pulsar-black hole
binaries, which provide a high-precision measurement on binary orbital phase
shifts induced by DM accretion into black holes, and can help rule out DM
models and study the nature of DM.

##  I Introduction

Dark matter (DM) consists around $26\%$ of the energy density in the Universe
[1], which deeply influences the cosmic evolution and shapes the large scale
structure. With decades of observations and studies, the Lambda-cold dark
matter ($\Lambda$CDM) model has become the standard model of Big Bang
cosmology [2, 3, 4]. However, there still exist unsolved DM problems between
theoretical predications and observations, such as the core-cusp problem [5,
6, 7], missing satellites problem [8, 9], etc., which push us to study the
nature of DM.

The difficulty in studying the nature of DM is its weak or absent interaction
with baryonic matter, we are not able to observe it directly. This left space
for a large number of theoretical models accounting for the gravitational
effects of DM, like weakly interacting massive particles (WIMPs) [10, 11],
ultralight DM [12], primordial black holes (PBHs) [13], modified gravity [14],
etc. Some of these models have restricted parameter ranges due to
observational constraints [15, 16, 17]. However, most of the models are not
completely ruled out and are waiting to be examined by future DM probes.

To develop a DM probe to study the nature of DM, the accretion of DM into a
black hole (BH) is a potential channel, where various DM models contribute
different accretion rates and result in distinct BH mass increment after long
enough timescale. Such an accretion effect of DM has been well studied in a
number of scenarios for understanding DM properties [18, 19, 20, 21, 22],
however, due to the slow accretion rate of DM, a measurable DM accretion
effect needs a supermassive host astrophysical object, which owns a complex
matter surrounding environment and causes difficulty in extracting the DM
information. The introduction of pulsar in BH accreting DM changes the story.
The pulsar can emit stable pulse signals, which provide a high time-resolution
measurement on its surrounding environment [23]. This high precision
measurement can even measure the mass loss of a pulsar due to electromagnetic
radiation [24]. When a pulsar rotates around an accreting BH, the cumulative
DM accretion slowly increases the BH mass and influences orbital evolution of
pulsar-black hole (PSR-BH) binaries. This deviation from a standard general-
relativistic orbital evolution would produce an orbital phase shift, which can
be detected in pulse timing after a long-time accretion cumulation. Since
accretion effects in various DM models are different, the detected orbital
phase shift induced by DM accretion is corresponding with DM models and model
parameters, which can work as a DM probe in ruling out DM models and studying
the nature of DM.

Although, we do not observe any PSR-BH binaries so far, numerous studies have
estimated the number of PSR-BH binaries inside the Milky Way is around
$\mathcal{O}(10)-\mathcal{O}(1000)$ [25, 26]. With the improvement of
sensitivity in radio telescopes, such as Five-Hundred Metre Aperture Spherical
Radio Telescope (FAST) [27], MeerKAT [28] and future Square Kilometre Array
(SKA) [29], the PSR-BH binaries are expected to be observed in the near
future. Also, gravitational waves (GWs) provide another window for detection.
The low frequency GW detectors, like Laser Interferometer Space Antenna (LISA)
[30], have high sensitivity in detecting GWs from small mass ratio binaries,
which could be PSR-BH binaries. A joint observation by radio telescopes and GW
detectors can enhance the possibility of the detection of PSR-BH binaries.

This paper is organized as follows. In Sec. II, we give a brief introduction
on the DM accretion rate into an astrophysical BH, including the accretion
rate of WIMPs in Sec. II.1, the accretion rate of a hot particle DM model in
Sec. II.2, the accretion rate of ultralight DM in Sec. II.3, a brief
estimation on the accretion rate of PBHs in Sec. II.4, and a comment on
baryonic matter accretion about how and in what cases it affects our proposal
in Sec. II.5. Then in Sec. III, we discuss about the pulse timing in PSR-BH
binaries, including using PSR-BH binaries to study the new phenomena in Sec.
III.1 and studying BH mass changing effect in PSR-BH binaries in Sec. III.2.
Finally, in Sec. IV, we numerically calculate the mass accretion effect in
PSR-BH binaries and use it to constrain the DM models, including constraining
WIMPs in Sec. IV.1, ultralight DM in Sec. IV.2, PBHs in Sec. IV.3.

##  II Dark Matter Accretion

Accretion into stellar objects, particularly BHs, has been studied for a long
time. In 1952, the Bondi accretion formula was derived for a non-relativistic
gas cloud [31]. Then in the framework of General Relativity, Michel derived a
formula for hot gas accretion [32]. Unruh, on the other hand, worked out the
case for scalar fields by solving the Klein-Gordon equation in a Schwarzschild
BH background geometry [33]. The reader will notice that different DM models
obey different accretion formulas, which will lead to different accretion
rates. For example, WIMPs obey Bondi’s formula, whereas hot particle DM, being
relativistic, will obey Michel’s formula, etc. In each of the coming
subsections we will briefly introduce one accretion formula. All the accretion
rates basically scale like the square of the BH mass $M_{B}^{2}$, but the
constant which multiplies it varies heavily from one DM model to another. And
the accreting BH that we are considering is a Schwarzschild BH.

###  II.1 Weakly Interacting Massive Particles

WIMPs were until recently a long time most popular DM candidate. WIMPs are
beyond standard model particles that interact only very weakly except for
their gravitational field. WIMPs mass can be anywhere between
$2\,\mathrm{GeV}$ and $100\,\mathrm{TeV}$ [34]. The accretion of WIMPs by BHs
is captured by the well known Bondi accretion formula [31]. The Bondi
accretion assumes a spherically symmetric stellar object stationarily
accreting a cloud of non-relativistic matter particles. For practical reasons,
we will refer to the form in which it is presented in [35] (for a detailed
derivation, see Appendix. A), where the equation reads,

$\displaystyle\frac{\text{d}M_{\rm B}}{\text{d}t}=4\pi\lambda_{B}(GM_{\rm
B})^{2}\frac{\rho_{\infty}}{\gamma^{\frac{3}{2}}\,\Theta_{\infty}^{\frac{3}{2}}c^{3}}~{},$
(1)

with

$\displaystyle\Theta=\frac{k_{B}T}{mc^{2}}=\frac{c_{s}^{2}}{\gamma c^{2}}~{},$
(2)

being the dimensionless temperature. Here, the physical constants $G$, $c$,
$k_{B}$ are Newton’s constant, the speed of light, and the Boltzmann constant,
respectively. $\rho_{\infty}$ is the density at infinity. Infinity in this
context being what is relevant for such an astrophysical scale rather than for
example cosmological scale. $M_{\rm B}$ is the mass of the accreting BH, $m$
is the WIMPs mass and $\gamma$ is the so called polytropic constant,
characterizing a polytropic fluid of pressure $P$ and density $\rho$ for which
$P\sim\rho^{\gamma}$. Moreover,
$\lambda_{B}=\frac{1}{4}\left(\frac{2}{5-3\gamma}\right)^{\frac{5-3\gamma}{2(\gamma-1)}}$.
$c_{s}$ is the sound speed of DM 111An ideal CDM model is collisionless with a
zero sound speed, however, CDM particles still have a non-zero velocity
dispersion (see [36] for details), which causes a free streaming away from
gravitational collapse and produces an effective sound speed., and it is
constrained by the rotation curve of the Milky Way as $c_{s}<10^{-4}c$ [37],
which puts an upper bound on dimensionless temperature
$\Theta<\mathcal{O}(10^{-8})$. As one would expect, the colder (smaller
dispersion) the DM is, the higher the accretion rate. Moreover, the heavier it
is (within the range of validity of the Bondi accretion), the higher its
accretion rate.

###  II.2 Hot Particle Dark Matter

For the hot DM case, general relativistic effects are quite important,
therefore we deal with Lorentz covariant quantities and have a full general
relativistic treatment. This is captured by the Michel accretion formula [32],
which has the same form as Bondi’s, with a different
$\lambda_{M}\neq\lambda_{B}$, a smaller mass and higher temperature, causing a
much smaller accretion rate (for detailed derivation, see Appendix. B).

$\displaystyle\frac{\text{d}M_{\rm B}}{\text{d}t}=4\pi\lambda_{M}(GM_{\rm
B})^{2}\frac{\rho_{\infty}}{\gamma^{\frac{3}{2}}\,\Theta_{\infty}^{\frac{3}{2}}c^{3}}~{}.$
(3)

Where $\lambda_{M}$, unlike in the Bondi accretion, depends also on the sound
speed of the medium but in general varies between $1$ and $2$ [35]. Since this
accretion rate is extremely small and cannot be observed in the method we
propose, we will not elaborate on it. But it is good to note that this model
is one of those that will be favoured by the experiment that we propose in the
case of the measurement showing there is no growth in the astrophysical BH
mass.

###  II.3 Ultralight Dark Matter

Recently, the scalar field models of DM have arguably become the most popular.
The ultralight DM model, on the small mass between $10^{-24}\,{\rm eV}$ and
$1\,{\rm eV}$, is believed to solve the small scale problems of CDM [12].

The accretion rate of the ultralight DM is derived from the Klein-Gordon
equation on a non-rotating BH background [33]. The mass accretion rate by a
non-rotating BH of mass $M_{\mathrm{B}}$ traveling with velocity $v$ through
uniformly distributed ultralight DM of mass $m_{\mathrm{ul}}$ and density
$\rho_{\mathrm{DM}}$ can be expressed as follows,

$\displaystyle\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}=\frac{32\pi^{2}(GM_{\mathrm{B}})^{3}m_{\mathrm{ul}}\rho_{\mathrm{DM}}}{\hbar
c^{3}v[1-\exp(-\xi)]}~{},$ (4)

where $\xi$ is defined as $\xi\equiv 2\pi GM_{\mathrm{B}}m_{\mathrm{ul}}/\hbar
v$ and $\hbar$ is the reduced Planck’s constant.

###  II.4 Primordial Black Holes

Apart from ultralight DM with mass smaller than $1\,\mathrm{eV}$ and WIMPs
with mass between $2\,\mathrm{GeV}$ and $100\,\mathrm{TeV}$, some ultraheavy
objects can also be the DM candidate, in particular, primordial black holes
(PBHs). PBHs were first introduced in [13], which proposed that primordial
perturbations could collapse to a BH, and PBHs was later used in accounting
for massive astrophysical compact halo objects (MACHOs) [38, 39].

As PBHs are much heavier than the other DM models that we have considered
above, the way they can be accreted by astrophysical BHs differs radically.
Consider the mass of PBH is above $1\,M_{\odot}$, the accretion of a PBH in an
astrophysical BH behaves like a binary evolution, and its accretion rate is
related with the merger rate of this binary system, which can be estimated
with the mean free path of the astrophysical BH and its moving velocity.

Given the amount of DM in our galaxy, if we distribute it into PBHs rather
than the other types of DM, those would be highly separated from each other,
due to their large mass. Here, we assume the PBH mass function is
monochromatic. Then, the mean free path $l$ in this system depends on the
cross section of BHs $\sigma\simeq 27\pi(GM_{\mathrm{B}}/c^{2})^{2}$ [40] and
number density of PBHs $n=\rho_{\mathrm{DM}}/M_{\mathrm{PBH}}$, which can be
expressed as follows,

$\displaystyle l_{f}=\frac{1}{\sigma
n}=\frac{1}{27\pi}\left(\frac{c^{2}}{GM_{\mathrm{B}}}\right)^{2}\frac{M_{\mathrm{PBH}}}{\rho_{\mathrm{DM}}}~{}.$
(5)

Then, the mean free time for this astrophysical BH can be approximated as
$t_{f}\sim l_{f}/v$, where $v$ is the relative velocity between the
astrophysical BH and PBHs. The PBHs accretion rate can be evaluated as
follows,

$\displaystyle\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}\simeq\frac{M_{\mathrm{PBH}}}{t_{f}}\simeq
27\pi(GM_{\mathrm{B}})^{2}\frac{\rho_{\mathrm{DM}}v}{c^{4}}~{}.$ (6)

###  II.5 Comment on Baryonic Matter Accretion

One important thing to take into account is the accretion of the baryonic
matter. A detailed and careful account is certainly needed. However, as we
keep the precise computation for a future project, we try here discuss the
magnitude of that accretion. First, one should notice that if the PSR-BH
binary is found outside the galactic disk, the baryonic matter there is very
scarce and baryonic accretion is to be comfortably ignored. On the other hand,
if the binary is in the galactic disk, we might expect a significant baryonic
matter accretion. However, as we are comparing the behaviour of the binary for
different DM models, and that baryonic matter accretion is the same in both
cases, this will decrease its importance in the process. It will contribute as
an added term to the BH mass, in both the cases that we will be comparing,
which decreases its significance. And finally, an important point is that in
our galaxy, the baryonic matter is in its great majority in stars (around
$90\%$, see [41] page 2). Therefore, we assume that if a whole star is
swallowed by a Milky Way BH we will be able to see that.

##  III Pulsar-Black Hole Binaries

###  III.1 Background

The PSR-BH binary system is the holy grail in radio astronomy, since the
stable pulse signals emitted from pulsar can provide high precision
measurements on the strong gravity field around the BH [42]. Due to the motion
of the pulsar around the BH, the pulse signals could be affected by the
periodic motion of the pulsar and the gravitational field in the binary, and
such effects contribute three different kinds of time delay in measuring the
orbital Time-of-Arrival (TOA) of receiving pulse signals [43] as follows,

$\displaystyle\Delta_{\mathrm{orb}}\mathrm{TOA}=\Delta_{R}+\Delta_{E}+\Delta_{S}~{}.$
(7)

Here, $\Delta_{R}$ is the Rømer delay, which describes the elapsed time of
light in passing the binary orbit. $\Delta_{E}$ is the Einstein delay that is
the general-relativistic time dilation in the PSR-BH gravitational field.
$\Delta_{S}$ is the Shapiro delay, which is the delayed time of light in
curved spacetime.

After measuring TOA of receiving pulses, we can use them to fit a given model
by minimizing the timing residual between data and model predictions (see
TEMPO2 program [44] for details). It gives five Post-Keplerian (PK)
parameters, which can help determine the mass of the binary components and the
orbital parameters of the binary. With these PK parameters, new phenomena can
be tested in a given binary model. The new phenomena produces a different
gravitational effect, which could influence the orbit evolution in a binary,
and such an orbit deviation would produce a significant orbital phase shift,
which could be detected after long enough observation time. In the general
relativistic background, the orbital phase shift $\Delta\phi$ can be
calculated as follows,

$\displaystyle\Delta\phi(t)=2\pi\int_{0}^{t}f(\tau)d\tau-2\pi\int_{0}^{t}f_{\mathrm{GR}}(\tau)d\tau~{},$
(8)

where $f_{\mathrm{GR}}$ is the orbital frequency in general relativity without
new phenomena and $f$ is the orbital frequency with new phenomena.

In order to ensure the detection of orbital phase shift induced by the new
phenomena, the measurement uncertainty of the orbital phase shift
$\sigma_{\Delta\phi}$ should be smaller than the measured $\Delta\phi$. The
uncertainty of orbital phase shift is determined by the single measurement
error and the number of independent measurements. Assuming the observation
time per day $t_{\mathrm{obs}}\simeq 10\,\mathrm{hrs}$, the orbital phase
error in a single continuous measurement is decided by the orbital phase
measurement error within one orbital period $\sigma_{\phi}$ divided by the
number of orbital periods within one continuous observation
$N_{P_{\mathrm{b}}}$, which is $\sigma_{\phi}/N_{P_{\mathrm{b}}}$. Here
$\sigma_{\phi}$ can be maximally estimated as a pulse period $P$ divided by
orbital period $P_{\mathrm{b}}$ and $N_{P_{\mathrm{b}}}$ equals
$t_{\mathrm{obs}}/P_{\mathrm{b}}$. In calculating the number of independent
measurements, we assume the observation is under running once per day, the
number of independent measurements within observation time $t$ is
$t/1\,\mathrm{days}$. Then the uncertainty of the orbital phase shift can be
calculated as follows (also see [45, 46]),

$\displaystyle\sigma_{\Delta\phi}=\frac{2\pi}{\sqrt{t/1\,\mathrm{day}}}\frac{P}{t_{\mathrm{obs}}}~{}.$
(9)

The detection of new phenomena in PSR-BH binaries requires two conditions, one
is that the orbital phase shift within observation time should be larger than
its measurement uncertainty and the other one is that the observation time for
the orbital phase shift cannot be longer than the duty time of the radio
telescope $T_{\mathrm{duty}}$ and the merger time of the binary
$T_{\mathrm{merger}}$. The conditions can be expressed as follows,

$\displaystyle|\Delta\phi(t)|>\sigma_{\Delta\phi}(t)~{},~{}~{}~{}t<\min(T_{\mathrm{duty}},T_{\mathrm{merger}})~{}.$
(10)

###  III.2 Mass Changing Binaries

The physical background of the PSR-BH binary is very complex, it includes the
baryonic matter and DM. Such a matter surrounding environment could cause
matter accretion effect around the BH, which slowly increases the mass of the
BH (the accretion effect around the pulsar is neglected, due to its relatively
small mass, the ratio of accretion rate between pulsar and BH would be smaller
than $\mathcal{O}(1\%)$, if the BH mass is larger than $10\,M_{\odot}$).
Although this matter accretion effect is extremely weak, it can still induce a
detectable orbital phase shift in PSR-BH binaries after a long enough
observation time. This orbital phase shift depends on the surrounding matter
density and matter properties, especially, DM properties. Various DM models
and model parameters would produce different orbital phase shifts in the PSR-
BH binary, which can be used to distinguish the DM models and constrain the
parameter regions of DM within a model.

In this scenario, studying the orbital evolution of the PSR-BH binary with a
changing BH mass is essential. During the orbital evolution of the PSR-BH
binary, the gravitational radiation takes the gravitational energy and angular
momentum away from the binary system. Follow [47] which is calculated as
follows,

$\displaystyle
P=\frac{G}{5c^{5}}\left(\frac{\text{d}^{3}Q_{ij}}{\text{d}t^{3}}\frac{\text{d}^{3}Q_{ij}}{\text{d}t^{3}}-\frac{1}{3}\frac{\text{d}^{3}Q_{ii}}{\text{d}t^{3}}\frac{\text{d}^{3}Q_{jj}}{\text{d}t^{3}}\right)~{},$
$\displaystyle\frac{\text{d}L_{i}}{\text{d}t}=-\frac{2G}{5c^{2}}\epsilon_{ijk}\frac{\text{d}^{2}Q_{mj}}{\text{d}t^{2}}\frac{\text{d}^{3}Q_{mk}}{\text{d}t^{3}}~{}.$
(11)

Here, $G$ is the Newton’s constant, $c$ is the speed of light,
$\epsilon_{ijk}$ is the three-dimensional Levi-Civita symbol, $Q_{ij}$ is a
tensor which is defined as $Q_{ij}=\sum_{\alpha}m_{\alpha}x_{\alpha
i}x_{\alpha j}$, its form in a binary system can be expressed as $Q_{xx}=\mu
d^{2}\cos^{2}\phi$, $Q_{yy}=\mu d^{2}\sin^{2}\phi$ and $Q_{xy}=Q_{yx}=\mu
d^{2}\sin\phi\cos\phi$, where $\mu$ is the reduced mass
$m_{1}m_{2}/(m_{1}+m_{2})$, $d$ is the distance between components in the
binary and $\phi$ is the orbital phase of the binary.

In Eq. (III.2), the magnitude of the gravitational radiation is determined by
the time derivative of $Q_{ij}$, where a changing BH mass contributes. Such a
changing BH mass not only contributes to the gravitational radiation, but also
influences the gravitational potential energy of the binary. Due to the
accretion of DM into BHs, the gravitational potential energy is transferred
from the DM to the PSR-BH binary, an overall contribution of gravitational
potential energy can be estimated as follows, a detailed derivation can be
found in Appendix. D.

$\displaystyle\frac{\text{d}E_{p}}{\text{d}t}=-\frac{Gm_{\mathrm{p}}}{a}\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}~{},$
(12)

where $a$ is the semi-major axis of the binary, $m_{\mathrm{p}}$ is the mass
of the pulsar and $M_{\mathrm{B}}$ is the BH mass. The total energy $E$ and
angular momentum $L$ of the PSR-BH binary are related with their orbital
parameters, semi-major axis $a$ and eccentricity $e$ as follows [47],

$\displaystyle
a=-\frac{Gm_{\mathrm{p}}M_{\mathrm{B}}}{2E}~{},~{}~{}~{}L^{2}=\frac{Gm_{\mathrm{p}}^{2}M_{\mathrm{B}}^{2}}{m_{\mathrm{p}}+M_{\mathrm{B}}}a(1-e^{2})~{}.$
(13)

Combining Eqs. (III.2–13), the time derivative of orbital parameters $da/dt$
and $de/dt$ can be numerically solved from the energy conservation and angular
momentum conservation, which can be expressed as follows,

$\displaystyle\frac{\text{d}E}{\text{d}t}=-P_{\mathrm{acc}}+\frac{\text{d}E_{p}}{\text{d}t}~{},~{}~{}~{}\frac{\text{d}L}{\text{d}t}=\frac{\text{d}L_{\mathrm{acc}}}{\text{d}t}~{}.$
(14)

Here, $P_{\mathrm{acc}}$ is the power of the gravitational radiation emitted
from the PSR-BH binary with a DM accretion into the BH.
$\text{d}L_{\mathrm{acc}}/\text{d}t$ is the time derivative of angular
momentum with the DM accretion effect. After obtaining the time derivative of
the semi-major axis $\text{d}a/\text{d}t$, the time derivative of orbital
frequency $\text{d}f/\text{d}t$ can be calculated from the Kepler’s third law
as follows,

$\displaystyle\frac{\text{d}f}{\text{d}t}=\frac{1}{4\pi}\frac{a^{-5/2}G^{1/2}}{(m_{\mathrm{p}}+M_{\mathrm{B}})^{1/2}}\left(a\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}-3(m_{\mathrm{p}}+M_{\mathrm{B}})\frac{\text{d}a}{\text{d}t}\right)~{}.$
(15)

Meanwhile, the orbital frequency without the DM accretion follows the standard
general relativistic evolution, which gives its time derivative as follows,

$\displaystyle\frac{\text{d}f_{\mathrm{GR}}}{\text{d}t}=-\frac{3}{4\pi}\frac{G^{1/2}(m_{\mathrm{p}}+M_{\mathrm{B}})^{1/2}}{a_{\mathrm{GR}}^{5/2}}\frac{\text{d}a_{\mathrm{GR}}}{\text{d}t}~{}.$
(16)

Here, subscript GR denotes that the evolution of physical quantity follows the
standard general relativitic evolution in [47]. Then the corresponding
frequency evolution $f(t)$ and $f_{\mathrm{GR}}(t)$ can be numerically solved,
following Eq. (8), orbital phase shift $\Delta\phi$ induced by the DM
accretion on the BH can be obtained. In order to make sure a detection of DM
accretion into the BH, Eq. (10) needs to be satisfied.

##  IV Dark Matter Accretion in PSR-BH Binaries

###  IV.1 Weakly Interacting Massive Particles

For WIMPs, we use Eq. (1) to calculate the effect of its accretion on PSR-BH
binaries. It shows that the DM accretion rate into BHs depends on the DM
density, BH mass and dimensionless temperature. Since the measurable pulsar
systems live in the Milky Way, the DM density follows a galactic DM density
profile, which we use the Navarro-Frenk-White (NFW) model to describe [48,
49],

$\displaystyle\rho(r)=\frac{\rho_{0}}{\frac{r}{r_{0}}(1+\frac{r}{r_{0}})^{2}}~{},$
(17)

where $\rho_{0}$ is the characteristic density and $r_{0}$ is the scale
length. By fitting the rotation curve of the Milky Way, we have the best fit
parameters are $\rho_{0}=0.052\,M_{\odot}/\mathrm{pc}^{3}$,
$r_{0}=8.1\,\mathrm{kpc}$ [50].

In a practical calculation, we take the position of the PSR-BH binary to be at
$10\,\mathrm{kpc}$ away from the center of the Milky Way. Following the
process introduced in Sec. III, the evolution of the orbital phase shift can
be obtained, as shown in Fig. 1.

Figure 1:  The evolution of orbital phase shift in the PSR-BH binary. We
assume that the mass of BH is $100\,M_{\odot}$, the mass of pulsar is
$1.6\,M_{\odot}$, the initial detected GW frequency is $0.01\,\mathrm{Hz}$ and
the eccentricity is $e=0$ ($e=0.6$) for solid curves (dashed curves). The gray
shadow regions are detectable orbital phase shift ranges by pulsar with pulse
period $1\,\mathrm{ms}$, $10\,\mathrm{ms}$ and $100\,\mathrm{ms}$,
respectively. The different colors of curves denote dimensionless temperature
of WIMPs are $10^{-11}$, $10^{-12}$ and $10^{-13}$, respectively.

We can find that a detectable timescale for WIMPs ($\Theta\sim 10^{-12}$)
accretion on $100\,M_{\odot}$ BHs in the PSR-BH binary system is
$\mathcal{O}(10)\,\mathrm{years}$. A smaller $\Theta$ represent a lower
temperature of WIMPs or heavier WIMPs mass, which could effectively increase
the accretion rate in Eq. (1) and enlarge the orbital phase shift in PSR-BH
binaries. Moreover, a small pulse period would effectively decrease the
detection timescale, due to the small measurement uncertainty in Eq. (9),
which indicates that millisecond PSR-BH binaries could be a potential
candidate for studying the nature of dark matter. Also, millisecond pulsars
are expected to form binaries with BHs via the “recycle” process [51].

In order to obtain a detectable parameter region for the BH mass and
dimensionless temperature, we follow the Eq. (10) in constraining parameters,
where we assume the pulsar mass is $1.6\,M_{\odot}$, pulse period is
$1\,\mathrm{ms}$ and the duty time of radio telescope is set as
$10\,\mathrm{years}$. The result is shown in Fig. 2.

Figure 2:  The detectable regions for dimensionless temperature $\Theta$ with
different BH masses in PSR-BH binaries. We assume the pulsar mass is
$1.6\,M_{\odot}$, pulse period is $1\,\mathrm{ms}$ and the duty time of radio
telescope is $10\,\mathrm{years}$. The blue (red) shadow regions represent the
result from PSR-BH binaries with GW frequency $10^{-2}\,\mathrm{Hz}$
($10^{-3}\,\mathrm{Hz}$). The solid (dashed) regions represent the result from
PSR-BH binaries with eccentricity $e=0$ ($e=0.6$).

We can find that a larger BH mass and higher orbital frequency could extend
the range of a detectable dimensionless temperature. This is because that a
larger BH mass and higher orbital frequency can increase the power of GW
radiation, which speeds up the shrinkage of orbit of the PSR-BH binary. Then,
orbital phase shift induced by the WIMPs accretion on BH would be enlarged.
Also, for a larger eccentricity (dashed curves in Fig. 2), stronger GW is
emitted which increases the orbital phase shift and extends the observable
$\Theta$ range. An optimal detectable range of dark matter properties depends
on various factors, namely, BH mass, orbital period, eccentricity, it means
the most advantageous binary parameters for studying the nature of dark matter
are different for various BH masses. The lower bound of shadow regions are
given by the condition
$\Delta\phi(T_{\mathrm{duty}})=\sigma_{\Delta\phi}(T_{\mathrm{duty}})$, while
the upper bound are given by the condition
$\Delta\phi(T_{\mathrm{merger}})=\sigma_{\Delta\phi}(T_{\mathrm{merger}})$.
Therefore, we can find that even though a large BH mass could improve the
detectable range of $\Theta$, its short merger time prohibits a large
cumulative orbital phase and decreases the detectable range of $\Theta$.

However, in usual recycle processes of millisecond pulsars, the mass of BH
companion can hardly exceed $100\,M_{\odot}$ [51, 26] and such a stellar mass
BH companion can only help constrain a small part of DM parameter range as
shown in Fig. 2. In order to probe a larger DM parameter range, a massive BH
companion (especially around $10^{4}\,M_{\odot}$ in Fig. 2) is needed in the
PSR-BH binary, which could be formed via several channels. The first channel
is having a millisecond PSR-BH binary in a dense stellar environment, then
through sequences of stellar interactions, where the BH encounters and merges
with a number of stars, and turns into an intermediate-mass black hole (IMBH).
Finally this millisecond PSR-BH binary becomes a millisecond PSR-IMBH binary,
[52, 53]. The second channel is having an IMBH formed near the center of a
Milky Way-like galaxy or globular clusters, then the strong gravitational
field of IMBH can help capture the single millisecond pulsar or millisecond
pulsar binaries. As a large population of pulsars, including millisecond
pulsar, is expected in globular clusters, such processes could finally lead to
the formation of millisecond PSR-IMBH binaries [54, 55]. The third possible
alternative channel is through the so-called hierarchical triple systems. The
hierarchical triple systems consist of a binary system (a pulsar and a low-
mass companion) orbiting an IMBH. The tidal interactions between the inner
IMBH and outer binary can disrupt the outer binary and kick out the low-mass
companion, forming a tight millisecond PSR-IMBH binary [56].

In this numerical calculation on the orbital phase shift, we use Eq. (15) and
(16), which is the orbital evolution in non-relativistic limit. This limit is
only valid in the inspiral phase of the PSR-BH binary, so an upper bound of
the orbital frequency during inspiral phase should be put in calculating the
orbital phase shift. This maximal inspiral phase frequency can be estimated as
$f_{\mathrm{insp}}=(a\eta^{2}+b\eta+c)/2\pi GM$ [57, 58, 59]. $\eta$ is the
symmetric mass ratio, which is defined as $\eta\equiv
m_{\mathrm{p}}M_{\mathrm{B}}/M^{2}$ in the PSR-BH binary, $M\equiv
m_{\mathrm{p}}+M_{\mathrm{B}}$ and the coefficients $a=0.29740$, $b=0.04481$,
$c=0.09556$ (see Table 1 of [60]).

For a larger detectable $\Theta$ range, we can extend the location of the PSR-
BH binary to a position with the higher DM density, where the mass accretion
rate is effectively enhanced, hence increases the orbital phase shift of the
binary and enlarges the detectable window of $\Theta$. We consider the PSR-BH
binary with $0.01\,\mathrm{Hz}$ GW frequency, circular orbit and
$1\,\mathrm{ms}$ pulse period. Within a $10\,\mathrm{years}$ observation, the
detectable range for $\Theta$ is shown in Fig. 3.

Figure 3:  The detectable regions for dimensionless temperature $\Theta$ with
different locations of PSR-BH binaries. We assume the pulsar mass is
$1.6\,M_{\odot}$. The red (blue) shadow regions represent the result from the
mass of BH is $100\,M_{\odot}$ ($1000\,M_{\odot}$).

It clearly shows a high DM density location near the center of the Milky Way
can extend the detectable range of $\Theta$, and the accretion effect on a
$1000\,M_{\odot}$ BH can extend the range of $\Theta$ up to a value close to
$10^{-8}$, which is an upper bound of $\Theta$ inside the Milky Way [37].

###  IV.2 Ultralight Dark Matter

To estimate the accretion rate of ultralight DM, we use Eq. (4), which gives
the mass accretion rate by a non-rotating BH of mass $M_{\mathrm{B}}$
traveling through a uniform distributed scalar field. Due to the accretion
rate of ultralight DM is relatively weak compared with WIMPs, we mainly focus
its accretion at the center of the Milky Way. Follow [61], we apply the
central density and virial velocity of the soliton in Eq. (4), which gives its
accretion rate,

$\displaystyle\frac{\text{d}M_{\rm
B}}{\text{d}t}=\frac{2.5\,M_{\odot}}{10^{17}\,{\rm yr}}\left(\frac{M_{\rm
B}}{\tilde{M}_{\mathrm{B}}}\right)^{2}\left(\frac{m_{\rm
ul}}{\tilde{m}_{\mathrm{ul}}}\right)^{6}\left(\frac{M_{\rm
sol}}{\tilde{M}_{\mathrm{sol}}}\right)^{4}~{},$ (18)

where the reference BH mass $\tilde{M}_{\mathrm{B}}=100\,M_{\odot}$, the
reference ultralight DM mass $\tilde{m}_{\mathrm{ul}}=10^{-22}\,\mathrm{eV}$
and the reference soliton mass $\tilde{M}_{\mathrm{sol}}=10^{10}\,M_{\odot}$.
The accretion rate of ultralight DM $\dot{M}_{\mathrm{B}}\propto
M_{\mathrm{B}}^{2}$, which is similar with the WIMPs’ accretion rate,
therefore, in a given parameter setting
$(M_{\mathrm{B}},m_{\mathrm{ul}},M_{\mathrm{sol}})$, the time evolution of
orbital phase shift is similar like the behavior in Fig. 1.

In order to find a detectable mass range of ultralight DM in PSR-BH binaries,
we set the soliton mass of the Milky Way is $10^{9}\,M_{\odot}$ [62], and we
follow the procedures introduced in Sec. III, using Eq. (10), where we assume
the pulsar mass is $1.6\,M_{\odot}$, pulse period is $1\,\mathrm{ms}$ and the
duty time of radio telescope is set as $10\,\mathrm{years}$. The result is
shown in Fig. 4.

Figure 4:  The detectable regions for the mass of ultralight DM
$m_{\mathrm{ul}}$ with different BH masses in PSR-BH binaries. The blue (red)
shadow regions represent the result from PSR-BH binaries with GW frequency
$10^{-2}\,\mathrm{Hz}$ ($10^{-3}\,\mathrm{Hz}$). The solid (dashed) regions
represent the result from PSR-BH binaries with eccentricity $e=0$ ($e=0.6$).

We can find that a larger BH mass and higher orbital frequency can extend the
detectable mass range of ultralight DM as previous discussions in WIMPs. And
the orbital phase shift induced by a larger ultralight DM mass could be
detected in PSR-BH binaries with a wider BH mass range, this is because a
larger ultralight DM mass can effectively enhance the accretion rate of
ultralight DM in Eq. (18). Also, contour in Fig. 4 is similar with that in
Fig. 2, the difference is only a rescale on the values of x-axis. The reason
is that the accretion rate in WIMPs and ultralight DM can be generalized as
$\dot{M}_{\mathrm{B}}=f(\alpha)M_{\mathrm{B}}^{2}$, where $\alpha$ is the
parameter of DM models, namely, $\Theta$ in WIMPs and $m_{\mathrm{ul}}$ in
ultralight DM. For each value of $f(\alpha)$, it relates with two different
values in $\Theta$ and $m_{\mathrm{ul}}$ and corresponds with a specific
accretion rate. Meanwhile, this accretion rate determines a BH mass range in
PSR-BH binaries, which is the same BH mass range in two DM models, it causes a
rescaled x-axis in their contours.

###  IV.3 Primordial Black Holes

In order to estimate the magnitude of PBHs accretion rate, we can compare PBHs
accretion rate in Eq. (6) with WIMPs accretion rate in Eq. (1). The parameters
$\lambda_{B}$, $\gamma$ are around $\mathcal{O}(1)$ in Eq. (1), then we have

$\displaystyle\frac{\dot{M}_{\mathrm{B}}^{\mathrm{P}}}{\dot{M}_{\mathrm{B}}^{\mathrm{W}}}\simeq\frac{27}{4}\frac{v}{c}\Theta^{3/2}~{}.$
(19)

Here, we use the velocity in rotation curve of the Milky Way [63] to
approximate the relative velocity between the BH and PBHs $v$, which is around
$200\,\mathrm{km/s}$. And we set $\Theta\sim 10^{-10}$, which is the largest
detectable value in a $10^{-2}\,\mathrm{Hz}$ PSR-BH binary. Then, the
accretion rate ratio between PBHs model and WIMPs can be evaluated as
$\dot{M}_{\mathrm{B}}^{\mathrm{P}}/\dot{M}_{\mathrm{B}}^{\mathrm{W}}\sim\mathcal{O}(10^{-17})$.
Such a small PBHs accretion rate can hardly be detected in our PSR-BH binary
systems.

Actually, the accretion rate in Eq. (6) is a time averaged accretion rate, and
the mass increment due to accreting PBHs only occur when the astrophysical BH
interacts with a PBH. This interaction happens within a very short time
interval, so a physical PBH accretion rate can be approximated as a summation
of Dirac delta functions

$\displaystyle\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}=M_{\mathrm{PBH}}\sum_{n=1}^{\infty}\delta(t-nt_{f})~{},$
(20)

where $t_{f}$ is the mean free time of each merger of the accreting BH and
PBH. Therefore, this mean free time cannot be too long that such a merger
event would not happen once during the observational duty time. We can give a
brief estimation on this mean free time by dividing the mean free path by the
relative velocity between the accreting BH and PBHs. In Eq. (5), we assume the
mass of accreting BH is $100\,M_{\odot}$, PBH mass is $1\,M_{\odot}$ and
$\rho_{\mathrm{DM}}=0.013\,M_{\odot}/\mathrm{pc}^{3}$, which is the DM density
at the location of the sun in NFW model [50]. As above, the relative velocity
between the accreting BH and PBHs is approximated by the velocity in rotation
curve $v\sim 200\,\mathrm{km/s}$. Then a mean free time is
$t_{f}=l_{f}/v\sim\mathcal{O}(10^{26})\,\mathrm{years}$, which is too long to
be detected in a PSR-BH binary.

##  V Conclusions

To summarize, we propose a DM probe that can be used to study the nature of
DM. This probe is based on the DM accretion in PSR-BH binaries. The DM
accretion could slowly increase the BH mass, and such a mass increment in a
PSR-BH binary would affect its gravitational radiation and decrease
gravitational potential energy, then induce an orbital phase shift in the
orbital evolution of the PSR-BH binary. Various DM models contribute different
mass accretion rates and induce distinct orbital phase shifts. After a long
observation time in pulse timing, such an orbital phase shift could be
detected. From observable orbital phase shifts, it can help us distinguish the
DM models and constrain the parameters in DM models.

In this work, we mainly focus on three DM models, WIMPs, ultralight DM, and
PBHs. For WIMPs, the accretion rate follows $\dot{M}_{\mathrm{B}}\sim
M_{\mathrm{B}}^{2}\rho_{\mathrm{DM}}\Theta^{-3/2}$, a larger BH mass, the
higher DM density and a lower dimensionless temperature, all can enhance the
WIMPs accretion rate. Consider WIMPs with $\Theta\sim\mathcal{O}(10^{-12})$,
accreted by a PSR-BH binary with a $100\,M_{\odot}$ BH mass at the position
$10\,\mathrm{kpc}$ away from the center of the Milky Way, the detectable
timescale of induced orbital phase shift is $\mathcal{O}(10)\,\mathrm{years}$.
Within a $10\,\mathrm{years}$ observation, a value of $\Theta$ can be detected
up to $10^{-8}$ in a PSR-BH binary with $10^{-2}\,\mathrm{Hz}$ GW frequency,
and the observable range of $\Theta$ can be extended with a larger BH mass and
higher orbital frequency. For ultralight DM, we mainly consider the accretion
inside the soliton of the Milky Way, due to their weak accretion rate, which
follows $\dot{M}_{\mathrm{B}}\sim M_{\mathrm{B}}^{2}m_{\mathrm{ul}}^{6}$, so a
larger mass of ultralight DM can effectively improve the detectability of
accretion effect. In our parameter setting, the orbital phase shift induced by
the accretion of ultralight DM with mass above
$\mathcal{O}(10^{-20})\,\mathrm{eV}$ can be detected, which could help
constrain the mass the ultralight DM. For PBHs, their number density inside
the Milky Way is small due to a large PBH mass, and this small number density
causes an extremely long mean free time of a BH encountering PBHs, around
$\mathcal{O}(10^{26})\,\mathrm{years}$, which is not a detectable timescale.
Therefore, a null detection result may indicate a possibility of PBHs as DM or
some other DM models with undetected parameter regions.

In above discussions, the mass accretion rate in different DM models basically
follows $\dot{M}_{\mathrm{B}}=f(\alpha)M_{\mathrm{B}}^{2}$, where $\alpha$ is
the DM parameter. Therefore, a detected DM accretion rate may correspond to
different DM parameters in their models, which could cause difficulty in
distinguishing them. In pining down a specific DM model, some other
constraints on DM [64] can be used in breaking this DM model degeneracy, then
a DM parameter in this model could be constrained from observed accretion
rate. In addition, some other effects such as baryonic matter accretion and
dynamical friction in a DM density spike [65], should be taken into
consideration in a real data analysis, which is neglected in above
calculations.

Apart from mass accretion, some other phenomena could also change the BH mass
and induce an orbital phase shift in PSR-BH binaries, such as superradiance
effect around a Kerr BH [66, 67, 68], which could extract $\mathcal{O}(10\%)$
mass of host BH to form a light boson cloud [69, 70]. Therefore, detailed
studies on these significant effects in PSR-BH binaries could shed light on
unknown physics and show a great potential of the PSR-BH binary.

## Acknowledgements

We would like to thank Lam Hui, Yi Wang, Henry Tye, and Leonardo Modesto for
the very helpful comments and advice.

##  Appendix A Weakly Interacting Massive Particles Accretion

Here we will introduce the Bondi accretion formula which applies to the WIMPs.
We assume spherical symmetry, and a steady state accretion, with the uniform
density $\rho_{\infty}$ and pressure $p_{\infty}$ at spatial infinity. Thus,
at infinity, the sound speed $c_{s}$ reads, $c_{s,\infty}=(\frac{\gamma
p_{\infty}}{\rho_{\infty}})^{\frac{1}{2}}$. The steady flow translates as

$\displaystyle\dot{M}=-4\pi r^{2}\rho u~{},$ (21)

where $u$ is the radial velocity. We can then write the equations for momentum
conservation,

$\displaystyle
u\frac{\text{d}u}{\text{d}r}+\frac{c_{s}^{2}}{\rho}\frac{\text{d}\rho}{\text{d}r}+\frac{GM}{r^{2}}=0~{},$
(22)

and the mass conservation,

$\displaystyle\frac{1}{\rho}\frac{\text{d}\rho}{\text{d}r}=-\frac{2}{r}-\frac{1}{u}\frac{\text{d}u}{\text{d}r}~{}.$
(23)

Substituting the latter into the former yields,

$\displaystyle\frac{1}{2}\left(1-\frac{c_{s}^{2}}{u^{2}}\right)\frac{\text{d}u^{2}}{\text{d}r}=\frac{-GM}{r^{2}}\left(1-\frac{2c_{s}^{2}r}{GM}\right)~{},$
(24)

which is the Bondi equation. We can contemplate the term between parentheses
on the LHS. As $r\to\infty$ and $c_{s}\to c_{s,\infty}$ that term is negative.
Whereas for $r\to 0$ it is positive again. Assuming continuity, there must be
a point where it vanishes. At that point the LHS vanishes simultaneously for
the $u^{2}(r_{s})\equiv u_{s}^{2}=c_{s}^{2}$. For our case, the physically
relevant solution is this so called transonic solution with the condition

$\displaystyle u^{2}\to 0\,\,\,{\rm as}\,\,\,r\to\infty~{}.$ (25)

From equation (24), the sonic point is at

$\displaystyle r_{s}=\frac{GM}{2c_{s}^{2}}~{}.$ (26)

Moreover, from equation (22), with the Polytropic condition

$\displaystyle P\sim\rho^{\gamma}~{},$ (27)

one simply derives

$\displaystyle\frac{u^{2}}{2}+\frac{c_{s}^{2}(r)}{\gamma-1}-\frac{GM}{r}=\frac{c_{s,\infty}^{2}}{\gamma-1}~{}.$
(28)

Then, evaluating this at the sonic point, we find,

$\displaystyle c_{s}(r_{s})=c_{s,\infty}\sqrt{\frac{2}{5-3\gamma}}~{}.$ (29)

This finally brings us to the Bondi accretion law,

$\displaystyle\dot{M}_{B}=4\pi\lambda_{B}(GM)^{2}\frac{\rho_{\infty}}{c_{s,\infty}^{3}}~{},$
(30)

with

$\displaystyle\lambda_{B}=\frac{1}{4}\left(\frac{2}{5-3\gamma}\right)^{\frac{5-3\gamma}{2(\gamma-1)}}~{}.$
(31)

Another way of expressing Eq.(30), as a function of the temperature per
particle mass, is using

$\displaystyle\frac{k_{B}T}{m}=\frac{P}{\rho}=\frac{c_{s}^{2}}{\gamma}~{},$
(32)

where $k_{B}$ is the Boltzman constant and $m$ is the WIMP mass. This,
inserted in (30), yields

$\displaystyle\dot{M}_{B}=4\pi\lambda_{B}(GM)^{2}\frac{\rho_{\infty}m^{\frac{3}{2}}}{\gamma^{\frac{1}{2}}\,(k_{B}T)^{\frac{3}{2}}}~{},$
(33)

Which is equivalent to Eq.(1).

##  Appendix B Hot Particle Dark Matter Accretion

We start from the continuity equation, and the vanishing of the divergence of
the stress-energy tensor,

$\displaystyle\nabla_{\mu}\left(\rho u^{\mu}\right)=0~{},$ (34)
$\displaystyle\nabla_{\mu}T^{\mu\nu}=0~{}.$ (35)

Here, $U^{\mu}$ is the DM fluid 4-velocity, and
$T^{\mu\nu}=\rho(1+h)u^{\mu}u^{\nu}+pg^{\mu\nu}$. We assume a Schwarzschild
metric background,

$\displaystyle\text{d}s^{2}=-(1-\frac{2M}{r})\text{d}t^{2}+\frac{1}{1-\frac{2M}{r}}\text{d}r^{2}+r^{2}\text{d}\Omega^{2}~{}.$
(36)

The corresponding metric determinant square root is thus

$\displaystyle\sqrt{-g}=r^{2}\sin\theta~{},$ (37)

and the four-velocity, defined in our spherically symmetric case as

$\displaystyle
u^{\mu}=\frac{\text{d}x^{\mu}}{\text{d}\tau}=(u^{t},u^{r},0,0)~{},$ (38)

with $u^{\mu}u_{\mu}=-1$. The components of the four-velocity thus relate to
each other as

$\displaystyle g_{tt}(u^{t})^{2}+g_{rr}(u_{r})^{2}=-1$ (39)
$\displaystyle(u^{t})^{2}=\frac{1-\frac{2M}{r}+(u^{r})^{2}}{(1-\frac{2M}{r})^{2}}~{}.$
(40)

Now we go back to Eq.(34) which yields,

$\displaystyle\frac{\text{d}}{\text{d}r}(r^{2}\rho u^{r})=0~{},$ (41)

and Eq.(35) yields,

$\displaystyle\frac{\text{d}}{\text{d}r}\left(r^{2}\rho(1+\frac{\gamma}{\gamma-1}\Theta)u^{t}u^{r}\right)=0~{}.$
(42)

Once integrated they reduce to

$\displaystyle\dot{M}=4\pi r^{2}\rho u^{r}~{},$ (43)

and

$\displaystyle(1+\frac{\gamma\Theta}{\gamma-1})\left(1-\frac{2M}{r}+|u^{r}|^{2}\right)=(1+\frac{\gamma\Theta_{\infty}}{\gamma-1})~{}.$
(44)

Combining the two equations, and following the same procedure as in Appendix.
A yields Eq. (3).

##  Appendix C Ultralight Dark Matter Accretion

In [33], Unruh considered the Klein-Gordon equation on a Schwarzschild BH
geometry,

$\displaystyle g^{\mu\nu}\phi_{,\mu;\nu}+m^{2}\phi=0~{}.$ (45)

With the separation of variables

$\displaystyle\phi(t,r,\theta,\varphi)={\rm e}^{-i\omega t}f_{\omega
l}(r)Y_{lm}(\theta,\varphi)~{},$ (46)

one is left with,

$\displaystyle\\!\\!\\!\\!\left[\frac{\rho^{2}}{r^{2}}\frac{\text{d}}{\text{d}r}(r^{2}\rho^{2}\frac{\text{d}}{\text{d}r})\\!+\\!\omega^{2}\\!-\\!\rho^{2}\\!\left(\\!m^{2}\\!\\!+\\!\frac{l(l+1)}{r^{2}}\\!\right)\\!\right]\\!\\!f_{\omega
l}(r)=0~{},$ (47)

where $\rho=(1-\frac{2M}{r})^{\frac{1}{2}}$. The equation is not solvable
analytically as it is, but, it can be approximated in three phases, far from
the BH, at intermediate distance, and very close. Each of these cases can be
analytically solved, then matched by imposing reasonable boundary conditions
at infinity and matching the boundary conditions until reaching the BH
surrounding. Having the solution near the BH, one can compute the current term

$\displaystyle
J^{\mu}=\frac{1}{2}i\phi^{*}\overleftrightarrow{\partial}^{\mu}\phi~{}.$ (48)

The number of particles absorbed by the BH will be a surface integral of the
radial current

$\displaystyle N=-\int_{\mathcal{S}}\sqrt{-g}J^{r}\text{d}\mathcal{S}~{}.$
(49)

Plugging the solution in we get the cross section

$\displaystyle\sigma(\omega,v)=\frac{N}{\omega v}=\frac{(4\pi
M)^{2}(1+v^{2})2M\omega}{v^{2}(1-{\rm e}^{-\frac{2\pi
M\omega(1+v^{2})}{v}})}~{}.$ (50)

##  Appendix D Gravitational Potential Energy Changing Rate

The gravitational potential energy in a PSR-BH binary system can be estimated
as

$\displaystyle E_{p}=-\frac{Gm_{\mathrm{p}}M_{\mathrm{B}}}{r}~{}.$ (51)

Due to the DM accretion into BH, the gravitational potential energy is
transferred from the DM to the BH, which results in the increasing in BH mass
and decreasing in gravitational potential energy, as follows,

$\displaystyle\frac{\text{d}E_{p}}{\text{d}t}=-\frac{Gm_{\mathrm{p}}}{r}\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}~{}.$
(52)

During the orbital motion of the binary, the distance $r$ between two
components can be expressed in polar coordinate as follows,

$\displaystyle r=\frac{a(1-e^{2})}{1+e\cos\psi}~{}.$ (53)

Then Eq. (52) can be written as

$\displaystyle\frac{\text{d}E_{p}}{\text{d}t}=-Gm_{\mathrm{p}}\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}\frac{1+e\cos\psi}{a(1-e^{2})}~{}.$
(54)

The average gravitational potential energy changing rate over one period is

$\displaystyle\left\langle\frac{\text{d}E_{p}}{\text{d}t}\right\rangle$
$\displaystyle=-\frac{1}{T}\int_{0}^{T}Gm_{\mathrm{p}}\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}\frac{1+e\cos\psi}{a(1-e^{2})}\text{d}t$
$\displaystyle=-\frac{1}{T}\int_{0}^{2\pi}Gm_{\mathrm{p}}\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}\frac{1+e\cos\psi}{a(1-e^{2})}\frac{\text{d}\psi}{\dot{\psi}}~{}.$
(55)

Here $T$ is orbital period and $\dot{\psi}$ is the angular velocity, which
follow

$\displaystyle T$ $\displaystyle=\frac{2\pi
a^{3/2}}{G^{1/2}(m_{\mathrm{p}}+M_{\mathrm{B}})^{1/2}}~{},$
$\displaystyle\dot{\psi}$
$\displaystyle=\frac{[G(m_{\mathrm{p}}+M_{\mathrm{B}})a(1-e^{2})]^{1/2}}{r^{2}}~{}.$
(56)

Then, the integral Eq. (D) gives the result

$\displaystyle\left\langle\frac{\text{d}E_{p}}{\text{d}t}\right\rangle=-\frac{Gm_{\mathrm{p}}}{a}\frac{\text{d}M_{\mathrm{B}}}{\text{d}t}~{}.$
(57)

## References

  * [1] Planck Collaboration, N. Aghanim et al., “Planck 2018 results. VI. Cosmological parameters,” Astron. Astrophys. 641 (2020) A6, arXiv:1807.06209 [astro-ph.CO]. [Erratum: Astron.Astrophys. 652, C4 (2021)].
  * [2] V. C. Rubin and W. K. Ford, Jr., “Rotation of the Andromeda Nebula from a Spectroscopic Survey of Emission Regions,” Astrophys. J. 159 (1970) 379–403.
  * [3] S. M. Faber and J. S. Gallagher, “Masses and mass-to-light ratios of galaxies,” Ann. Rev. Astron. Astrophys. 17 (1979) 135–183.
  * [4] M. S. Turner, G. Steigman, and L. M. Krauss, “The Flatness of the Universe: Reconciling Theoretical Prejudices with Observational Data,” Phys. Rev. Lett. 52 (1984) 2090–2093.
  * [5] R. A. Flores and J. R. Primack, “Observational and theoretical constraints on singular dark matter halos,” Astrophys. J. Lett. 427 (1994) L1–4, arXiv:astro-ph/9402004.
  * [6] B. Moore, “Evidence against dissipationless dark matter from observations of galaxy haloes,” Nature 370 (1994) 629.
  * [7] F. C. van den Bosch, B. E. Robertson, J. J. Dalcanton, and W. J. G. de Blok, “Constraints on the structure of dark matter halos from the rotation curves of low surface brightness galaxies,” Astron. J. 119 (2000) 1579–1591, arXiv:astro-ph/9911372.
  * [8] A. A. Klypin, A. V. Kravtsov, O. Valenzuela, and F. Prada, “Where are the missing Galactic satellites?,” Astrophys. J. 522 (1999) 82–92, arXiv:astro-ph/9901240.
  * [9] B. Moore, S. Ghigna, F. Governato, G. Lake, T. R. Quinn, J. Stadel, and P. Tozzi, “Dark matter substructure within galactic halos,” Astrophys. J. Lett. 524 (1999) L19–L22, arXiv:astro-ph/9907411.
  * [10] H.-Y. Chiu, “Symmetry between particle and anti-particle populations in the universe,” Phys. Rev. Lett. 17 (1966) 712.
  * [11] B. W. Lee and S. Weinberg, “Cosmological Lower Bound on Heavy Neutrino Masses,” Phys. Rev. Lett. 39 (1977) 165–168.
  * [12] E. G. M. Ferreira, “Ultra-light dark matter,” Astron. Astrophys. Rev. 29 no. 1, (2021) 7, arXiv:2005.03254 [astro-ph.CO].
  * [13] S. Hawking, “Gravitationally collapsed objects of very low mass,” Mon. Not. Roy. Astron. Soc. 152 (1971) 75.
  * [14] T. Clifton, P. G. Ferreira, A. Padilla, and C. Skordis, “Modified Gravity and Cosmology,” Phys. Rept. 513 (2012) 1–189, arXiv:1106.2476 [astro-ph.CO].
  * [15] C. Boehm and R. Schaeffer, “Constraints on dark matter interactions from structure formation: Damping lengths,” Astron. Astrophys. 438 (2005) 419–442, arXiv:astro-ph/0410591.
  * [16] DES Collaboration, E. O. Nadler et al., “Milky Way Satellite Census. III. Constraints on Dark Matter Properties from Observations of Milky Way Satellite Galaxies,” Phys. Rev. Lett. 126 (2021) 091101, arXiv:2008.00022 [astro-ph.CO].
  * [17] B. Carr, K. Kohri, Y. Sendouda, and J. Yokoyama, “Constraints on primordial black holes,” Rept. Prog. Phys. 84 no. 11, (2021) 116902, arXiv:2002.12778 [astro-ph.CO].
  * [18] S. Peirani and J. A. de Freitas Pacheco, “Dark Matter Accretion into Supermassive Black Holes,” Phys. Rev. D 77 (2008) 064023, arXiv:0802.2041 [astro-ph].
  * [19] R. Brito, V. Cardoso, and H. Okawa, “Accretion of dark matter by stars,” Phys. Rev. Lett. 115 no. 11, (2015) 111301, arXiv:1508.04773 [gr-qc].
  * [20] N. J. Outmezguine, O. Slone, W. Tangarife, L. Ubaldi, and T. Volansky, “Accretion of Dissipative Dark Matter onto Active Galactic Nuclei,” JHEP 11 (2018) 005, arXiv:1807.04750 [astro-ph.GA].
  * [21] K. Malhan, M. Valluri, and K. Freese, “Probing the nature of dark matter with accreted globular cluster streams,” Monthly Notices of the Royal Astronomical Society 501 no. 1, (2021) 179–200, arXiv:2005.12919 [astro-ph.CO].
  * [22] R. D’Agostino, R. Giambò, and O. Luongo, “Constraining primordial black holes as a fraction of dark matter through accretion disk luminosity,” Phys. Rev. D 107 no. 4, (2023) 043032, arXiv:2204.02098 [gr-qc].
  * [23] D. C. Backer, S. R. Kulkarni, C. Heiles, M. Davis, and W. Goss, “A millisecond pulsar,” Nature 300 no. 5893, (1982) 615–618.
  * [24] M. Kramer et al., “Strong-Field Gravity Tests with the Double Pulsar,” Phys. Rev. X 11 no. 4, (2021) 041050, arXiv:2112.06795 [astro-ph.HE].
  * [25] Y. Shao and X.-D. Li, “Black Hole/Pulsar Binaries in the Galaxy,” Mon. Not. Roy. Astron. Soc. 477 no. 1, (2018) L128–L132, arXiv:1804.06014 [astro-ph.HE].
  * [26] D. Chattopadhyay, S. Stevenson, J. R. Hurley, M. Bailes, and F. Broekgaarden, “Modelling neutron star–black hole binaries: future pulsar surveys and gravitational wave detectors,” Mon. Not. Roy. Astron. Soc. 504 no. 3, (2021) 3682–3710, arXiv:2011.13503 [astro-ph.HE].
  * [27] R. Nan, D. Li, C. Jin, Q. Wang, L. Zhu, W. Zhu, H. Zhang, Y. Yue, and L. Qian, “The five-hundred-meter aperture spherical radio telescope (fast) project,” International Journal of Modern Physics D 20 no. 06, (2011) 989–1024.
  * [28] R. Booth, W. De Blok, J. Jonas, and B. Fanaroff, “Meerkat key project science, specifications, and proposals,” arXiv:0910.2935 [astro-ph.IM].
  * [29] P. E. Dewdney, P. J. Hall, R. T. Schilizzi, and T. J. L. Lazio, “The square kilometre array,” Proceedings of the IEEE 97 no. 8, (2009) 1482–1496.
  * [30] P. Amaro-Seoane, H. Audley, S. Babak, J. Baker, E. Barausse, P. Bender, E. Berti, P. Binetruy, M. Born, D. Bortoluzzi, et al., “Laser interferometer space antenna,” arXiv:1702.00786 [astro-ph.IM].
  * [31] H. Bondi, “On spherically symmetrical accretion,” Mon. Not. Roy. Astron. Soc. 112 (1952) 195.
  * [32] F. Michel, “Accretion of matter by condensed objects,” Astrophys Space Sci 15 (1972) 153–160.
  * [33] W. G. Unruh, “Absorption Cross-Section of Small Black Holes,” Phys. Rev. D 14 (1976) 3251–3259.
  * [34] L. Roszkowski, E. M. Sessolo, and S. Trojanowski, “WIMP dark matter candidates and searches—current status and future prospects,” Rept. Prog. Phys. 81 no. 6, (2018) 066201, arXiv:1707.06277 [hep-ph].
  * [35] A. Aguayo-Ortiz, E. Tejeda, O. Sarbach, and D. López-Cámara, “Spherical accretion: Bondi, Michel and rotating black holes,” arXiv:2102.12529 [astro-ph.HE].
  * [36] C. Armendariz-Picon and J. T. Neelakanta, “How Cold is Cold Dark Matter?,” JCAP 03 (2014) 049, arXiv:1309.6971 [astro-ph.CO].
  * [37] P. P. Avelino and V. M. C. Ferreira, “Constraints on the dark matter sound speed from galactic scales: the cases of the Modified and Extended Chaplygin Gas,” Phys. Rev. D 91 no. 8, (2015) 083508, arXiv:1502.07583 [astro-ph.CO].
  * [38] M. Kawasaki, N. Sugiyama, and T. Yanagida, “Primordial black hole formation in a double inflation model in supergravity,” Phys. Rev. D 57 (1998) 6050–6056, arXiv:hep-ph/9710259.
  * [39] K. Jedamzik, “Could MACHOS be primordial black holes formed during the QCD epoch?,” Phys. Rept. 307 (1998) 155–162, arXiv:astro-ph/9805147.
  * [40] A. Zakharov, “Effective particle capture cross section of a schwarzschild black-hole,” Soviet Astronomy 32 (1988) 456.
  * [41] J. Binney and S. Tremaine, Galactic Dynamics: Second Edition.  Princeton Series in Astrophysics. Princeton University Press, 2011.  https://books.google.at/books?id=6mF4CKxlbLsC.
  * [42] K. Liu, R. P. Eatough, N. Wex, and M. Kramer, “Pulsar–black hole binaries: prospects for new gravity tests with future radio telescopes,” Mon. Not. Roy. Astron. Soc. 445 no. 3, (2014) 3115–3132, arXiv:1409.3882 [astro-ph.GA].
  * [43] J. H. Taylor, “Binary pulsars and relativistic gravity,” Rev. Mod. Phys. 66 (1994) 711–719.
  * [44] G. Hobbs, R. Edwards, and R. Manchester, “Tempo2, a new pulsar timing package. 1. overview,” Mon. Not. Roy. Astron. Soc. 369 (2006) 655–672, arXiv:astro-ph/0603381.
  * [45] Q. Ding, X. Tong, and Y. Wang, “Gravitational Collider Physics via Pulsar-Black Hole Binaries,” Astrophys. J. 908 no. 1, (2021) 78, arXiv:2009.11106 [astro-ph.HE].
  * [46] X. Tong, Y. Wang, and H.-Y. Zhu, “Gravitational Collider Physics via Pulsar–Black Hole Binaries II: Fine and Hyperfine Structures Are Favored,” Astrophys. J. 924 no. 2, (2022) 99, arXiv:2106.13484 [astro-ph.HE].
  * [47] P. C. Peters, “Gravitational Radiation and the Motion of Two Point Masses,” Phys. Rev. 136 (1964) B1224–B1232.
  * [48] J. F. Navarro, C. S. Frenk, and S. D. M. White, “The Structure of cold dark matter halos,” Astrophys. J. 462 (1996) 563–575, arXiv:astro-ph/9508025.
  * [49] J. F. Navarro, C. S. Frenk, and S. D. M. White, “A Universal density profile from hierarchical clustering,” Astrophys. J. 490 (1997) 493–508, arXiv:astro-ph/9611107.
  * [50] H.-N. Lin and X. Li, “The Dark Matter Profiles in the Milky Way,” Mon. Not. Roy. Astron. Soc. 487 no. 4, (2019) 5679–5684, arXiv:1906.08419 [astro-ph.GA].
  * [51] D. Clausen, S. Sigurdsson, and D. F. Chernoff, “Dynamically formed black hole+millisecond pulsar binaries in globular clusters,” Mon. Not. Roy. Astron. Soc. 442 no. 1, (2014) 207–219, arXiv:1404.7502 [astro-ph.HE].
  * [52] D. Merritt, T. Alexander, S. Mikkola, and C. M. Will, “Stellar Dynamics of Extreme-Mass-Ratio Inspirals,” Phys. Rev. D 84 (2011) 044024, arXiv:1102.3180 [astro-ph.CO].
  * [53] G. Fragione, I. Ginsburg, and B. Kocsis, “Gravitational Waves and Intermediate-mass Black Hole Retention in Globular Clusters,” Astrophys. J. 856 no. 2, (2018) 92, arXiv:1711.00483 [astro-ph.GA].
  * [54] B. Devecchi, M. Colpi, M. Mapelli, and A. Possenti, “Millisecond pulsars around intermediate-mass black holes in globular clusters,” Mon. Not. Roy. Astron. Soc. 380 (2007) 691, arXiv:0706.1656 [astro-ph].
  * [55] F. Zhang, Y. Lu, and Q. Yu, “On the Existence of Pulsars in the Vicinity of the Massive Black Hole in the Galactic Center,” Astrophys. J. 784 (2014) 106, arXiv:1402.2505 [astro-ph.GA].
  * [56] S. Toonen, A. Hamers, and S. Portegies Zwart, “The evolution of hierarchical triple star-systems,” Computational Astrophysics and Cosmology 3 no. 1, (2016) 1–36.
  * [57] D. F. Chernoff and L. S. Finn, “Gravitational radiation, inspiraling binaries, and cosmology,” Astrophys. J. Lett. 411 (1993) L5–L8, arXiv:gr-qc/9304020.
  * [58] P. Ajith et al., “Inspiral-merger-ringdown waveforms for black-hole binaries with non-precessing spins,” Phys. Rev. Lett. 106 (2011) 241101, arXiv:0909.2867 [gr-qc].
  * [59] X.-J. Zhu, E. Howell, T. Regimbau, D. Blair, and Z.-H. Zhu, “Stochastic Gravitational Wave Background from Coalescing Binary Black Holes,” Astrophys. J. 739 (2011) 86, arXiv:1104.3565 [gr-qc].
  * [60] P. Ajith et al., “A Template bank for gravitational waveforms from coalescing binary black holes. I. Non-spinning binaries,” Phys. Rev. D 77 (2008) 104017, arXiv:0710.2335 [gr-qc]. [Erratum: Phys.Rev.D 79, 129901 (2009)].
  * [61] L. Hui, J. P. Ostriker, S. Tremaine, and E. Witten, “Ultralight scalars as cosmological dark matter,” Phys. Rev. D 95 no. 4, (2017) 043541, arXiv:1610.08297 [astro-ph.CO].
  * [62] I. De Martino, T. Broadhurst, S. H. H. Tye, T. Chiueh, and H.-Y. Schive, “Dynamical Evidence of a Solitonic Core of $10^{9}M_{\odot}$ in the Milky Way,” Phys. Dark Univ. 28 (2020) 100503, arXiv:1807.08153 [astro-ph.GA].
  * [63] P. Bhattacharjee, S. Chaudhury, and S. Kundu, “Rotation Curve of the Milky Way out to $\sim$ 200 kpc,” Astrophys. J. 785 (2014) 63, arXiv:1310.2659 [astro-ph.GA].
  * [64] J. Billard et al., “Direct detection of dark matter—APPEC committee report*,” Rept. Prog. Phys. 85 no. 5, (2022) 056201, arXiv:2104.07634 [hep-ex].
  * [65] M. H. Chan and C. M. Lee, “Indirect Evidence for Dark Matter Density Spikes around Stellar-mass Black Holes,” Astrophys. J. Lett. 943 no. 2, (2023) L11, arXiv:2212.05664 [astro-ph.HE].
  * [66] Y. B. Zel’Dovich, “Generation of Waves by a Rotating Body,” Soviet Journal of Experimental and Theoretical Physics Letters 14 (Aug., 1971) 180.
  * [67] Y. B. Zel’Dovich, “Amplification of Cylindrical Electromagnetic Waves Reflected from a Rotating Body,” Soviet Journal of Experimental and Theoretical Physics 35 (Jan., 1972) 1085.
  * [68] R. Brito, V. Cardoso, and P. Pani, “Superradiance: New Frontiers in Black Hole Physics,” Lect. Notes Phys. 906 (2015) pp.1–237, arXiv:1501.06570 [gr-qc].
  * [69] K. K. Y. Ng, S. Vitale, O. A. Hannuksela, and T. G. F. Li, “Constraints on Ultralight Scalar Bosons within Black Hole Spin Measurements from the LIGO-Virgo GWTC-2,” Phys. Rev. Lett. 126 no. 15, (2021) 151102, arXiv:2011.06010 [gr-qc].
  * [70] L. Hui, Y. T. A. Law, L. Santoni, G. Sun, G. M. Tomaselli, and E. Trincherini, “Black hole superradiance with (dark) matter accretion,” arXiv:2208.06408 [gr-qc].


# Inverse Design Method Starting for Aperture Field from Phaseless Radiation
Pattern Based on Spectral Analysis

Chuan-Sheng Chen, Ren Wang, , Jin-Pin Liu, and Bing-Zhong Wang,  _This work
has been submitted to the IEEE for possible publication. Copyright may be
transferred without notice, after which this version may no longer be
accessible._ Manuscript received on Nov 27, 2022. This work was supported in
part by the National Natural Science Foundation of China under Grants 62171081
and 61901086. (Corresponding authors: Bing-Zhong Wang.) All four authors are
with the Institute of Applied Physics, University of Electronic Science and
Technology of China (UESTC), Chengdu, 611731, China, (e-mail:
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>and
bzwang@uestc.edu.cn).

###### Abstract

Existing electromagnetic inverse design methods are often established in the
spacial domain. This communication presents an inverse design method, which
can design aperture field for the desired phaseless radiation pattern, from
the spectral domain perspective. In addition, it naturally adapts to the
polarization constraint. Specifically, the inverse design can be converted
into solving the first kind of Fredholm integral equation, using the spectral
domain method. To deal with the ill-posedness of such integral equation, we
apply modal expansion to the integrand. To cope with the non-linearity
introduced by phaseless, we use a multi-objective optimization algorithm to
obtain the coefficients in the modal expansion. Finally, we use this method to
break the cosine rule that the directivity of 2D arrays drops as cosine of the
angle, which is a puzzle in wide-angle scanning. The numerical simulation
results meet expectations and illustrate the feasibility of the method.

###### Index Terms:

inverse problem; spectral domain method; aperture field; radiation pattern

##  I Introduction

The radiation pattern of antennas are always guided by the requirements of
electromagnetic systems. However, existing design methods are either only
applicable to specific antennas, which do not have good universality [1, 2],
or use optimization to achieve simple requirements which lack flexibility[3,
4]. A general design method for arbitrary radiation pattern, especially
phaseless pattern is worth studying. The main difficulty is how to effectively
handle the non-linearity introduced by phaseless. In addition, the
polarization criteria of the radiation pattern should not be ignored.
Obviously, the design problem starting from phaseless arbitrary pattern is a
typical inverse design problem that aims to find an unknown cause from its
known effect [5].

In fact, inverse design has gradually been favored by scholars [6, 7, 8, 9,
10]. The inverse source framework presented in [9] is based on the well-known
electric field integral equation (EFIE), and belongs to the spacial domain
analysis, as shown in (1). A conjugate gradient (CG) algorithm [11] is used to
iteratively update $\boldsymbol{J}$ and $\boldsymbol{M}$ to solve the integral
equation. Based on this framework, [10] used the three-layer admittance model
to design a metasurface for one desired radiation pattern.

$\displaystyle\mathcal{G}(\boldsymbol{J},\boldsymbol{M})=$
$\displaystyle-j\eta_{0}k_{0}\int\limits_{\Sigma^{+}}{\left[\boldsymbol{J}(\boldsymbol{r}^{\prime})+\frac{1}{k_{0}^{2}}\nabla\nabla_{s}^{\prime}\cdot\boldsymbol{J}(\boldsymbol{r}^{\prime})\right]}g(\boldsymbol{r},\boldsymbol{r}^{\prime})ds^{\prime}$
(1)
$\displaystyle-\nabla\times\int\limits_{\Sigma^{+}}{\boldsymbol{M}}(\boldsymbol{r}^{\prime})g(\boldsymbol{r},\boldsymbol{r}^{\prime})ds^{\prime}$

The spacial domain framework proposed in [9] has good generality. However,
when the design target is phaseless pattern, some calculations in (1) may be
redundant. This part of the computational burden is not negligible during the
iterative solution. In addition, the applicability of CG algorithm to non-
linear problems remains to be discussed. When dealing with dual polarization
constraints, it is irregular but important to balance the weights of different
polarizations in the objective function.

This communication proposes an inverse design method which can deal with the
phaseless pattern well. It reduces the computational burden by combining the
far-field asymptotics of the spectral domain method, and is naturally adapt
different polarization constraint by combining with the coordinate
transformation. We use multi-objective optimization algorithms to deal with
the non-linearity introduced by phaseless. To illustrate the significance of
our method, we use it to break the cosine rule that the directivity of 2D
arrays drops as cosine of the scanning angle, and the numerical results show
that we succeeded.

The rest of this communication is presented as follows. Section II introduces
the basic forward problem. Section III introduces the inverse design
methodology, which ultimately translates the inverse design problem into
solving the first kind of two dimensional Fredholm integral equation. Section
IV analyses how to solve the integral equations, and Section V demonstrates
the capability of the method through an example breaking the cosine rule.

##  II Forward Problem Statement

The general aperture radiation problem is calculating the radiation field from
the known aperture field $\boldsymbol{E}_{\mathrm{ape}}(x,y)|_{z=0}$. Its
schematic diagram is shown in Fig. 1. The analysis below consider only the
upper half-space ($z>0$), that is $\theta\in[0,{\pi}/{2}]$.
$\boldsymbol{E}_{\mathrm{ape}}$, which has a definition domain
$\left\\{(x,y)|x\in[-{a}/{2},{a}/{2}],y\in[-{b}/{2},{b}/{2}]\right\\}$, locate
in the infinite PEC plane ($z=0$ plane). According to the uniqueness theorem,
as long as the tangential components of $\boldsymbol{E}_{\mathrm{ape}}(x,y)$
($E_{\mathrm{ape},x}$ and $E_{\mathrm{ape},y}$) are known,
$\boldsymbol{E}(x,y,z)|_{z>0}$ is uniquely determined. This forward problem
can be analyzed from the spatial or spectral domains.

Figure 1: Schematic diagram of the general forward aperture problem

In the spatial domain analysis, it can be transformed into a problem of
magnetic current
$(\boldsymbol{M}_{s}=-2\boldsymbol{n}\times\boldsymbol{E}_{a})$ radiation in
free space[12]. $\boldsymbol{E}(x,y,z)|_{z>0}$ can be calculated by using the
Green’s function method shown in (1) or by introducing the auxiliary potential
functions $\boldsymbol{A}$ and $\boldsymbol{F}$.

In the spectral domain, as written in (2), the monochromatic wave field
$\boldsymbol{E}(x,y,z)$ radiated by the aperture can be regarded as a
superposition of plane waves of the form
$\boldsymbol{f}\left(\theta,\phi\right)\exp({-j\boldsymbol{k}\cdot\boldsymbol{r}})$[13].

$\displaystyle\boldsymbol{E}(x,y,z)=\frac{1}{4\pi^{2}}\int_{-\infty}^{\infty}{\int_{-\infty}^{\infty}{\boldsymbol{f}\left(\theta,\phi\right)e^{-j\boldsymbol{k}\cdot\boldsymbol{r}}}}\mathrm{d}k_{x}\mathrm{d}k_{y}$
(2)

where $k_{x}=k\sin\theta\cos\phi$, $k_{y}=k\sin\theta\sin\phi$, and
$k_{z}^{2}=k^{2}-k_{x}^{2}-k_{y}^{2}$. $f_{x}$ and $f_{y}$ are given by (3)
and are the components of $\boldsymbol{f}$.

$\displaystyle f_{x}(\theta,\phi)$
$\displaystyle=\int_{-\frac{b}{2}}^{\frac{b}{2}}\int_{-\frac{a}{2}}^{\frac{a}{2}}E_{\mathrm{ape},x}(x^{\prime},y^{\prime})\cdot
e^{j(k_{x}x^{\prime}+k_{y}y^{\prime})}\mathrm{d}x^{\prime}\mathrm{d}y^{\prime}$
(3a) $\displaystyle f_{y}(\theta,\phi)$
$\displaystyle=\int_{-\frac{b}{2}}^{\frac{b}{2}}\int_{-\frac{a}{2}}^{\frac{a}{2}}E_{\mathrm{ape},y}(x^{\prime},y^{\prime})\cdot
e^{j(k_{x}x^{\prime}+k_{y}y^{\prime})}\mathrm{d}x^{\prime}\mathrm{d}y^{\prime}$
(3b)

The far-field asymptotic evaluation of (2) is (4), by using the Stationary
Phase method.

$\boldsymbol{E}(r,\theta,\phi)\approx\ j\frac{ke^{-jkr}}{2\pi
r}\boldsymbol{e}_{\mathrm{far}}(\theta,\phi)$ (4a)
$\displaystyle\boldsymbol{e}_{\mathrm{far}}(\theta,\phi)$
$\displaystyle=\boldsymbol{\hat{a}}_{\theta}e_{\mathrm{far},\theta}(\theta,\phi)+\boldsymbol{\hat{a}}_{\phi}e_{\mathrm{far},\phi}(\theta,\phi)$
(4b)
$\displaystyle=\boldsymbol{\hat{a}}_{\theta}\left(f_{x}\cos\phi+f_{y}\sin\phi\right)$
$\displaystyle\quad+\boldsymbol{\hat{a}}_{\phi}[\cos\theta\left(-f_{x}\sin\phi+f_{y}\cos\phi\right)]$

In most cases, we pay more attention to the pattern
$\boldsymbol{e}_{\mathrm{far}}(\theta,\phi)$ than $j\frac{ke^{-jkr}}{2\pi r}$,
the coefficient on distance. Obviously, the $\boldsymbol{e}_{\mathrm{far}}$
here includes both polarization and phase information.

##  III Inverse Design Methodology

This section focuses on how to introduce phaseless pattern into the existing
formulations, so as to convert phaseless pattern targets into constraints on
aperture field. It is important to emphasize that, during the inverse design
process, the given is the phaseless pattern
$\left|e_{\mathrm{far},\theta}\right|$ and
$\left|e_{\mathrm{far},\phi}\right|$, which can not be directly brought into
(5).

$\displaystyle f_{x}(\theta,\phi)$
$\displaystyle=e_{\mathrm{far},\theta}(\theta,\phi)\cos\theta-
e_{\mathrm{far},\phi}(\theta,\phi){\sin\phi}/{\cos\theta}$ (5a) $\displaystyle
f_{y}(\theta,\phi)$
$\displaystyle=e_{\mathrm{far},\theta}(\theta,\phi)\sin\phi+e_{\mathrm{far},\phi}(\theta,\phi){\cos\phi}/{\cos\theta}$
(5b)

So we obtain (6) by taking modulus of both sides of (5). $\gamma$ in (6) is a
function of $\theta$ and $\phi$, which represents the phase difference of two
polarization components at a certain direction $(\theta,\phi)$.

$\displaystyle\begin{split}|f_{x}(\theta,\phi)|=&\sqrt{|e_{\mathrm{far},\theta}|^{2}|\cos\theta|^{2}+|e_{\mathrm{far},\phi}|^{2}|\sin\phi/\cos\theta|^{2}}\\\
&\quad\overline{-2|e_{\mathrm{far},\theta}||e_{\mathrm{far},\phi}||\sin\phi|\cos\gamma}\end{split}$
(6a)
$\displaystyle\begin{split}|f_{y}(\theta,\phi)|=&\sqrt{|e_{\mathrm{far},\theta}|^{2}|\sin\phi|^{2}+|e_{\mathrm{far},\phi}|^{2}|\cos\phi/\cos\theta|^{2}}\\\
&\quad\overline{+2|e_{\mathrm{far},\theta}||e_{\mathrm{far},\phi}||\sin\phi\cos\phi/\cos\theta|\cos\gamma}\end{split}$
(6b)

Combined with engineering applications, we also give another expression of (6)
for horizontal and vertical polarization ($|e_{\mathrm{far},h}|$ and
$|e_{\mathrm{far},v}|$), by using Ludwig3 coordinate system, as shown in (7).

$\displaystyle|f_{x}(\theta,\phi)|=\sqrt{|l_{1}|^{2}+|l_{2}|^{2}+2|l_{1}||l_{2}|\cos\gamma}$
(7a)
$\displaystyle|f_{y}(\theta,\phi)|=\sqrt{|l_{3}|^{2}+|l_{4}|^{2}+2|l_{3}||l_{4}|\cos\gamma}$
(7b)

where

$\displaystyle l_{1}$
$\displaystyle=e_{\mathrm{far},h}\left(\cos^{2}\phi+\sin^{2}\phi/\cos\theta\right);$
$\displaystyle l_{2}$
$\displaystyle=e_{\mathrm{far},v}\sin\phi\cos\phi\left(1-1/\cos\theta\right);$
$\displaystyle l_{3}$
$\displaystyle=e_{\mathrm{far},h}\sin\phi\cos\phi\left(1-1/\cos\theta\right);$
$\displaystyle l_{4}$
$\displaystyle=e_{\mathrm{far},v}\left(\sin^{2}\phi+\cos^{2}\phi/\cos\theta\right);$

Several common cases regarding (7) are discussed next.

  * •

Case I: Linear Polarization

When only a single polarization is needed, another orthogonal polarization can
be set to 0. For example, when only the vertical polarization
$e_{\mathrm{far},v}$ is considered, (7) degenerates to (8).

$|f_{x}(\theta,\phi)|=|l_{2}|;\ |f_{y}(\theta,\phi)|=|l_{4}|$ (8)

  * •

Case II: Circular Polarization

Circular polarization requires that the modulus values of the two polarization
components are equal, and the phase difference is $90^{\circ}$. That is
$|e_{\mathrm{far},h}|=|e_{\mathrm{far},v}|$, and $\gamma=90^{\circ}$ in any
direction $(\theta,\phi)$. Then, (7) degenerates to (9).

$|f_{x}(\theta,\phi)|=\sqrt{|l_{1}|^{2}+|l_{2}|^{2}};\
|f_{y}(\theta,\phi)|=\sqrt{|l_{3}|^{2}+|l_{4}|^{2}}$ (9)

  * •

Case III: Elliptical Polarization

Once $\gamma$ is determined, equation (7) itself corresponds to elliptical
polarization.

The $|f_{x}|$ and $|f_{y}|$ in (6) and (7) act as bridges to deliver the
phaseless pattern constraints to the aperture field constraints. Hence, the
problem of solving the aperture field can be summarized in two steps:

  1. 1.

Get $|f_{x}|$ and $|f_{y}|$ from the desired $|e_{\mathrm{far},\theta}|$ and
$|e_{\mathrm{far},\phi}|$ or $|e_{\mathrm{far},h}|$ and $|e_{\mathrm{far},v}|$
by using (6) or (7).

  2. 2.

Find a suitable set of solutions($E_{\mathrm{ape},x}$ and
$E_{\mathrm{ape},y}$) from $|f_{x}|$ and $|f_{y}|$ according to (10).

$\displaystyle|f_{x}(\theta,\phi)|$
$\displaystyle=\left|\int_{-\frac{b}{2}}^{\frac{b}{2}}\int_{-\frac{a}{2}}^{\frac{a}{2}}E_{\mathrm{ape},x}(x^{\prime},y^{\prime})\cdot
e^{j(k_{x}x^{\prime}+k_{y}y^{\prime})}\mathrm{d}x^{\prime}\mathrm{d}y^{\prime}\right|$
(10a) $\displaystyle|f_{y}(\theta,\phi)|$
$\displaystyle=\left|\int_{-\frac{b}{2}}^{\frac{b}{2}}\int_{-\frac{a}{2}}^{\frac{a}{2}}E_{\mathrm{ape},y}(x^{\prime},y^{\prime})\cdot
e^{j(k_{x}x^{\prime}+k_{y}y^{\prime})}\mathrm{d}x^{\prime}\mathrm{d}y^{\prime}\right|$
(10b)

During Step 2, the inverse Fourier transform can not be used in (10). The root
cause is the non-linearity introduced by the phaseless. Therefore, we can only
treat (10) as the first kind of two-dimensional Fredholm integral equations
after the modulo operation, which is difficult to obtain analytical
solutions[14].

$\displaystyle\left|f_{x}(\theta,\phi)\right|=\left|\int_{-\frac{b}{2}}^{\frac{b}{2}}{\int_{-\frac{a}{2}}^{\frac{a}{2}}{}}\right.$
$\displaystyle\left\\{\sum_{n=1}^{N}\left[\left(\sum_{m=1}^{M}\alpha_{mn}^{x}\sin\left(\frac{m\pi}{a}(x^{\prime}+\frac{a}{2})\right)+\sum_{m=0}^{M}\beta_{mn}^{x}\cos\left(\frac{m\pi}{a}(x^{\prime}+\frac{a}{2})\right)\right)\sin\left(\frac{n\pi}{b}(y^{\prime}+\frac{b}{2})\right)\right]\right\\}$
(11a)
$\displaystyle\left.e^{j(k_{x}x^{\prime}+k_{y}y^{\prime})}\mathrm{d}x^{\prime}\mathrm{d}y^{\prime}\right|$
$\displaystyle|f_{y}(\theta,\phi)|=\left|\int_{-\frac{b}{2}}^{\frac{b}{2}}{\int_{-\frac{a}{2}}^{\frac{a}{2}}{}}\right.$
$\displaystyle\left\\{\sum_{m=1}^{M}\left[\left(\sum_{n=1}^{N}\alpha_{mn}^{y}\sin\left(\frac{n\pi}{b}(y^{\prime}+\frac{b}{2})\right)+\sum_{n=0}^{N}\beta_{mn}^{y}\cos\left(\frac{n\pi}{b}(y^{\prime}+\frac{b}{2})\right)\right)\sin\left(\frac{m\pi}{a}(x^{\prime}+\frac{a}{2})\right)\right]\right\\}$
(11b)
$\displaystyle\left.e^{j(k_{x}x^{\prime}+k_{y}y^{\prime})}\mathrm{d}x^{\prime}\mathrm{d}y^{\prime}\right|$

##  IV Solving Integral Equation for Aperture Field

This section describes how to solve the integral equation (10), in this
electromagnetic inverse problem.

###  IV-A Modal Expansion on Aperture Field

Since directly solving the integral equation is difficult, we expand the
unknown integrand $\boldsymbol{E}_{\mathrm{ape}}$ using the modal expansion
method. Specifically, $\boldsymbol{E}_{\mathrm{ape}}$ satisfies the boundary
conditions (12), as can be seen from Fig. 1.

$\left\\{\begin{array}[]{l}E_{\mathrm{ape},x}|_{y=-\frac{b}{2},\frac{b}{2}}=0\\\
E_{\mathrm{ape},y}|_{x=-\frac{a}{2},\frac{a}{2}}=0\\\ \end{array}\right.$ (12)

Equation (11) convert the problem of solving
$\boldsymbol{E}_{\mathrm{ape}}(x,y)$ to solving the unknown complex modal
coefficients $\alpha_{mn}^{i}$ and $\beta_{mn}^{i}$, $(i=x,y)$. The selection
of the basis function is not unique. The trigonometric functions were chosen
because they are the harmonic functions of the waveguide harmonic equation, if
we regard the modal shown in Fig. 1 as a waveguide connected to the infinite
PEC plane. Due to the existence of the analytic expression, (11) can use an
analytical integration instead of numerical integration, which improves
computational efficiency. The analytical integration of (11) is given in
Appendix A.

###  IV-B Obtain Modal Coefficients by Using MOO

We obtain $\alpha_{mn}^{i}$ and $\beta_{mn}^{i}$ $(i=x,y)$ by using the multi-
objective optimization algorithm(MOO). Here we use the PSO toolbox, and define
(13) as the objective function of MOO,

$obj=1-\frac{\mathrm{cov}\left(\mathcal{N}\left(|\boldsymbol{f}_{\mathrm{tar}}(\theta,\phi)|\right),\mathcal{N}\left(|\boldsymbol{f}_{\mathrm{cal}}(\theta,\phi)|\right)\right)}{\sigma_{\mathcal{N}\left(|\boldsymbol{f}_{\mathrm{tar}}(\theta,\phi)|\right)}\cdot\sigma_{\mathcal{N}\left(|\boldsymbol{f}_{\mathrm{cal}}(\theta,\phi)|\right)}}$
(13)

where $\mathrm{cov}$ denotes the covariance of two vectors, the operator
$\mathcal{N}$ represents unfolding a two-dimensional matrix into a one-
dimensional vector, $|\boldsymbol{f}_{\mathrm{tar}}|$ is known from (6) and
(7), and $\boldsymbol{f}_{\mathrm{cal}}$ is calculated from the optimized
$\alpha_{mn}^{y}$ and $\beta_{mn}^{y}$. The choice of objective function is
not unique and can be selected according to the design requirements.

##  V Numerical Simulation

This section shows numerical experiments using the above method to implement a
specific $|\boldsymbol{e}_{\mathrm{tar}}|$. To illustrate the significance of
our method, we use it to an example breaking the cosine rule.

###  V-A Specially selected design targets

It is well known that when scanning to 60 degrees, the directivity of the 2D
array factor is half of that at broadside. This is the cosine rule in wide-
angle scanning[15].

Take a uniform amplitude broadside 16*16 dipole array with elements spaced at
$\lambda/{2}$ distributed in the $xoy$ plane as an example, its three-
dimensional radiation pattern is shown in Fig. 2(a), denoted as
$|\boldsymbol{e}_{\mathrm{tar},0^{\circ}}|$. The directivity of the dipole
array is 25.29 dBi and decreases as the scan angle increases, as shown by the
solid line in Fig. 3. The dashed line in Fig. 3 represents the results
according to the cosine rule.

(a)  (b)

(c)  (d)

Figure 2: Target radiation pattern scanning to different angles. (a) Scanning
to $(0^{\circ},0^{\circ})$;(b) Scanning to $(30^{\circ},0^{\circ})$; (c)
Scanning to $(45^{\circ},0^{\circ})$;(d) Scanning to $(60^{\circ},0^{\circ})$.
Figure 3: Comparison of directivity of the cosine rule, dipole array, and the
aperture field we designed.

First, we need a set of design targets. We obtain it by rotating the
coordinate system. And this set of target pattern has different beam direction
angle but the same directivity of 25.29 dBi, which does not conform to the
cosine rule, as shown in Fig. 2(b, c, d). Assume that the polarization is
linear and only vertical, set (14) according to case I in section III.

$\displaystyle|{e}_{\mathrm{tar},\text{scan\\_degree},v}|$
$\displaystyle=|\boldsymbol{e}_{\mathrm{tar},\text{scan\\_angle}}|$ (14a)
$\displaystyle{e}_{\mathrm{tar},\text{scan\\_angle},h}$ $\displaystyle=0$
(14b)

Once targets are determined, the first step is to calculate $|f_{x}|$ and
$|f_{y}|$ by using (8). Here, we only consider $|f_{y}|$, this is because the
calculated $|f_{x}|$ are approximately 0. Thus, the aperture field has only
$E_{\mathrm{ape},y}$, but no $E_{\mathrm{ape},x}$.

###  V-B Numerical Simulation Settings

Next step is solving $E_{\mathrm{ape},y}$ (more precisely, obtaining
$\alpha_{mn}^{y}$ and $\beta_{mn}^{y}$) from $|f_{y}|$ according to (15b) by
using MOO.

In the following numerical simulation, the size of the aperture is set to
$7.5\lambda\times 7.5\lambda$, as the 16*16 array. That is, the dimension of
aperture $a$ and $b$ are set to be $7.5\lambda$. The mode numbers $M$ and $N$
are set to be $8$ at small scanning angles and are set to be $16$ at more than
$45^{\circ}$. Fewer modes would introduce low design freedom, resulting in
failure to achieve targets. Too many modes increase the computational burden
and make the calculated aperture field vary dramatically. Such oscillation is
not expected and is difficult to achieve. There has to be a tradeoff, and an
acceptable approach is that the shortest modal wavelength is less than half of
the operating wavelength. This makes the oscillation easy to realize in
physics. $\theta$ and $\phi$ need to be divided into appropriate steps
according to the fluctuation of the target pattern. Here, $\theta$ from 0 to
$\pi/2$ and $\phi$ from 0 to $2\pi$ are divided into 45 and 180 parts.

###  V-C Numerical Simulation Results

The normalized amplitude and phase distribution of aperture field calculated
by this method are shown in Fig. 4. To verify the correctness, we use a script
to construct $E_{\mathrm{ape},y}$ and excite it in CST Microwave Studio. The
full-wave simulation results are shown in Fig.5, the directivity of the four
sub-pictures is 25.29 dBi, 25.29 dBi, 25.24 dBi, and 24.14 dBi, respectively.

Figure 4: The normalized amplitude distribution and phase distribution of
aperture field obtained by inverse design.

(a)  (b)

(c)  (d)

Figure 5:  The radiation pattern of aperture field designed by this method
with different scanning angles. (a) Aiming at $(0^{\circ},0^{\circ})$;(b)
Aiming at $(30^{\circ},0^{\circ})$; (c) Aiming at $(45^{\circ},0^{\circ})$;(d)
Aiming at $(60^{\circ},0^{\circ})$.

We have also done numerical simulation for other scanning angles and
summarized their directivity in the dot-dash line in Fig. 3. So far, we have
successfully inverse designed a series of aperture field, the directivity of
which break the cosine rule. This fully demonstrates the significance of the
method.

##  VI Conclusion

A method for designing aperture field from an arbitrary phaseless pattern is
proposed in this communication. This method has strong universality and can
deal with the non-linearity introduced by phaseless. We give the derivation of
the method, the analytical expression of the integral equation and the
parameter values. And successfully use this method to design a set of aperture
field that could break the cosine rule.

There is still much follow-up work to be done, such as what the bounds of
achievable radiation pattern complexity with finite apertures are, how to
physically implement the aperture field and so on.

##  Appendix A Analytical expression of Equation (11)

Equation (15) is the analytic expression of Equation (11).

$\displaystyle f_{x}(\theta,\phi)=\sum_{n=1}^{N}$
$\displaystyle\left(\sum_{m=1}^{M}\alpha_{mn}^{x}\prod_{i=1}^{3}{g_{1,i}}\left(m,n,\theta,\phi\right)\right.$
(15a)
$\displaystyle\left.+\sum_{m=0}^{M}\beta_{mn}^{x}\prod_{i=1}^{3}{g_{2,i}}\left(m,n,\theta,\phi\right)\right)$
$\displaystyle f_{y}(\theta,\phi)=\sum_{m=1}^{M}$
$\displaystyle\left(\sum_{n=1}^{N}\alpha_{mn}^{y}\prod_{i=1}^{3}{g_{1,i}}\left(m,n,\theta,\phi\right)\right.$
(15b)
$\displaystyle\left.+\sum_{n=0}^{N}\beta_{mn}^{y}\prod_{i=1}^{3}{g_{3,i}}\left(m,n,\theta,\phi\right)\right)$

where the $g_{j,i}$ are defined in (16).

$\displaystyle\chi_{1}=\sin\theta\sin\phi$ (16a)
$\displaystyle\chi_{2}=\sin\theta\cos\phi$ (16b)
$\displaystyle\begin{split}&g_{1,1}=\left(4abe^{-\frac{1}{2}jk\cos\phi\left(b\cos\theta+a\sin\theta\right)}\right)\\\
&\quad\times\Big{[}\left(4n^{2}\pi^{2}+b^{2}k^{2}\left(\cos\left(2\phi\right)+2\cos\left(2\theta\right)\sin^{2}\left(\phi\right)-1\right)\right)\Big{.}\\\
&\quad\Big{.}\times\left(a^{2}k^{2}\chi_{2}^{2}-m^{2}\pi^{2}\right)\Big{]}^{-1}\end{split}$
(16c) $\displaystyle
g_{1,2}=\left(-m\pi+e^{jka\chi_{2}}\left(m\pi\cos\left(m\pi\right)-jak\sin\left(m\pi\right)\chi_{2}\right)\right)$
(16d)
$\displaystyle\begin{split}&g_{1,3}=\left(e^{\frac{1}{2}jkb\cos\left(\phi+\theta\right)}n\pi\right.\\\
&\qquad\quad\left.+e^{\frac{1}{2}jkb\cos\left(\phi-\theta\right)}\left(-n\pi\cos\left(n\pi\right)+jbk\sin\left(n\pi\right)\chi_{1}\right)\right)\end{split}$
(16e) $\displaystyle g_{3,1}=g_{2,1}=g_{1,1}$ (16f)

$\displaystyle
g_{2,2}=\left(-ak\chi_{2}+e^{jka\chi_{2}}\left(-jm\pi\sin\left(m\pi\right)+ak\cos\left(m\pi\right)\chi_{2}\right)\right)$
(16g) $\displaystyle g_{2,3}=-jg_{1,3}$ (16h) $\displaystyle
g_{3,2}=\left(-jm\pi+e^{jka\chi_{2}}\left(jm\pi\cos\left(m\pi\right)+ak\sin\left(m\pi\right)\chi_{2}\right)\right)$
(16i)
$\displaystyle\begin{split}&g_{3,3}=\left(-bk\chi_{1}e^{\frac{1}{2}jkb\cos\left(\phi+\theta\right)}\right.\\\
&\qquad\quad\
\left.+e^{\frac{1}{2}jkb\cos\left(\phi-\theta\right)}\left(-jn\pi\sin\left(n\pi\right)+bk\cos\left(n\pi\right)\chi_{1}\right)\right)\end{split}$
(16j)

## References

  * [1] L. A. Bronckers, A. Roc’h, and A. B. Smolders, “A new design method for frequency- reconfigurable antennas using multiple tuning components,” _IEEE Transactions on Antennas and Propagation_ , vol. 67, no. 12, pp. 7285-7295, Dec. 2019.
  * [2] Y. Liu, Y. Jia, W. Zhang, Y. Wang, S. Gong, and G. Liao, “An integrated radiation and scattering performance design method of low-RCS patch antenna array with different antenna elements,” _IEEE Transactions on Antennas and Propagation_ , vol. 67, no. 9, pp. 6199-6204, Sep. 2019.
  * [3] X. Huang, Y. Liu, P. You, M. Zhang, and Q. H. Liu, “Fast linear array synthesis including coupling effects utilizing iterative FFT via least-squares active element pattern,” _IEEE Antennas and Wireless Propagation Letters_ , vol. 16, pp. 804-807, Aug. 2017.
  * [4] R. Bera, D. Mandal, R. Kar, and S. P. Ghoshal, “Non-uniform single-ring antenna array design using wavelet mutation based novel particle swarm optimization technique,” _Computers & Electrical Engineering_, vol. 61, pp. 151-172, Jul. 2017. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0045790617312569
  * [5] N. K. Nikolova, _Introduction to Microwave Imaging_ , ser. EuMA High Frequency Technologies Series. Cambridge University Press, Jul. 2017.
  * [6] S. Molesky, Z. Lin, A. Y. Piggott, W. Jin, J. Vucković, and A. W. Rodriguez, “Inverse design in nanophotonics,” _Nature Photonics_ , vol. 12, no. 11, pp. 659-670, Nov. 2018.
  * [7] Z. Wei and X. Chen, “Physics-inspired convolutional neural network for solving full-wave inverse scattering problems,” _IEEE Transactions on Antennas and Propagation_ , vol. 67, no. 9, pp. 6138-6148, Sep. 2019.
  * [8] M. Salucci, A. Gelmini, G. Oliveri, N. Anselmi, and A. Massa, “Synthesis of shaped beam reflectarrays with constrained geometry by exploiting nonradiating surface currents,” _IEEE Transactions on Antennas and Propagation_ , vol. 66, no. 11, pp. 5805-5817, Nov. 2018.
  * [9] T. Brown, C. Narendra, Y. Vahabzadeh, C. Caloz, and P. Mojabi, “On the use of electromagnetic inversion for metasurface design,” _IEEE Transactions on Antennas and Propagation_ , vol. 68, no. 3, pp. 1812-1824, Mar. 2020.
  * [10] C. Narendra, T. Brown, and P. Mojabi, “Gradient-based electromagnetic inversion for metasurface design using circuit models,” _IEEE Transactions on Antennas and Propagation_ , vol. 70, no. 3, pp. 2046-2058, Mar. 2022.
  * [11] R. Fletcher and C. M. Reeves, “Function minimization by conjugate gradients,” _The Computer Journal_ , vol. 7, no. 2, pp. 149-154, Jan. 1964\. [Online]. Available: https://doi.org/10.1093/comjnl/7.2.149
  * [12] A. Taaghol and T. Sarkar, “Near-field to near/far-field transformation for arbitrary near-field geometry, utilizing an equivalent magnetic current,” _IEEE Transactions on Electromagnetic Compatibility_ , vol. 38, no. 3, pp. 536–542, Aug 1996.
  * [13] M. Sørensen, H. Kajbaf, V. V. Khilkevich, L. Zhang, and D. Pommerenke, “Analysis of the effect on image quality of different scanning point selection methods in sparse esm,” _IEEE Transactions on Electromagnetic Compatibility_ , vol. 61, no. 6, pp. 1823–1831, Dec 2019.
  * [14] D. Yuan and X. Zhang, “An overview of numerical methods for the first kind Fredholm integral equation,” _SN Applied Sciences_ , vol. 1, no. 10, p. 1178, Oct. 2019. [Online]. Available: https://doi.org/10.1007/s42452-019-1228-3 http://link.springer.com/10.1007/s42452-019-1228-3
  * [15] N. Liu, S.-W. Qu, Y. Ma, and S. Yang, “Irregular phased array architecture with small directivity drop for wide-angle scanning application,” _IEEE Transactions on Antennas and Propagation_ , pp. 1–1, 2022.