#  Highly-confined and tunable plasmonics based on two-dimensional solid-state
defect lattices

Ali Ghorashi1 Nicholas Rivera1,4 Bowen Shi1,6 Ravishankar Sundararaman2
Efthimios Kaxiras3,4 John Joannopoulos1,5 Marin Soljačić1 1 Department of
Physics, Massachusetts Institute of Technology, Cambridge, MA 02139, USA  2
Rensselaer Polytechnic Institute, 110 8th Street, Troy, New York 12180, USA  3
School of Engineering and Applied Sciences, Harvard University, Cambridge, MA
02134, USA  4 Department of Physics, Harvard University, Cambridge, MA 02134,
USA  5 Institute for Soldier Nanotechnologies, Massachusetts Institute of
Technology, Cambridge, MA 02139, USA  6 State Key Laboratory for Mesoscopic
Physics and Department of Physics, Peking University, Beijing 100871, P. R.
China

(April 2023)

###### Abstract

Plasmons, collective excitations of electrons in solids, are associated with
strongly confined electromagnetic fields, with wavelengths far below the
wavelength of photons in free space. This strong confinement promises the
realization of optoelectronic devices that could bridge the size difference
between photonic and electronic devices. However, despite decades of research
in plasmonics, many applications remain limited by plasmonic losses, thus
motivating a search for new engineered plasmonic materials with lower losses.
A promising pathway for low-loss plasmonic materials is the engineering of
materials with flat and energetically isolated metallic bands, which can
strongly limit phonon-assisted optical losses, a major contributor to short
plasmonic lifetimes. Such electronic band structures may be created by
judiciously introducing an ordered lattice of defects in an insulating host
material. Here, we explore this approach, presenting several low-loss, highly-
confined, and tunable plasmonic materials based on arrays of carbon
substitutions in hexagonal boron nitride (hBN) monolayers. From our first-
principles calculations based on density functional theory (DFT), we find
plasmonic structures with mid-infrared plasmons featuring very high
confinements ($\lambda_{\text{vacuum}}/\lambda_{\text{plasmon}}$ exceeding
2000) and quality factors in excess of 1000. We provide a systematic
explanation of how crystal structure, electronic bandwidth, and many-body
effects affect the plasmonic dispersions and losses of these materials. The
results are thus of relevance to low-loss plasmon engineering in other flat
band systems.

The optoelectronic properties of two-dimensional (2D) materials are of major
interest due to the qualitatively different physics of electron-photon
interactions in reduced dimensions. Accordingly, the search for stable 2D
materials with specific optoelectronic properties has been a topic of intense
research in recent years Huang et al. (2022); Lončarić et al. (2018); Gjerding
et al. (2021). In particular, the expansion of the repertoire of stable 2D
materials has made strides in three directions:
(i) The tuning of geometric properties of van der Waals (vdW) bilayers and
trilayers, most notably by creating moiré patterns Lewandowski and Levitov
(2019);
(ii) the advent or prediction of 2D analogues of naturally occurring three
dimensional metals Sundararaman et al. (2020); Novoselov et al. (2005); and
(iii) the introduction of defects in common 2D materials Huang et al. (2022)
to induce desired optoelectronic properties.

A major impetus has been to enable collective excitations, such as plasmons
Boriskina et al. (2017), phonon-polaritons Dai et al. (2019); Rivera et al.
(2019) and exciton-polaritons Novko et al. (2021), with tailored dispersions,
high confinements, and/or low losses Khurgin (2015a), which would enable a
wide range of new light-matter interaction effects Rivera et al. (2016);
Rivera and Kaminer (2020). In the case of plasmons specifically, whose
electromagnetic fields can be confined far below the free-space wavelength of
photons, many of the promising applications envisioned for the field decades
ago are still hindered to this day by loss. Strongly-confined and low-loss
plasmonic excitations could lead to major advances for most envisioned
applications of plasmonics in fields spanning photovoltaics Atwater and Polman
(2010), spectroscopy Langer et al. (2019), biosensing Homola (2003), and
ultrahigh resolution lasers Noginov et al. (2009).

Our study is motivated by the goal of eliminating, or significantly reducing,
the ubiquitous losses intrinsic to plasmonic materials Khurgin (2015b);
Gjerding et al. (2017). In particular, we focus on 2D materials whose metallic
character is _induced_ by the presence of defects (an example of case (iii)
above). To avoid loss channels such as interband transitions, we restrict our
search to a host material with a large band gap, namely hexagonal boron
nitride (hBN) (6 eV bandgap, see Cassabois et al. (2016)), and defects that
produce moderately flat bands near the middle of the band gap 111A truly flat
band does not support plasmonic excitations and a highly dispersive band
increases intraband losses.. For the defects, we choose carbon atoms, since
their atomic size, similar to that of B and N, minimizes defect-induced
lattice strain. The substitutional defect structures we consider are denoted
as $C_{X}^{n\times n}$, indicating an $n\times n$ supercell of the primitive
unit cell of hBN, in which one atom labeled $X$ is replaced by a C atom (an
example of which is shown in Fig. 1 (b)). We study supercells with
$n=\sqrt{3},2,3,4$ and $X$=B or N, that is, a total of eight defect-containing
structures. Typically, larger supercells can host flatter bands, which are
more conducive to low-loss plasmonics.

Figure 1:  Electronic bandstructure of carbon substitutional defect lattices
in hexagonal boron nitride. (a) Electronic bands for $C_{B}^{n\times n}$ (left
column) and $C_{N}^{n\times n}$ (right column) structures. Energy is measured
from the Fermi level. For the six lattices that are insulating and
ferromagnetic, we shift the Fermi energy to half filling of the defect band
closest to the hBN midgap. We report the bandwidth, $W$, and conductivity,
$\sigma$, for each structure ($\sigma_{g}$ is graphene’s conductivity at 0.5
eV doping). The labels (1) and (3) denote intraband and interband Landau
damping processes, respectively. The labels (2) and (4) denote interband and
interband phonon-assisted damping processes, respectively. Black dashed lines
denote plasmons and pink dashed lines denote phonons. (b) The Wigner-Seitz
cells for the 4 structures we considered, within the largest, $4\times 4$
structure. The brown circle denotes the C defect at the center of the cell.
White circles denote N atoms for the B substitutional lattices and B atoms for
the N substitutional lattices. Green circles denote B atoms for the B
substitutional lattices and N atoms for the N substitutional lattices. The
color scheme is the same as in (a) with purple, orange, red, and light blue
corresponding to $\sqrt{3}\times\sqrt{3}$, $2\times 2$, $3\times 3$, and
$4\times 4$, respectively.

The eight structures and corresponding electronic band structures are shown in
Fig. 1. Our DFT calculations indicate that the relaxed structures with N
substitutional atoms are perfectly planar whereas those with B are slightly
buckled; see Supplementary Information (SI). While this is a minor structural
difference, we show later that it has significant impact on the magnitude of
the electron-phonon interaction. All of the structures, other than the two
$C^{\sqrt{3}\times\sqrt{3}}_{X}$ ones, are fully spin polarized, meaning that
they have to be doped in order to support plasmonic excitations. As previously
predicted Weston et al. (2018); Liu et al. (2022), we find that the structures
with larger supercells are ferromagnetic with spin gaps on the order of $1$
eV. To verify their ferromangetic order, we calculate the ground state of a
$2\times 2$ supercell of $C_{B}^{2\times 2}$ and find that the ferromagnetic
state is preferred, with magnetization $4\mu_{B}$. The structures we report
have bandwidths ranging from 3.3 eV ($C_{B}^{\sqrt{3}\times\sqrt{3}}$) to 0.13
eV($C_{B}^{4\times 4}$). This decrease in the bandwidth as lattice size is
increased is accompanied by a commensurate decrease in the Fermi velocity and
the onset of ferromagnetism (lifting of spin degeneracy), both of which
contribute to a lowering of the Drude conductivity as seen in Fig. 1 (a). This
leads to lower frequency plasmons for the larger defect supercells.

Limiting plasmonic losses has been an active area of research for several
decades Giuliani and Quinn (1983). In all proposed candidates for low-loss
plasmonics, losses are mitigated by suppressing the phase space for direct and
indirect (phonon-assisted) transitions into the electron-hole continuum. In
the case of surface plasmons, for instance, one early proposal for avoiding
losses was by engineering semiconducting superlattices that energetically
separate the surface plasmonic band from the electron-hole continuum Giuliani
and Quinn (1983). In the case of our structures that include substitutional
defects, a similar energetic separation should exist, as the structures with
larger periodicity host flat bands that are well separated from the hBN
valence and conduction bands. Accordingly, we calculate the transverse
magnetic (TM) polarized plasmonic dispersions and associated losses in the
proposed structures. For the $C_{B}^{\sqrt{3}\times\sqrt{3}}$ and
$C_{N}^{\sqrt{3}\times\sqrt{3}}$ structures, which are metallic at charge
neutrality, we calculate the plasmonic properties without imposing any changes
in band occupation. For the structures with larger supercells, we move the
Fermi level to half filling of the defect band closest to the hBN midgap and
calculate plasmonic properties within the rigid band approximation Lee and
Mahanti (2012), which neglects changes to the band structure due to doping (we
analyze the validity of this approximation below). Experimentally, this tuning
of the Fermi level could be implemented by gating or chemical doping.

To obtain the plasmonic dispersion for each case we evaluate the nonlocal,
frequency dependent conductivity. Energies and wavefunctions used to calculate
the conductivity are obtained from an ab initio tight binding model, derived
from Wannier interpolation of the DFT band structures Sundararaman et al.
(2020). Specifically, we calculate the plasmon dispersion through the poles of
the inverse dielectric function, $\varepsilon^{-1}(q,\omega)$, which at finite
temperature is given by Mahan (2013):

$\frac{1}{\varepsilon(q,i\omega_{n})}=1-\frac{V_{q}}{\Omega}\int_{0}^{\beta}e^{i\omega_{n}\tau}\braket{{\cal
T}\rho(q,\tau)\rho(-q,0)}d\tau,$ (1)

where $\Omega$ is the unit cell area, ${\cal T}$ is the imaginary time
ordering operator, $q$ is the wavevector and
$V_{q}=\frac{e^{2}}{2\epsilon_{0}q}$ is the Coulomb interaction in 2D (with
$\epsilon_{0}$ the vacuum permittivity), $\beta$ is the inverse temperature,
and $\rho(q,\tau)$ is the density operator in the Heisenberg representation.
The above expression yields the inverse dielectric function at a bosonic
Matsubara frequency, $\omega_{n}$. The retarded inverse dielectric function is
then calculated by analytically continuing
$i\omega_{n}\rightarrow\omega+i\delta$.

Electron-electron interactions are included through the Random Phase
Approximation (RPA) Mahan (2013), which gives the standard result (neglecting
local-field effects) Agarwal et al. (2014) 222In the present work, we
disregarded local field effects Adler (1962); Wiser (1963), which are likely
substantial only when $|q|=|q+G|$, where $G$ is a reciprocal lattice vector. :

$\varepsilon(q,\omega)=1-\frac{e^{2}}{2\varepsilon_{0}q\Omega}{\cal
F}(q,\omega)$ (2) ${\cal
F}(q,\omega)=\sum_{k}\frac{f(\epsilon^{n}_{k+q})-f(\epsilon^{m}_{k})}{\epsilon^{n}_{k+q}-\epsilon^{m}_{k}-\hbar\omega-i\delta}|\braket{k+q,n}{k,m}|^{2},$
(3)

where $\langle k+q,n|$, $|k,m\rangle$ denote the cell periodic components of
the Kohn-Sham eigenstates with corresponding energy eigenvalues
$\epsilon^{n}_{k+q},\epsilon^{m}_{k}$. $f(\epsilon)$ is the Fermi occupation,
with $\delta$ a small quantity to avoid singularities at the poles. We note
that in Eq. 3 the spin indices are subsumed into the band indices, $m,n$. We
calculate $\varepsilon(q,\omega)$ from Eq. 2 and Eq. 3 by first evaluating the
imaginary part of this function and subsequently exploiting the Kramers-Kronig
relations to compute the corresponding real part. We report plasmonic
dispersions and confinements in Fig. 2 (a) as compared to the most well-
established 2D plasmonic platform, graphene, at 0.5 eV doping from the Dirac
point (0.5 eV doping puts the plasmons in the mid-IR). The plasmons in the
proposed structures cover a frequency range from 0 to $\sim 1.3$ eV, with
plasmons below 1 eV (see SI) immune to interband and intraband losses. These
plasmons have confinements in the plotted wavevector range of up to $\sim 8$
times that of graphene (for the $C_{B}^{4\times 4}$ structure). We note also
that the generically small plasmonic group velocities observed for most of the
proposed structures is a consequence of interband screening, a nearly
universal phenomenon in 2D da Jornada et al. (2020). This is not the case for
the $C_{B}^{\sqrt{3}\times\sqrt{3}}$ and $C_{N}^{\sqrt{3}\times\sqrt{3}}$
structures, as in those two cases the plasmon approaches the interband
continuum at large wavevectors $q$.

To take into account the effect of phonons on the imaginary part of the
dielectric function or, equivalently, the real part of the conductivity, we
follow the prescription given in Stauber and Peres (2008); Mahan (2013); Allen
(2015) and evaluate the current-current response in the presence of electron-
phonon interactions. We then translate this into a decay time in the Drude
limit, using the relation:

$\sigma(\omega)=\frac{e^{2}v_{F}^{2}g(\epsilon_{F})\tau(\omega)}{2\Omega(1-i\omega\tau(\omega))},$
(4)

where $v_{F}$ is the Fermi velocity, $g(\epsilon_{F})$ is the density of
states per unit cell at the Fermi level and $\tau(\omega)$ is the transport
decay time (see SI). Through a Feynman diagrammatic expansion of the
conductivity, the decay time can be rewritten in a style reminiscent of
Fermi’s golden rule Allen (1971); Brown et al. (2016), giving an expression
that is essentially a Fermi surface average of the frequency-dependent carrier
decay rate:

$\tau^{-1}(\omega)=\frac{2\pi}{N_{k}N_{k^{\prime}}\hbar^{2}\omega
g(\epsilon_{F})}\sum_{k,k^{\prime}j\pm}|g_{k,k^{\prime}}^{j}|^{2}\Big{(}N^{j,\mp}_{k-k^{\prime}}f_{k}-\\\
N_{k-k^{\prime}}^{j,\pm}f_{k^{\prime}}\pm
f_{k}f_{k^{\prime}}\Big{)}\delta(\epsilon_{k}+\hbar\omega\pm\hbar\omega^{j}_{k-k^{\prime}}-\epsilon_{k^{\prime}})\Bigg{(}1-\frac{v_{k}\cdot
v_{k^{\prime}}}{|v_{k}||v_{k^{\prime}}|}\Bigg{)},$ (5)

where $v_{k}$, $v_{k^{\prime}}$ are the electronic velocities (corresponding
to the intraband momentum matrix elements) at wavevectors $k$ and
$k^{\prime}$, respectively, $\omega_{q}^{j}$ are the phonon frequencies of
branch $j$ at wavevector $q$, with corresponding $N^{j}_{q}$ Bose occupation
factors. For convenience, in the above equation we have defined the quantities

$N_{q}^{j,\pm}\equiv\frac{1}{2}+N_{q}^{j}\pm\frac{1}{2},$

where the plus sign corresponds to phonon emission and the minus sign
corresponds to phonon absorption. We sum over phonon bands, indexed by $j$,
but we include only the defect-related electronic band in the evaluation of
the decay rate from Eq. 5, which is an exact expression in the frequency
regime of interest (0-1 eV). We note that the combination of velocity factors
in our expression, in the language of Feynman diagrams, is a consequence of a
vertex correction to the polarization bubble. In addition, we note that the
combination of occupation factors in Eq. 5 correctly takes into account
reverse processes, as may be seen by the fact that the same combination of
occupation factors appears in the time derivative of the phonon occupation
factor, which vanishes at equilibrium due to detailed balance Pines and
Schrieffer (1962). Lastly, we mention that Eq. 5 simplifies to the Fermi
surface averaged electron-phonon decay rate in the limit
$\hbar\omega\rightarrow 0$.

Figure 2: Plasmonic dispersions, confinements and losses of carbon
substitutional lattices in hexagonal boron nitride. (a) Plasmonic dispersions
(solid) and confinements (dashed) lines compared to graphene. Wavevectors were
chosen to lie on the $\Gamma-M$ direction. (b) Decay times for $C_{B}^{n\times
n}$ ($C_{N}^{n\times n}$) in solid (dashed) lines and graphene at 0.5 eV
doping (dashed grey line). The point at which a solid (dashed) line turns into
a dotted line indicates the highest plasmon frequency. The termination points
of the dotted lines indicates the highest theoretical plasmon frequency (see
discussion surrounding Eq. 6).

In Fig. 2 (b) we present plasmonic decay times, $\tau(\omega)$. The
$C_{N}^{\sqrt{3}\times\sqrt{3}}$ and $C_{N}^{2\times 2}$lattices yield better
plasmonic lifetimes than graphene at high frequencies, with
$C_{N}^{\sqrt{3}\times\sqrt{3}}$ plasmons having lifetimes about four times
those of graphene plasmons at $\sim 1$ eV. The fact that the decay times are
lower for the B substitutional structures is due to enhancement of the
electron-phonon interaction through the aformentioned buckling of the $C_{B}$
structures. We verified this by explicitly comparing the decay times for
buckled and unbuckled (non-relaxed) $C_{B}$ structures (see SI). Though two
$C_{N}$ lattices yield a small frequency range at which plasmonic quality is
improved, it is clear that in general the proposed structures have shorter
plasmonic decay times than graphene at most frequencies. This may be
attributed to the fact that in the carrier decay rate, the density of
electronic states effectively shows up twice in the numerator and only once in
the denominator (SI). Thus, as a flat band hosts a tightly confined (in
frequency) region with a high density of states, our observation of a small
decay time at low frequencies is to be expected. In addition, the fact that
flatter bands in general host lower frequency plasmons results in plasmons
being pushed into the regime of high loss. This fact, however, offers insight
that can lead to a solution of the loss problem, posed as follows: If flatter
bands yield a frequency range of low plasmonic loss that is above the actual
plasmon frequency, how can the plasmon frequency be pushed into this low-loss
region? We investigate this question by determining, for each structure, their
maximum theoretical plasmon frequency, defined, for each $q$, as the highest
possible frequency, $\omega$, for which the dielectric function
$\varepsilon(q,\omega)$, defined in Eq. 2, vanishes if we introduce a factor
of $\alpha$, a positive constant, multiplying the quantity ${\cal
F}(q,\omega)$ defined in Eq. 3. Note that in the case of an infinite
superlattice of structures stacked in the direction perpendicular to the
planes ($z$-axis), this parameter takes the form

$\alpha=\frac{\sinh(ql)}{\cosh(ql)-1}$ (6)

where $l$ is the spacing of layers Giuliani and Quinn (1983). In Fig. 2 (b) we
denote by dotted lines the regime beyond the plasmon frequency and below the
maximum theoretical plasmon frequency. Interestingly, for the structures with
larger periodicity, where interband polarization plays a dominant role, the
plasmon frequency is not enhanced significantly in the stacked superlattices.
Only for the case of $C_{N}^{2\times 2}$ and $C_{B}^{2\times 2}$ stacked
superlattices is the plasmon dispersion enhanced significantly. To determine
more accurately how much of an enhancement is possible, we calculate the
plasmon dispersion of a $C_{N}^{2\times 2}$ superlattice with $l=1$ nm and
determine that the plasmon frequency could reach $\approx 0.9$eV with a
quality factor of 1100. Note that this result comes with a caveat: The reason
that a large value for $l$ is needed to obtain this effect, is that simply
stacking the doped structures at their equilibrium interlayer distance
$l\approx 0.3$ nm drastically changes the electronic structure. Thus, to
actually achieve this effect experimentally, the doped layers need to be
separated by a spacer of low dielectric constant, such as pristine hBN. A
systematic study of how the electronic structure would change in these many
layered systems is beyond the scope of the present work.

Figure 3: Effect of doping on the spin-dependent bandstructures and plasmonic
dispersion. (a) Band structures near the Fermi level as a function of doping
for $C^{2\times 2}_{B}$ (left column) and $C_{B}^{4\times 4}$ (right column).
(b) Accuracy of our spin splitting model given by Eq. 7. Each blue dot
represents a DFT calculation for which we plot our predicted spin splitting
against the DFT result. The proximity of each point to the $y=x$ line (shown
in black) indicates the validity of the perturbation theory model. (c)
Plasmonic dispersion for $C_{B}^{2\times 2}$ at various values of doping
calculated in the rigid band approximation (dashed lines) and from the
explicit charged DFT calculations (solid lines)

We turn next to the validity of the rigid band approximation. Interestingly,
we find that electron and hole doping in the structures we consider has the
effect of tuning the spin-gap at the Fermi level (see Fig. 3 (a)), which is at
odds with the bands being rigidly doped. We attribute this result to changes
in the exchange potential with respect to doping. In particular, the exchange
potential in DFT is a functional of the spin-resolved electron densities. From
first order perturbation theory in the local density approximation (LDA), the
doping results in a spin gap given by (see SI for details):

$\begin{split}\Delta(\rho_{\uparrow}(r),\rho_{\downarrow}(r))=\Big{(}\frac{6}{\pi}\Big{)}^{1/3}\int(\rho_{\uparrow}^{1/3}(r)-\rho_{\downarrow}^{1/3}(r))\times\\\
\frac{(|\psi_{\uparrow}(r)|^{2}+|\psi_{\downarrow}(r)|^{2})}{2}d^{3}r,\end{split}$
(7)

where $\rho_{\uparrow},\rho_{\downarrow}$ and
$\psi_{\uparrow},\psi_{\downarrow}$ are the spin up (down) densities and
wavefunctions, respectively. We evaluate the reliability of this model by
performing DFT calculations on the doped structures as shown in Fig. 3 (b).
For each material, we use three functionals: PBE Perdew et al. (1996), Slater
LDA without correlation Dirac (1930) and LDA with correlation Perdew et al.
(1981). For each value of doping, we calculate the spin splitting from DFT at
the $\Gamma$ point (center of the Brillouin Zone in reciprocal space). We also
use the self consistent density and the Kohn-Sham wavefunctions for the defect
bands to calculate the predicted spin splitting value given by Eq. 7. In Fig.
3(b), we plot our prediction against the exact DFT values. The plots in Fig.
3(b) are obtained by using the PBE-GGA exchange correlation functional for the
structures $C_{B}^{2\times 2}$ and $C_{B}^{4\times 4}$. As is evident, the fit
of our model to the numerical DFT results is very accurate. We present similar
plots for the other two exchange-correlation functions and for all other
lattices in the SI. Note that the accurate fit we obtain in Fig. 3(b)
indicates that the effect of correlation and gradient terms in the Kohn-Sham
Hamiltonian largely cancel out. Lastly, we note that a purely analytic
estimate of the spin splitting may be carried out by decomposing the densities
in terms of $1s$, $2s$, and $2p$ orbitals. We perform this calculation in the
SI and obtain good agreement between this purely analytical calculation and
our first principles results.

In Fig. 3 (c) we show plasmonic dispersions as a function of doping. We find
that the explicitly doped plasmons differ qualitatively from those obtained
through the rigid band approximation. At first this may seem like a
consequence of the changing spin gap; however, we systematically investigated
this issue and found that the major factor is actually a significant
enhancement in interband wavefunction overlaps. Thus, the discrepancies in
Fig. 3 (c) are primarily a consequence of the neutralizing charge background.
To better understand how doping affects the plasmonic bands, other physically
motivated scenarios for doping, such as through lithium intercalation Profeta
et al. (2012), need to be investigated.

In conclusion, we have introduced a set of candidate 2D materials with novel
plasmonic properties, consisting of doped hBN through C substitution at either
B or N sites. In particular, we predict these structures to host plasmons with
confinements up to eight times the maximum achievable in graphene, with decay
times that can also surpass that of graphene, for frequency ranges exceeding
$\sim 0.4$ eV. We have shown that stacking these structures could yield
exceptionally high quality factors. We expect that imperfections in the
periodicity of these materials to not have qualitative consequences as long as
the density of defects is similar to the superlattices we investigated.
However, this issue and the effect of impurities Peres et al. (2008) is a
topic for future investigation. As most of the proposed structures have low
Fermi velocities, it is possible that their electron-electron interactions
necessitate an approach beyond RPA. Investigation of the validity of the
random phase approximation in treating these flat band systems and their
collective excitations warrants further investigation. Similar work which
explored beyond RPA diagrams in the case of graphene Gangadharaiah et al.
(2008), has already been done but not applied to defect structures of the type
considered here. The effect of the electron-plasmon interaction on the carrier
lifetimes should also be considered for a better assessment of the Drude decay
time Polini et al. (2008). In addition, even in the absence of doping, exciton
polaritons may exist in the proposed structures Henriques et al. (2022), and
needs to be further investigated. We note, in closing, that there has recently
been a flurry of work in creating databases of 2D materials Haastrup et al.
(2018); Gjerding et al. (2021). While a “blind” enumeration of the plasmonic
properties of all tabulated materials would be an overly demanding task, the
results presented here suggest a simpler approach. Namely, we expect that
filtering the available databases for materials with isolated flat bands at
the Fermi level and high structural and thermal stability would be a first
step in identifying the most promising candidates for low plasmonic losses.

### Computational Methods

Density functional theory calculations were carried out with the use of the
JDFTx package Sundararaman et al. (2017) with norm-conserving pseudopotentials
Schlipf and Gygi (2015) and Coulomb truncation for 2D materials Sundararaman
and Arias (2013). The ab initio tight binding models used to calculate
plasmonic properties were obtained by mapping onto a maximally localized
Wannier basis Souza et al. (2001). Ground state properties were calculated
with a Fermi-level smearing of $10^{-5}$ Hartree, corresponding to a
temperature of approximately 3 K. A higher smearing was used for the two
metallic lattices due to the existence of a Fermi surface. Unless otherwise
stated, the PBE exchange-correlation functional was used Perdew et al. (1996).
Comparisons between different functionals implemented in our modeling of
doping dependent spin splitting additionally made use of pure Slater exchange
Dirac (1930) and the Perdew-Zunger local density approximation Perdew et al.
(1981). Ground state properties were found through either the self-consistent
field method Kresse and Furthmüller (1996) or through electronic minimization
Freysoldt et al. (2009). We show plots indicating convergence of the decay
time results with respect to the phonon supercell size in SI.

### Acknowledgments

This material is based upon work supported by the Air Force Office of
Scientific Research under the award number FA9550-21-1-0299, as well as in
part by the U. S. Army Research Office through the Institute for Soldier
Nanotechnologies at MIT, under Collaborative Agreement Number
W911NF-18-2-0048. A.G. thanks the National Science Foundation Graduate
Research Fellowship for financial support during the preparation of this
paper. N.R. acknowledges the support of a Junior Fellowship from the Harvard
Society of Fellows, as well as earlier support from a Computational Science
Graduate Fellowship of the Department of Energy (DE-FG02-97ER25308), and a
Dean’s Fellowship from the MIT School of Science. E.K. is supported in part by
an Army Research Office grant under Cooperative Agreement Number
W911NF-21-2-0147. The authors would also like to thank Ali Fahimniya, Cyprian
Lewandowski, Thomas Christensen, Marinko Jablan, and Jennifer Coulter for
useful discussions.

##  I Supplementary Information for: Highly-confined and tunable plasmonics
based on two-dimensional solid-state defect lattices

##  Appendix A Summary of Main Equations and Results

Herein, we derive the main equations used in the main text and provide
additional sanity checks, further motivating our results. First, we rigorously
prove our RPA formulation for the dielectric function of our spinful systems.
Next, we derive the expressions used to calculate the plasmonic dispersions,
Eq. S9, and confinements. We next compare our decay time formula with an
alternative expression , Eq. S12, reported in the literature Kumar et al.
(2022). Lastly, we lay out our formulation of doping dependent spin splitting,
culminating in Eq. S17 in the local density approximation. We also include a
figure depicting the convergence of our decay time calculations with respect
to phonon supercell size (Fig. S8) and show how buckling of the B substitution
lattices induces enhanced electron-phonon coupling in Fig. S7.

##  Appendix B The Random Phase Approximation for Spinful Systems

As the RPA is usually derived in a spin agnostic manner, we include a
derivation below that retains the full spin structure. We remind the reader of
the exact inverse dielectric function:

$\small\epsilon^{-1}(q,i\omega_{n})=1-\frac{V_{q}}{\Omega}\sum_{k,k^{\prime}}\int_{0}^{\beta}\braket{T(c_{k+q,\uparrow}^{\dagger}(\tau)c_{k,\uparrow}(\tau)+c_{k+q,\downarrow}^{\dagger}(\tau)c_{k,\downarrow}(\tau))(c_{k^{\prime}-q,\uparrow}^{\dagger}(\tau^{\prime})c_{k^{\prime},\uparrow}(\tau^{\prime})+c_{k^{\prime}-q,\downarrow}^{\dagger}(\tau^{\prime})c_{k^{\prime},\downarrow}(\tau^{\prime}))}e^{i\omega_{n}(\tau-\tau^{\prime})}d(\tau-\tau^{\prime}),$
(S1)

where $V_{q}=\frac{e^{2}}{2\epsilon_{0}q}$ is the Coulomb interaction in 2D, q
is the wavevector, $\Omega$ is the unit cell area, $i\omega_{n}$ is a bosonic
Matsubara frequency (we take $i\omega_{n}\rightarrow\omega+i\delta$ to get the
retarded inverse dielectric function), and the $c_{k,(\downarrow,\uparrow)}$,
$c^{\dagger}_{k,(\downarrow,\uparrow)}$ are spin-resolved Bloch state
annihilation/creation operators. In the above, all operators are in the
Heisenberg picture. We may write the above in terms of diagonal and off-
diagonal (in spin indices) (interacting) density-density correlation
functions:

$\epsilon^{-1}(q,\omega)=1+\frac{V_{q}}{\Omega}(\chi_{\uparrow\uparrow}(q,\omega)+\chi_{\downarrow\downarrow}(q,\omega)+\chi_{\downarrow\uparrow}(q,\omega)+\chi_{\uparrow\downarrow}(q,\omega))$
(S2)

We write the RPA Dyson series for the diagonal up-up correlation function as
(subsuming the area of the sample into the correlation function):

$\chi_{\uparrow\uparrow}(q,\omega)=\chi_{\uparrow\uparrow}^{0}(q,\omega)+V_{q}(\chi_{\uparrow\uparrow}^{0}(q,\omega))^{2}+V_{q}^{2}(\chi_{\uparrow\uparrow}^{0}(q,\omega))^{2}(\chi_{\uparrow\uparrow}^{0}(q,\omega)+\chi_{\downarrow\downarrow}^{0}(q,\omega))+...$
(S3)

At every order (after the zeroth order), we have two factors of the non-
interacting up-up correlation function and one factor of
$(\chi_{\uparrow\uparrow}^{0}(q,\omega)+\chi_{\downarrow\downarrow}^{0}(q,\omega))$.
In the language of Feynman diagrams, this simply means that the end points of
the RPA expansion have restricted spin indices whereas the intermediate
vertices do not. The down-down correlation function may be similarly expanded
as a Dyson series, and we do not reproduce it here. The up-down correlation
function may be written as:

$\chi_{\uparrow\downarrow}(q,\omega)=V_{q}\chi^{0}_{\uparrow\uparrow}(q,\omega)\chi^{0}_{\downarrow\downarrow}(q,\omega)+V_{q}^{2}\chi^{0}_{\uparrow\uparrow}(q,\omega)\chi^{0}_{\downarrow\downarrow}(q,\omega)(\chi^{0}_{\uparrow\uparrow}(q,\omega)+\chi^{0}_{\downarrow\downarrow}(q,\omega))+...$
(S4)

By the same logic as above, this Dyson series may be motivated by the fact
that the endpoint vertices of the Feynman diagrammatic expansion have
restricted spin indices (up and down, respectively). We may now readily sum up
the Dyson series:

$\chi_{\uparrow\uparrow}(q,\omega)=\chi^{0}_{\uparrow\uparrow}(q,\omega)+\frac{V_{q}\chi^{0}_{\uparrow\uparrow}(q,\omega)^{2}}{1-V_{q}(\chi^{0}_{\uparrow\uparrow}(q,\omega)+\chi^{0}_{\downarrow\downarrow}(q,\omega))},\chi_{\uparrow\downarrow}(q,\omega)=\frac{V_{q}\chi^{0}_{\uparrow\uparrow}\chi^{0}_{\downarrow\downarrow}}{1-V_{q}(\chi^{0}_{\uparrow\uparrow}(q,\omega)+\chi^{0}_{\downarrow\downarrow}(q,\omega))}$
(S5)

Therefore, we have:

$\epsilon^{-1}(q,\omega)=1+V_{q}\frac{\chi^{0}_{\uparrow\uparrow}(q,\omega)(1-V_{q}\chi^{0}_{\downarrow\downarrow}(q,\omega))+\chi^{0}_{\downarrow\downarrow}(q,\omega)(1-V_{q}\chi^{0}_{\uparrow\uparrow}(q,\omega))+2\chi^{0}_{\uparrow\uparrow}(q,\omega))\chi_{\downarrow\downarrow}(q,\omega))}{1-V_{q}(\chi_{\uparrow\uparrow}(q,\omega)+\chi_{\downarrow\downarrow}(q,\omega))}$
(S6)

Which gives us for the dielectric function,

$\epsilon(q,\omega)=1-V_{q}(\chi^{0}_{\uparrow\uparrow}(q,\omega)+\chi^{0}_{\downarrow\downarrow}(q,\omega))$
(S7)

This is the same result as derived in Agarwal et al. (2014), although in that
work the Dyson series is written in matrix form:

$\chi^{-1}_{ij}=(\chi^{0})^{-1}_{ij}-V_{q}$ (S8)

##  Appendix C Plasmonic Dispersions

We calculate TM polarized plasmons with mode profiles given as follows for the
electric field (we choose the propagation direction of the plasmon to be in
the x direction- other propagation directions follow trivially):

$E(x,y,z)=\begin{cases}e^{iqx-Qz}(e_{x}+\frac{iq}{Q}e_{z}),&\text{if }z>0\\\
e^{iqx+Qz}(e_{x}-\frac{iq}{Q}e_{z}),&\text{if }z<0\end{cases}$

To establish a relation between the frequency, $\omega$ and the wavevector,
$q$, we calculate the current, given by
$J_{x}(q,\omega)=\sigma(q,\omega)E_{x}(q,\omega)$. This current carries a
corresponding charge density
$\rho(q,\omega)=\frac{q}{\omega}\sigma(q,\omega)E_{x}(q,\omega)$. Using
Gauss’s theorem, we obtain:

$\frac{2iq}{Q}=\frac{q}{\omega\epsilon_{0}}\sigma(q,\omega)\rightarrow\frac{2i}{\sqrt{q^{2}-\omega^{2}/c^{2}}}=\frac{\sigma(q,\omega)}{\omega\epsilon_{0}}$
(S9)

In principle, one must calculate the full non-local and frequency dependent
conductivity- which is the way by which we obtain the reported plasmonic
dispersions in the main text. To give an order of magnitude estimate of the
plasmonic dispersion, however, we may calculate the Drude conductivity, given
by:

$\frac{ig_{s}e^{2}}{N_{k}\omega\Omega}\sum_{k,n}\delta(\epsilon_{k}^{n}-\mu)|v^{x}_{n,k}|^{2},$
(S10)

where $v_{n,k}^{i}$ is the velocity of band $n$ at wavevector $k$ in direction
$i$, $\Omega$ is the unit cell area, $g_{s}$ is the spin degeneracy- which is
1 in the case of our defect lattices larger than $\sqrt{3}\times\sqrt{3}$ and
2 in the case of the two $\sqrt{3}\times\sqrt{3}$ lattices- and $N_{k}$ is the
number of k points sampled in the Brillouin zone. We find that the local,
Drude, conductivities of $C_{B}^{\sqrt{3}\times\sqrt{3}}$, $C_{B}^{2\times
2}$, $C_{B}^{3\times 3}$, $C_{B}^{4\times 4}$ are, respectively, $278\%$,
$86\%$, $18.5\%$, $7.5\%$ that of graphene at $0.5$ eV doping, and the Drude
conductivites for$C_{N}^{\sqrt{3}\times\sqrt{3}}$, $C_{N}^{2\times 2}$,
$C_{N}^{3\times 3}$, $C_{N}^{4\times 4}$ are, respectively, $187\%$, $75\%$,
$13\%$, $8.7\%$ that of graphene at 0.5 eV doping. Physically, this trend
reflects the fact that larger superlattices host progressively flatter bands-
with lower Fermi velocities- at the chemical potential. We note, also, that
our lattices have a background effective dielectric contrast given by:

$\epsilon(q)\approx 1+c\alpha_{BN}q$ (S11)

where $\alpha_{BN}$ is 12.69 in atomic units Lončarić et al. (2018) and $c$
ranges from $1.6$ to $3.5$ depending on the particular superlattice. We
confirm the value given for $\alpha_{BN}$ through a DFT calculation of pure
hBN and obtain the range of values for $c$ by calculating how interband
transitions alter the plasmonic dispersion in Fig. S5.

##  Appendix D Decay Times

The lowest order of plasmon decay in our defect lattices is through the
electron-phonon interaction. This is because all six plasmon dispersions are
Landau undamped, as verified in Fig. S4. In order to solidify the validity of
the decay time formula presented in the main text, we compare decay times
obtained through our formula with those obtained through:

$\tau^{-1}_{FS}(\omega)=\frac{1}{\hbar
g(\epsilon_{F})}\frac{2\pi}{N_{k^{\prime}}N_{k}}\sum_{k,k^{\prime}}|g_{k,k^{\prime}}^{n,m,\alpha}|^{2}\Big{(}1-\frac{v_{k}^{n}.v_{k^{\prime}}^{m}}{|v_{k}^{n}||v_{k^{\prime}}^{m}|}\Big{)}\delta(\epsilon_{k}^{n}-\epsilon_{F})\delta(\epsilon_{k^{\prime}}^{m}-\epsilon_{F})\sum_{\pm}\frac{\pm
b(\omega\pm\omega_{q}^{\alpha})}{(e^{\pm\beta\hbar\omega_{q}^{\alpha}}-1)b(\omega)}$
(S12)

Where $\beta(\omega)=\frac{\beta\hbar\omega}{1-e^{-\beta\hbar\omega}}$,
$g_{k,k^{\prime}}^{n,m,\alpha}$ is the electron-phonon matrix element
corresponding to electronic bands $n$ and $m$ at wavevectors $k$ and
$k^{\prime}$, respectively, and a phonon in branch $\alpha$ at wavevector
$k^{\prime}-k$ (of frequency $\omega_{k^{\prime}-k}^{\alpha}$).
$\epsilon_{k}^{n}$, $\epsilon_{k^{\prime}}^{m}$ are the electronic energies.
This formula has been used previously Kumar et al. (2022) to calculate the
frequency dependent carrier relaxation time. We compare decay times obtained
through our equation with those obtained through Eq. S12 in Fig. S1 and obtain
good agreement. The discrepancy for argentene at high frequencies is
attributed to additional interband transitions that are not captured in Eq.
S12.

Figure S1: Comparison of different decay times. Dashed lines are obtained
through the formula for $\tau(\omega)$ given in the main text, while the solid
lines are given by the alternative decay time formula $\tau_{FS}(\omega)$ in
Eq. S12. The discrepancy at high frequencies for Argentene is attributed to an
interband transition that invalidates the Fermi surface approximation.

A natural question to ask is how the above decay time is related to the
plasmonic decay time as defined in Wunsch et al. (2006). To make this
connection, we write the plasmonic decay rate as Wunsch et al. (2006):

$\frac{\Re(\sigma(\omega)/\omega)}{\frac{\partial\Im(\sigma(\omega)/\omega}{\partial\omega}}=1/2\tau(\omega)$
(S13)

##  Appendix E Doping Dependent Spin Splitting Model

We start with the expression for the exchange energy in the local density
approximation (LDA)

$E_{x}[\rho_{\uparrow},\rho_{\downarrow}]=\frac{1}{2}\Big{(}E_{x}[2\rho_{\uparrow}]+E_{x}[2\rho_{\downarrow}]\Big{)},$
(S14)

where $\rho_{\uparrow},\rho_{\downarrow}$ are the spin-resolved electronic
densities. We write the total electronic density as $\rho$ in the following.
In the absence of spin polarization, we have
$\rho_{\uparrow}=\rho_{\downarrow}=\rho/2$, and we obtain the same exchange
energy as when we neglect the existence of the two spin species. To determine
how the Kohn-Sham eigenvalues change- as a function of doping- we must examine
the exchange part of the Kohn-Sham potential. The matrix of exchange
potentials is given by Giustino (2014):

$V_{x}^{\alpha\beta}=\frac{\delta
E_{x}}{\delta\rho_{\alpha\beta}}=\delta_{\alpha\beta}\frac{\delta
E_{x}}{\delta\rho_{\alpha}},$ (S15)

Where we define:

$\rho_{\alpha\beta}=\sum_{i}\psi_{i,\alpha}^{*}\psi_{i,\beta}$ (S16)

Where $i$ indexes all non-spin quantum numbers (such as band and crystal
momentum) and $\alpha,\beta$ are spin indices (note also that, by definition,
$\rho_{\uparrow\uparrow}=\rho_{\uparrow}$ and
$\rho_{\downarrow\downarrow}=\rho_{\downarrow}$. By first order perturbation
theory, the spin splitting is given by (in atomic units):

$\Delta=\int\Big{(}V^{\uparrow}(r)-V^{\downarrow}(r)\Big{)}\frac{|\psi_{\uparrow}(r)|^{2}+|\psi_{\downarrow}(r)|^{2}}{2}d^{3}r=\Big{(}\frac{3}{\pi}\Big{)}^{1/3}2^{1/3}\int(\rho_{\uparrow}^{1/3}(r)-\rho_{\downarrow}^{1/3}(r))\frac{|\psi_{\uparrow}(r)|^{2}+|\psi_{\downarrow}(r)|^{2}}{2}d^{3}r$
(S17)

We first provide a fully analytical model for the spin splitting by using the
fact that the electronic density interacting with the 2pz orbital at the Fermi
level may be well described by a sum of densities contributed by the 1s, 2s
and the 2pz orbital itself.

$\Big{(}\frac{6}{\pi}\Big{)}^{1/3}\int\Big{[}\Big{(}|\psi_{1s}(r)|^{2}+|\psi_{2s}|^{2}+(1-n)|\psi_{2pz}|^{2}\Big{)}^{1/3}-\Big{(}|\psi_{1s}(r)|^{2}+|\psi_{2s}|^{2}\Big{)}^{1/3}\Big{]}|\psi_{2pz}|^{2}d^{3}r,$
(S18)

where n is the absolute value of the excess charge. We numerically integrate
this using the Hydrogenic wavefunctions for the 1s, 2s, and 2pz orbitals and
obtain good agreement with our DFT results as shown in Fig. S2.

In Fig. S3, we plot, for all six defect lattices, the result of Eq. S17
against the DFT spin splitting result. We do this for three exchange-
correlation functionals. Unsurpringly, LDA without correlation has the best
fit with our model. This is motivated by the fact that our perturbation theory
model uses the LDA exchange and neglects the effect of correlation. What is
more surprising, however, is that the more accurate PBE exchange also yields
good agreement with our model.

##  Appendix F Hubbard Model For Spin Splitting

We next show how the spin splitting may be predicted from the semianalytic
evaluation of the Hubbard $U$ parameter. Explicitly, the Hubbard $U$ is given
by:

$U=\frac{e^{2}}{4\pi\epsilon_{0}}\int|\psi_{pz}(r_{1})|^{2}\frac{1}{|r_{1}-r_{2}|}|\psi_{pz}(r_{2})|^{2}d^{3}r_{1}d^{3}r_{2}$
(S19)

However, this model overestimates the actual $U$ parameter. This is because we
have not taken into account screening. To take screening into account, we
consider the case where we have full spin polarization, and, therefore, an
insulator. In this case, the screened Coulomb interaction is given by (using
Eq. S11):

$W(\rho,\rho^{\prime},z,z^{\prime})=\frac{1}{(2\pi)^{2}}\int\frac{e^{2}e^{-iq|\rho-\rho^{\prime}|\cos(\theta)}e^{-q|z-z^{\prime}|}}{2\epsilon_{0}q(1+\alpha
q)}qdqd\theta$ (S20)

We perform the integral over $\theta$, giving us:

$\frac{e^{2}}{4\pi\epsilon_{0}}\int\frac{J_{0}(q|\rho-\rho^{\prime}|)e^{-q|z-z^{\prime}|}}{(1+\alpha
q)}dq$ (S21)

Therefore, the screened version of equation S19 is:

$\frac{e^{2}}{4\pi\epsilon_{0}}\int|\psi_{pz}(r_{1})|^{2}\frac{J_{0}(q|\rho-\rho^{\prime}|)}{1+\alpha
q}e^{-q|z-z^{\prime}|}|\psi_{pz}(r_{2})|^{2}d^{3}r_{1}d^{3}r_{2}dq$ (S22)

Evaluating for $\alpha\approx 20$ angstroms (a typical value for our defect
lattices), gives us $U\approx 1.1-1.2$ eV, which is in good agreement with our
ab-initio results of $\approx 1$ eV.

Figure S2: Accuracy of analytic spin splitting model. Analytic spin splitting
model plotted against the result obtained from DFT through PBE exchange-
correlation. Proximity to line at 45 degrees (black line) indicates accuracy
of the model. Figure S3: Accuracy of first order perturbation theory in
describing spin splitting. Each point defines a unique DFT calculation. x axis
defines the spin splitting as obtained through DFT, whereas the y axis is the
result obtained through first order perturbation theory with the densities and
wavefunctions obtained through DFT. The accuracy of the model is determined by
how close a given point is to the solid line at 45 degrees. Figure S4: Landau
undamped plasmons. Plasmonic dispersions (blue lines) sandwiched between the
lower and upper bounds of the Landau damping region. In all cases, the lower
(upper) bound corresponds to intraband (interband) electron-hole transitions.
Figure S5: Effective dielectric contrast. Plasmonic dispersions derived purely
from intraband transitions (dashed lines) and from both intraband and
interband transitions (solid lines).

##  Appendix G Buckling in Boron Substitution Lattices

In Fig. S6 we show the phonon dispersions of $C_{B}^{n\times n}$ for the case
where lattice optimization is constrained to the x-y plane (top row) and for
the case where the atoms are allowed to move out of plane (bottom row). The
imaginary phonon frequencies in the top row subplots (shown as negative
frequencies for convenience), indicate that all $C_{B}^{n\times n}$ lattices
are buckled. We show in Fig. S7 that the buckling leads to much lower decay
times.

Figure S6: Phonon dispersions and stability. (a) Phonons in planar systems of
$C_{B}^{n\times n}$. The imaginary phonon modes shown as dashed lines indicate
that a planar system is not stable. (b) Phonons of buckled $C_{B}^{n\times n}$
lattices. The lack of imaginary phonon frequencies indicates structural
stability. Note that $C_{B}^{\sqrt{3}\times\sqrt{3}}$ has small imaginary
phonon frequencies near $\Gamma$ on the order of 1meV. Since these are at
incommensurate points in the Brillouin zone and their mangitudes get smaller
as the phonon supercell is enlarged, they do not indicate a true instability.
(c) Same as (a) but for $C_{N}^{n\times n}$ Lack of imaginary phonon
frequencies indicates that planar $C_{N}^{n\times n}$ lattices are dynamically
stable.

##  Appendix H Plasmon Dispersion in Layered Systems

We consider an infinite number of equally spaced planes of defect lattices
(with lattice constant in the layered direction being $a$). To find the
plasmon dispersion, it is sufficient to consider one interface, between
regions 1 and 2. We take the ansatz $E_{2}=e^{ika}E_{1}$, where $E_{1,2}$ are
the electric fields in regions 1 and 2, respectively. We take as the electric
field in region 1:

$e^{iqx-i\omega t}(e^{-iQz}(Q,0,q)+Ee^{iQz}(Q,0,-q))$ (S23)

Continuity of the transverse part of the electric field at the interface
requires:

$e^{ika}(1+E)=e^{-iQa}+Ee^{iQa}\rightarrow
E=\frac{e^{ika}-e^{-iQa}}{e^{iQa}-e^{ika}}$ (S24)

The H field in region 1 is given by:

$H=-\frac{1}{\omega\mu}e^{iqx-i\omega
t}\frac{\omega^{2}}{c^{2}}(e^{-iQz}+Ee^{iQz})e_{y}=-\omega\epsilon
e^{iqx-i\omega t}(e^{-iQz}+Ee^{iQz})e_{y}$ (S25)

Therefore, the change of the transverse part of the H field gives us:

$\epsilon(e^{ika}-e^{-iQa}+E(e^{ika}-e^{iQa}))=-\frac{iQ}{q^{2}}\Pi(q,\omega)(1+E)e^{ika},$
(S26)

where $\Pi(q,\omega)$ is the polarizability Jablan et al. (2009). Substituting
for E and setting $Q=iq$, we get:

$2\epsilon(e^{ika}-e^{qa})=\frac{1}{q}\Pi(q,\omega)\frac{e^{-qa}-e^{qa}}{e^{-qa}-e^{ika}}e^{ika}\rightarrow
1=\frac{1}{2\epsilon
q}\frac{e^{-qa}-e^{qa}}{e^{-qa}-e^{ika}-e^{-ika}+e^{qa}}\Pi(q,\omega)$ (S27)

Therefore, we arrive at the final equation:

$1+\frac{\sinh(qa)}{\cosh(qa)-\cos(ka)}\frac{1}{2\epsilon q}\Pi(q,\omega)=0$
(S28)

Figure S7: Influence of buckling on electron-phonon coupling: We calculate the
decay time with the electron-phonon matrix elements of the relaxed buckled and
unbuckled $C_{B}^{n\times n}$ systems. The generically higher decay times for
the unbuckled systems indicate the enhancement of the electron-phonon coupling
through buckling.

##  Appendix I Convergence with Respect to Phonon Supercell Size

In Fig. S8 we show the convergence with regards to supercell size for two of
the defect lattices we studied.

Figure S8: Convergence of decay times with respect to phonon supercells: We
show the difference in decay times for two $C_{N}$ substitution lattices. As
seen, the results do not change qualitatively with phonon supercell size.

## References

  * Huang et al. (2022) P. Huang, M. Grzeszczyk, K. Vaklinova, K. Watanabe, T. Taniguchi, K. Novoselov, and M. Koperski, Physical Review B 106, 014107 (2022).
  * Lončarić et al. (2018) I. Lončarić, Z. Rukelj, V. M. Silkin, and V. Despoja, npj 2D Materials and Applications 2, 1 (2018).
  * Gjerding et al. (2021) M. N. Gjerding, A. Taghizadeh, A. Rasmussen, S. Ali, F. Bertoldo, T. Deilmann, N. R. Knøsgaard, M. Kruse, A. H. Larsen, S. Manti, et al., 2D Materials 8, 044002 (2021).
  * Lewandowski and Levitov (2019) C. Lewandowski and L. Levitov, Proceedings of the National Academy of Sciences 116, 20869 (2019).
  * Sundararaman et al. (2020) R. Sundararaman, T. Christensen, Y. Ping, N. Rivera, J. D. Joannopoulos, M. Soljačić, and P. Narang, Physical Review Materials 4, 074011 (2020).
  * Novoselov et al. (2005) K. S. Novoselov, D. Jiang, F. Schedin, T. Booth, V. Khotkevich, S. Morozov, and A. K. Geim, Proceedings of the National Academy of Sciences 102, 10451 (2005).
  * Boriskina et al. (2017) S. V. Boriskina, T. A. Cooper, L. Zeng, G. Ni, J. K. Tong, Y. Tsurimaki, Y. Huang, L. Meroueh, G. Mahan, and G. Chen, Advances in Optics and Photonics 9, 775 (2017).
  * Dai et al. (2019) S. Dai, W. Fang, N. Rivera, Y. Stehle, B.-Y. Jiang, J. Shen, R. Y. Tay, C. J. Ciccarino, Q. Ma, D. Rodan-Legrain, et al., Advanced materials 31, 1806603 (2019).
  * Rivera et al. (2019) N. Rivera, T. Christensen, and P. Narang, Nano Letters 19, 2653 (2019).
  * Novko et al. (2021) D. Novko, K. Lyon, D. J. Mowbray, and V. Despoja, Physical Review B 104, 115421 (2021).
  * Khurgin (2015a) J. B. Khurgin, Faraday discussions 178, 109 (2015a).
  * Rivera et al. (2016) N. Rivera, I. Kaminer, B. Zhen, J. D. Joannopoulos, and M. Soljačić, Science 353, 263 (2016).
  * Rivera and Kaminer (2020) N. Rivera and I. Kaminer, Nature Reviews Physics 2, 538 (2020).
  * Atwater and Polman (2010) H. A. Atwater and A. Polman, Nature materials 9, 205 (2010).
  * Langer et al. (2019) J. Langer, D. Jimenez de Aberasturi, J. Aizpurua, R. A. Alvarez-Puebla, B. Auguié, J. J. Baumberg, G. C. Bazan, S. E. Bell, A. Boisen, A. G. Brolo, et al., ACS nano 14, 28 (2019).
  * Homola (2003) J. Homola, Analytical and bioanalytical chemistry 377, 528 (2003).
  * Noginov et al. (2009) M. Noginov, G. Zhu, A. Belgrave, R. Bakker, V. Shalaev, E. Narimanov, S. Stout, E. Herz, T. Suteewong, and U. Wiesner, Nature 460, 1110 (2009).
  * Khurgin (2015b) J. B. Khurgin, Nature nanotechnology 10, 2 (2015b).
  * Gjerding et al. (2017) M. N. Gjerding, M. Pandey, and K. S. Thygesen, Nature communications 8, 1 (2017).
  * Cassabois et al. (2016) G. Cassabois, P. Valvin, and B. Gil, Nature photonics 10, 262 (2016).
  * Weston et al. (2018) L. Weston, D. Wickramaratne, M. Mackoit, A. Alkauskas, and C. Van de Walle, Physical Review B 97, 214104 (2018).
  * Liu et al. (2022) W. Liu, N.-J. Guo, S. Yu, Y. Meng, Z. Li, Y.-Z. Yang, Z.-A. Wang, X.-D. Zeng, L.-K. Xie, J.-F. Wang, et al., Materials for Quantum Technology (2022).
  * Giuliani and Quinn (1983) G. F. Giuliani and J. Quinn, Physical review letters 51, 919 (1983).
  * Lee and Mahanti (2012) M.-S. Lee and S. D. Mahanti, Physical Review B 85, 165149 (2012).
  * Mahan (2013) G. D. Mahan, _Many-particle physics_ (Springer Science & Business Media, 2013).
  * Agarwal et al. (2014) A. Agarwal, M. Polini, G. Vignale, and M. E. Flatté, Physical Review B 90, 155409 (2014).
  * da Jornada et al. (2020) F. H. da Jornada, L. Xian, A. Rubio, and S. G. Louie, Nature communications 11, 1013 (2020).
  * Stauber and Peres (2008) T. Stauber and N. Peres, Journal of Physics: Condensed Matter 20, 055002 (2008).
  * Allen (2015) P. B. Allen, Physical Review B 92, 054305 (2015).
  * Allen (1971) P. Allen, Physical Review B 3, 305 (1971).
  * Brown et al. (2016) A. M. Brown, R. Sundararaman, P. Narang, W. A. Goddard III, and H. A. Atwater, ACS nano 10, 957 (2016).
  * Pines and Schrieffer (1962) D. Pines and J. R. Schrieffer, Physical Review 125, 804 (1962).
  * Perdew et al. (1996) J. P. Perdew, K. Burke, and M. Ernzerhof, Physical review letters 77, 3865 (1996).
  * Dirac (1930) P. A. Dirac, in _Mathematical proceedings of the Cambridge philosophical society_ (Cambridge University Press, 1930), vol. 26, pp. 376–385.
  * Perdew et al. (1981) J. Perdew, E. McMullen, and A. Zunger, Physical Review A 23, 2785 (1981).
  * Profeta et al. (2012) G. Profeta, M. Calandra, and F. Mauri, Nature physics 8, 131 (2012).
  * Peres et al. (2008) N. Peres, T. Stauber, and A. C. Neto, EPL (Europhysics Letters) 84, 38002 (2008).
  * Gangadharaiah et al. (2008) S. Gangadharaiah, A. Farid, and E. Mishchenko, Physical review letters 100, 166802 (2008).
  * Polini et al. (2008) M. Polini, R. Asgari, G. Borghi, Y. Barlas, T. Pereg-Barnea, and A. MacDonald, Physical Review B 77, 081411 (2008).
  * Henriques et al. (2022) J. Henriques, B. Amorim, R. Ribeiro, and N. Peres, Physical Review B 105, 115421 (2022).
  * Haastrup et al. (2018) S. Haastrup, M. Strange, M. Pandey, T. Deilmann, P. S. Schmidt, N. F. Hinsche, M. N. Gjerding, D. Torelli, P. M. Larsen, A. C. Riis-Jensen, et al., 2D Materials 5, 042002 (2018).
  * Sundararaman et al. (2017) R. Sundararaman, K. Letchworth-Weaver, K. A. Schwarz, D. Gunceler, Y. Ozhabes, and T. Arias, SoftwareX 6, 278 (2017).
  * Schlipf and Gygi (2015) M. Schlipf and F. Gygi, Computer Physics Communications 196, 36 (2015).
  * Sundararaman and Arias (2013) R. Sundararaman and T. Arias, Physical Review B 87, 165122 (2013).
  * Souza et al. (2001) I. Souza, N. Marzari, and D. Vanderbilt, Physical Review B 65, 035109 (2001).
  * Kresse and Furthmüller (1996) G. Kresse and J. Furthmüller, Physical review B 54, 11169 (1996).
  * Freysoldt et al. (2009) C. Freysoldt, S. Boeck, and J. Neugebauer, Physical Review B 79, 241103 (2009).
  * Kumar et al. (2022) S. Kumar, C. Multunas, and R. Sundararaman, Physical Review Materials 6, 125201 (2022).
  * Wunsch et al. (2006) B. Wunsch, T. Stauber, F. Sols, and F. Guinea, New Journal of Physics 8, 318 (2006).
  * Giustino (2014) F. Giustino, _Materials modelling using density functional theory: properties and predictions_ (Oxford University Press, 2014).
  * Jablan et al. (2009) M. Jablan, H. Buljan, and M. Soljačić, Physical review B 80, 245435 (2009).
  * Adler (1962) S. L. Adler, Physical Review 126, 413 (1962).
  * Wiser (1963) N. Wiser, Physical Review 129, 62 (1963).


# Enumerating coprime permutations

Ashwin Sah  and  Mehtaab Sawhney  Department of Mathematics, Massachusetts
Institute of Technology, Cambridge, MA 02139, USA<EMAIL_ADDRESS>
###### Abstract.

Define a permutation $\sigma$ to be coprime if $\gcd(m,\sigma(m))=1$ for
$m\in[n]$. In this note, proving a recent conjecture of Pomerance, we prove
that the number of coprime permutations on $[n]$ is $n!\cdot(c+o(1))^{n}$
where

$c=\prod_{p\text{ prime }}\frac{(p-1)^{2(1-1/p)}}{p\cdot(p-2)^{(1-2/p)}}.$

The techniques involve entropy maximization for the upper bound, and a mixture
of number-theoretic bounds, permanent estimates, and the absorbing method for
the lower bound.

Sah and Sawhney were supported by NSF Graduate Research Fellowship Program
DGE-1745302. Sah was supported by the PD Soros Fellowship.

##  1\. Introduction

######  Definition 1.1.

Define a permutation $\sigma\colon[n]\to[n]$ to be _coprime_ if
$\operatorname{gcd}(m,\pi(m))=1$ for $m\in[n]$. Let $C(n)$ be the number of
coprime permutations on $[n]$.

In recent work, Pomerance [10] provided nontrivial bounds on $C(n)$ proving
that $n!/3.73^{n}<C(n)<n!/2.5^{n}$ for $n$ sufficiently large. Pomerance then
conjectured the existence of a constant $c$ such that $C(n)=n!(c+o(1))^{n}$
(see [10, Section 5]); our main objective is to provide such an asymptotic.

######  Theorem 1.2.

$C(n)=n!\bigg{(}\prod_{p\emph{ prime
}}\frac{(p-1)^{2(1-1/p)}}{p\cdot(p-2)^{(1-2/p)}}+\exp(-\Omega(\sqrt{\log
n\log\log n}))\bigg{)}^{n}.$

######  Remark.

The $p=2$ term is assumed to be $1/2$. The same method demonstrates that for
fixed $k$, $C_{k}(n)$, which counts permutations $\sigma\colon[n]\to[n]$ with
$\gcd(\ell,\sigma(\ell),k!)=1$ (see [10, Section 5]), satisfies

$C_{k}(n)=n!\bigg{(}\prod_{p\leq
k}\frac{(p-1)^{2(1-1/p)}}{p\cdot(p-2)^{(1-2/p)}}\bigg{)}^{n}\cdot\exp(O(\log
n)),$

though we forgo the details. McNew [8] has independently achieved similar
results for fixed $k$.

Theorem 1.2 is certainly not the only result regarding finding mappings
between sets of integers such that the labels of the matched edges are
coprime. One such example is work of Pomerance and Selfridge [12], answering a
question of Newman, demonstrating that such a coprime matching exists between
any set of $n$ consecutive integers $A$ and $[n]$. Another is work of Bohman
and Peng [2], in connection with the lonely runner conjecture, proving that
any two intervals $A,B$ of length $2m$ in $[n]$ with $m=2^{\Omega((\log\log
n)^{2})}$ have a coprime matching. This was strengthened in later work of
Pomerance [11] to require only $m=\Omega((\log n)^{2})$. We also mention that
enumerating $C(n)$ was also considered by Jackson [6] and is listed as A005326
on the OEIS [14].

###  1.1. Proof techniques

For the upper bound we study $C(n)$ by looking at permutations where
$\gcd(j,\pi(j))$ has no prime divisor at most a growing function $W$. Defining
$\tau=\prod_{p\leq W}p$ and examining the portion of integers $j$ which are
$a\bmod\tau$ and map to $b\bmod\tau$, we find that the problem essentially
immediately reduces to an _entropy maximization_ program. Using the concavity
and subadditivity of entropy as well as the Chinese Remainder Theorem, we can
obtain the upper bound, which is done in Section 2. Note that as written, this
approach essentially necessarily requires $W=O(\log n)$ as otherwise $\tau>n$
and the integers $\\{1,\ldots,n\\}$ are not equidistributed in residue classes
$\bmod\tau$. In order to do better, we (a) restrict attention only to whether
any prime $p\leq W$ divides or does not divide an integer $j\in[n]$, and (b)
focus only on possible prime divisor patterns which produce large buckets
(specifically, buckets which have at most $k\approx(\log n)^{1/2}$ prime
divisors), and treat the rest as an error term. Since $\prod_{p\leq
W}(1-p^{-1})=\Theta(1/\log W)$ is a negligible factor, the condition on $W$ we
obtain is roughly that $W^{k}\leq n^{1/100}$, which allows for consideration
of significantly more prime divisors.

For the proof of the lower bound in Section 3, the role of larger prime
factors immediately becomes crucial. Let $B_{S}$ be the set of integers in
$[n]$ with $S$ being the precise set of prime divisors at most $W$. For the
sake of this discussion $W=\exp(c(\log n)^{1/2}(\log\log n)^{1/2})$ and
$k=(\log n)^{1/2}(\log\log n)^{-1/2}$. Our proof is based on first taking a
“template”, assigning the appropriate “entropy-maximizing” fraction of $B_{S}$
to map to $B_{S^{\prime}}$. Noting the size of $W$ it follows immediately that
(say)

$\sum_{\begin{subarray}{c}p|n\\\ p>W\end{subarray}}1/p\leq W^{-1/2}$

as $n$ has at most $2\log n/(\log\log n)$ distinct prime divisors. This
suggests that the effect of “large primes” is minimal and therefore the upper
bound proof should be reversible.

To make this argument rigorous, one chooses random subsets of $B_{S}$ to map
to $B_{S^{\prime}}$ in accordance with the entropy maximizing distribution.
One then proves that with high probability the associated coprimality graph
corresponding to $B_{S},B_{S^{\prime}}$, if $S$ and $S^{\prime}$ have size at
most $k$, has at most approximately a $W^{-1/2}$ fraction of edges missing at
any vertex. This is done via the Chernoff bound and the prior estimate
regarding large prime factors. Combined with estimates regarding the number of
perfect matchings in graph of high minimum degree due to Alon, Rödl, and
Ruciński [1], one can (essentially) handle matching all integers with at most
$k$ prime factors at most $W$.

For the remaining integers, we “absorb” this small subset into the set of
integers divisible by no primes at most $W$, which is a very large set, namely
of size $\Omega(n/\log W)$. On the other hand, the number of integers with at
least $k$ prime factors less than $W$ is $n\exp(-\Omega(k\log k))$. See [15]
for a more extensive discussion of the absorbing method and its history. Our
approach here can be thought of as analogous to a “robust absorption
strategy”, which has appeared for instance in work of Montgomery [9], in which
bipartite graphs with a “resilient matching” property are utilized. Finally,
we note that we have omitted certain special considerations at the prime
$p=2$, which are handled via a reduction of Pomerance [10] (Lemma 3.1).

###  1.2. Further directions

We end by mentioning that we believe obtaining sharper asymptotics for $C(n)$
is an interesting open problem; our methods are fundamentally limited to error
terms that are of the form $n^{-o(1)}$ as we crucially rely on considering
buckets whose shared prime divisors multiply to less than $n^{1/2}$ (say) in
order to guarantee that the interval has the expected number of elements in
such buckets. Thus either we cannot consider many primes, or cannot consider
buckets with many prime divisors, either of which is a serious issue. Given
this, error terms of the form $O(n^{-\Omega(1)})$ would require a substantial
modification of our methods and even the truth of such a statement is not
obvious. Additionally, obtaining an asymptotic expansion for $C(n)$ appears
unlikely (as speculated by Jackson [6]).

We note that a reduction of Pomerance [10] shows it is enough to count coprime
bijections $f\colon\\{1,3,\ldots,2n-1\\}\to[n]$. This is the number of
matchings in a bipartite graph $i\in\\{1,\ldots,2n-1\\}$ on the left and
$j\in[n]$ on the right are connected if $\gcd(i,j)=1$. Our proof can be
interpreted as considering the non-edges that occur due to small primes and
noting that the “remaining non-edges” are sparse. For sharper bounds, it may
be possible to analyze an edge-removal process that targets the maximum
entropy distribution described in Section 1.1: at each step, uniformly
randomly choose an edge $(i,j)$ in the current graph with probability say
proportional to some $g(i,j)$ and then delete $i$ from the left and $j$ from
the right. Here we might take

$g(i,j)=\prod_{2<p\leq n}g_{p}(\mathbbm{1}_{p|i},\mathbbm{1}_{p|j})$

with $g_{p}(0,0)=(p-2)/p$, $g_{p}(0,1)=g_{p}(1,0)=1/p$, and $g_{p}(1,1)=0$.
Depending on how long one can control this process, the remaining vertices can
be handled by the absorbing method. Our proof can be thought of as very
similar to analyzing such a process where we truncate the weights at $p\leq
W$.

Techniques involving “weighted removal processes” or that target an “entropy-
maximizing distribution” along with the absorbing method have been recently
used to give sharp estimates on the $n$-queens problem [13], and may be
applicable to further situations where one wishes to count objects with
structured biases (such as number-theoretic considerations, as in the present
work). This in turn builds on recent powerful techniques introduced for
counting in quasirandom situations such as counting combinatorial designs [7],
in which polynomially strong error terms are known. It remains interesting to
see the limits of such techniques in number-theoretic settings.

### Notation

Throughout this paper we will reserve $p$ for variables which range over the
primes. We will also implicitly assume $n$ is sufficiently large for certain
inequalities to hold. We write $a=b\pm c$ to mean that $a\in[b-c,b+c]$, and
$f=O(g)$ means $|f|\leq Cg$ for some absolute constant $C$ while $f=\Omega(g)$
means $f\geq C|g|$ for some absolute $C$. We use $\log^{(k)}n$ to denote the
$k$-times iterated logarithm, e.g. $\log^{(5)}n=\log\log\log\log\log n$. We
may sometimes omit floors and ceilings without comment if not essential to the
situation. We say an event occurs whp (with high probability) if it occurs
with probability tending to $1$ as $n\to\infty$. Given a set $S$, we write
$\binom{S}{k},\binom{S}{\leq k},\binom{S}{<k}$ for the collection of subsets
of $S$ of size $k$, size at most $k$, and size less than $k$, respectively.
Finally, as above, we often refer to a bijective function or permutation as a
matching between the domain and codomain, especially when dealing with extra
constraints encoded via a graph.

### Acknowledgements

We thank Carl Pomerance for comments motivating our proof of Lemma 2.1 and for
communicating the results of McNew [8] regarding $C_{k}(n)$.

##  2\. Upper Bound

###  2.1. Setup

Fix $W=\exp(2^{-10}(\log n\log\log n)^{1/2})$ and $k=2^{-5}(\log n/\log\log
n)^{1/2}$. Let $\alpha=\exp(-k\log k)$, which will roughly correspond to the
error term we derive. Our upper bound on the number of coprime permutations
$\sigma\colon[n]\to[n]$ comes via enumerating permutations such that
$p\nmid\gcd(j,\sigma(j))$ for all $p\leq W$. Let $\mathcal{P}_{W}=\\{p\colon
p\leq W\\}$ be the primes which are at most $W$.

We first define a series of sets with respect to the divisibility structure of
$\mathcal{P}_{W}$. For each $S\subseteq\mathcal{P}_{W}$, let

$A_{S}=[n]\cap\bigcap_{p\in S}\\{m\colon
p|m\\}\cap\bigcap_{p\in\mathcal{P}_{W}\setminus S}\\{m\colon p\nmid m\\}$

and

$A_{\geq
k}=[n]\setminus\bigg{(}\bigcup_{S\in\binom{\mathcal{P}_{W}}{<k}}A_{S}\bigg{)}.$

To start, we prove that $A_{S}$ has the “predicted size” from an independent
heuristic on small primes if $|S|<k$, while $A_{\geq k}$ is a vanishing small
proportion of the set. This will allow us to efficiently encode the
divisibility structure of all primes at most $W$. Note that this is nontrivial
as a more naive bound based on the Chinese Remainder Theorem (CRT) would
restrict us to having $W$ such that $\prod_{p\leq W}p\leq n^{1/2}$ (say) and
hence $W\leq\log n$.

######  Lemma 2.1.

For all $S\subseteq\mathcal{P}_{W}$ with $|S|<k$ we have

$|A_{S}|=(1\pm\alpha^{2})n\prod_{p\in
S}p^{-1}\prod_{p\in\mathcal{P}_{W}\setminus S}(1-p^{-1}).$

Furthermore,

$|A_{\geq k}|\leq n\alpha^{1/2}.$

###### Proof.

We first handle $A_{\geq k}$. Let $Y$ be a uniformly random positive integer
in $[n]$ and $Y_{p}$ be the indicator of $p|Y$ for $1\leq p\leq W$. Note that

$\displaystyle\mathbb{P}[Y\in A_{\geq k}]$ $\displaystyle\leq\sum_{1\leq
p_{1}<\cdots<p_{k}\leq W}\mathbb{P}\bigg{[}\bigcap_{1\leq j\leq
k}Y_{p_{j}}\bigg{]}=\sum_{1\leq p_{1}<\cdots<p_{k}\leq
W}\mathbb{P}\bigg{[}\bigwedge_{1\leq j\leq k}(p_{j}|Y)\bigg{]}$
$\displaystyle\leq 2\sum_{1\leq p_{1}<\cdots<p_{k}\leq W}\prod_{1\leq j\leq
k}\frac{1}{p_{j}}\leq\frac{2}{k!}\bigg{(}\sum_{1\leq p\leq
W}\frac{1}{p}\bigg{)}^{k}\leq\exp(-k\log k/2)$ (2.1)

where we have used that $W^{k}\leq n$ to derive the probability that $Y$ is
divisible by $\prod_{1\leq j\leq k}p_{j}$ and used that $\sum_{1\leq p\leq
W}1/p=O(\log\log n)$. To derive the necessary bound on $A_{S}$, we proceed via
the Principle of Inclusion-Exclusion (PIE), in particular with the use of
Bonferroni’s inequalities [3]. In particular, we have

$\displaystyle\mathbb{P}[Y\in A_{S}]$
$\displaystyle=\mathbb{E}\bigg{[}\prod_{p\in
S}Y_{p}\prod_{p\in\mathcal{P}_{W}\setminus S}(1-Y_{p})\bigg{]}$
$\displaystyle=\sum_{0\leq j\leq
4k}(-1)^{j}\sum_{T\subseteq\binom{\mathcal{P}_{W}\setminus
S}{j}}\mathbb{E}\bigg{[}\prod_{p\in S}Y_{p}\prod_{p\in
T}Y_{p}\bigg{]}\pm\sum_{T\subseteq\binom{\mathcal{P}_{W}\setminus
S}{4k}}\mathbb{E}\bigg{[}\prod_{p\in S}Y_{p}\prod_{p\in T}Y_{p}\bigg{]}$
$\displaystyle=\sum_{0\leq j\leq
4k}(-1)^{j}\sum_{T\subseteq\binom{\mathcal{P}_{W}\setminus S}{j}}\prod_{p\in
S\cup T}p^{-1}\pm\bigg{(}\bigg{(}\prod_{p\in S}p^{-1}\bigg{)}\exp(-3k\log
k)+(2W)^{4k}/n\bigg{)}$ $\displaystyle=\prod_{p\in
S}p^{-1}\prod_{p\in\mathcal{P}_{W}\setminus
S}(1-p^{-1})\pm\bigg{(}\bigg{(}\prod_{p\in S}p^{-1}\bigg{)}\exp(-3k\log
k)+(2W)^{4k}/n\bigg{)}$ $\displaystyle=(1\pm\exp(-2k\log k))\prod_{p\in
S}p^{-1}\prod_{p\in\mathcal{P}_{W}\setminus S}(1-p^{-1}).$

The first error term in the third line comes from a bound similar to 2.1,
while the second error term in the third line comes from the slight deviations
in probability that a $Y$ is divisible by some product of primes at most
$W^{k}\leq n$. The fourth line comes from bounding the tail terms of the
product expansion of $\prod_{p\in S}p^{-1}\prod_{p\in\mathcal{P}_{W}\setminus
S}(1-p^{-1})$ (e.g. by use of the Bonferroni inequalities once again) and the
fifth line follows upon converting additive error to multiplicative error,
using $W^{k}\leq n^{1/8}$ and
$\prod_{p\in\mathcal{P}_{W}}(1-p^{-1})=\Omega(1/\log W)$. This completes the
proof. ∎

Given a coprime permutation $\sigma\colon[n]\to[n]$ and
$S_{1},S_{2}\subseteq\mathcal{P}_{W}$ we let

$\beta_{S_{1},S_{2}}(\sigma)=\\#\\{m\in[n]\colon m\in
A_{S_{1}},~{}\sigma(m)\in A_{S_{2}}\\}$

and $\beta(\sigma)$ be the tuple of all such counts. Note there are at most
$(2W)^{k}\leq\exp(k\log k)\leq n^{1/4}$ such sets $S$ with $|S|<k$ and that by
Lemma 2.1, we have $|A_{\geq k}|\leq n\alpha^{1/2}$. This implies that there
are at most

$(n+1)^{\sqrt{n}}\cdot(2n^{2})^{n\alpha^{1/2}}\leq e^{n\alpha^{1/3}}$

possible choices for $\beta(\sigma)$: first, for every pair $S_{1},S_{2}$ with
$|S_{1}|,|S_{2}|<k$ we choose up to $n+1$ possibilities, and then we bound the
ways to simultaneously choose all $\beta_{S,S^{\prime}}(\sigma)$ for
$S\subseteq\mathcal{P}_{W}$ and $|S^{\prime}|\geq k$ (and similar for
$\beta_{S^{\prime},S}(\sigma)$), which is a total of $t\leq n^{2}$ values
which must sum to $|A_{\geq k}|$.

Therefore we have that

$C(n)\leq e^{n\alpha^{1/3}}\cdot\sup_{\beta}\\{\sigma\colon[n]\to[n]\text{
coprime}\colon\beta(\sigma)=\beta\\}$ (2.2)

Finally we define the following parameters:

$\rho_{S_{1},S_{2}}=\beta_{S_{1},S_{2}}/n,~{}\beta_{S_{1},\cdot}=\sum_{S_{2}}\beta_{S_{1},S_{2}}=|A_{S_{1}}|,~{}\beta_{\cdot,S_{2}}=\sum_{S_{1}}\beta_{S_{1},S_{2}}=|A_{S_{2}}|.$

###  2.2. Completing the upper bound

We now bound the number of choices of $\sigma$ given $\beta(\sigma)$. This
will complete the proof of the upper bound due to 2.2. We bound the number of
permutations $\sigma$ satisfying $\beta(\sigma)=\beta$ as follows:

  * •

For each pair of subsets $S_{1},S_{2}\subseteq\mathcal{P}_{W}$, choose which
$\beta_{S_{1},S_{2}}$ elements of $A_{S_{1}}$ which have image in $A_{S_{2}}$.
Call such elements $V_{S_{1}\to S_{2}}$.

  * •

For each pair of subsets $S_{1},S_{2}\subseteq\mathcal{P}_{W}$, choose which
$\beta_{S_{1},S_{2}}$ elements of $A_{S_{2}}$ which have preimage in
$A_{S_{1}}$. Call such elements $V_{S_{1}\leftarrow S_{2}}$.

  * •

Determine a precise mapping between $V_{S_{1}\to S_{2}}$ and
$V_{S_{1}\leftarrow S_{2}}$.

Note also that if $\rho_{S_{1},S_{2}}\neq 0$ where $S_{1}\cap
S_{2}\neq\emptyset$, then the expression within the supremum in 2.2 is
actually $0$. Thus we can restrict attention to values of $\beta$ where
$\rho_{S_{1},S_{2}}=0$ when $S_{1}\cap S_{2}\neq\emptyset$.

Given this procedure, we see that the number of $\sigma$ such that
$\beta(\sigma)=\beta$ is bounded by

$\displaystyle\prod_{S_{1},S_{2}\subseteq\mathcal{P}_{W}}\beta_{S_{1},S_{2}}!\prod_{S_{1}\subseteq\mathcal{P}_{W}}\frac{\beta_{S_{1},\cdot}!}{\prod_{S_{2}\subseteq\mathcal{P}_{W}}\beta_{S_{1},S_{2}}!}\prod_{S_{2}\subseteq\mathcal{P}_{W}}\frac{\beta_{\cdot,S_{2}}!}{\prod_{S_{1}\subseteq\mathcal{P}_{W}}\beta_{S_{1},S_{2}}!}$
$\displaystyle=\prod_{S_{1},S_{2}\subseteq\mathcal{P}_{W}}(\beta_{S_{1},S_{2}}!)^{-1}\prod_{S_{1}\subseteq\mathcal{P}_{W}}\beta_{S_{1},\cdot}!\prod_{S_{2}\subseteq\mathcal{P}_{W}}\beta_{\cdot,S_{2}}!$
$\displaystyle\leq|A_{\geq
k}|!^{2}\prod_{|S_{1}|,|S_{2}|<k}(\beta_{S_{1},S_{2}}!)^{-1}\prod_{|S_{1}|<k}\beta_{S_{1},\cdot}!\prod_{|S_{2}|<k}\beta_{\cdot,S_{2}}!$
$\displaystyle\leq
e^{n\alpha^{1/3}}\cdot\prod_{|S_{1}|,|S_{2}|<k}(\beta_{S_{1},S_{2}}!)^{-1}\prod_{|S|<k}|A_{S}|!^{2}$
(2.3)

where we used that factorials are log-convex in the third line to combine the
terms for $|S|\geq k$, and used Lemma 2.1 and Stirling’s approximation in the
fourth.

We now bound each term in the above product. Write $n_{S}=n\prod_{p\in
S}p^{-1}\prod_{p\in\mathcal{P}_{W}\setminus S}(1-p^{-1})$. Noting that $n!\leq
3(n+1)(n/e)^{n}$ for all $n\geq 1$ and $(2W)^{k}\leq n^{1/4}$, and applying
Lemma 2.1 we find that

$\displaystyle\prod_{|S|\leq k}|A_{S}|!$ $\displaystyle\leq
e^{n\alpha^{1/3}}\cdot\prod_{|S|<k}((1+\alpha^{2})n_{S}/e)^{(1+\alpha^{2})n_{S}}$
$\displaystyle\leq
e^{2n\alpha^{1/3}}\cdot\prod_{S\subseteq\mathcal{P}_{W}}\bigg{(}(n/e)\prod_{p\in
S}p^{-1}\prod_{p\in\mathcal{P}_{W}\setminus S}(1-p^{-1})\bigg{)}^{n\prod_{p\in
S}p^{-1}\prod_{p\in\mathcal{P}_{W}\setminus S}(1-p^{-1})}$
$\displaystyle=e^{2n\alpha^{1/3}}\cdot(n/e)^{n}\prod_{p\in\mathcal{P}_{W}}(1-1/p)^{n(1-1/p)}(1/p)^{n/p}$

Next noting that $n!\geq(n/e)^{n}$ for all $n\geq 1$, we find that

$\displaystyle\prod_{|S_{1}|,|S_{2}|<k}(\beta_{S_{1},S_{2}}!)^{-1}$
$\displaystyle\leq
e^{n}\prod_{|S_{1}|,|S_{2}|<k}{\beta_{S_{1},S_{2}}}^{-\beta_{S_{1},S_{2}}}\leq
e^{n\alpha^{1/3}}\cdot
e^{n}\prod_{S_{1},S_{2}\subseteq\mathcal{P}_{W}}{\beta_{S_{1},S_{2}}}^{-\beta_{S_{1},S_{2}}}$
$\displaystyle=e^{n\alpha^{1/3}}\cdot
e^{n}\exp\Big{(}-n\sum_{S_{1},S_{2}}\rho_{S_{1},S_{2}}\log\beta_{S_{1},S_{2}}\Big{)}$
$\displaystyle=e^{n\alpha^{1/3}}\cdot(e/n)^{n}\exp\Big{(}-n\sum_{S_{1},S_{2}}\rho_{S_{1},S_{2}}\log\rho_{S_{1},S_{2}}\Big{)}.$

Let $Z$ be the random variable which takes on the value $(S_{1},S_{2})$ with
probability $\rho_{S_{1},S_{2}}$. Treating each set $S_{i}$ as a vector in
$\\{0,1\\}^{\mathcal{P}_{W}}$, we may equivalently view $Z=((\mathbbm{1}_{p\in
S_{1}},\mathbbm{1}_{p\in
S_{2}}))_{p\in\mathcal{P}_{W}}\in\\{0,1\\}^{2|\mathcal{P}_{W}|}$. Let $Z_{p}$
be $(\mathbbm{1}_{p\in S_{1}},\mathbbm{1}_{p\in S_{2}})$ when $(S_{1},S_{2})$
is drawn according to $Z$, and note that $Z$ is determined upon knowing
$Z_{p}$ for all $p\in\mathcal{P}_{W}$. By the subadditivity of entropy we
therefore find that

$\exp\Big{(}-n\sum_{S_{1},S_{2}}\rho_{S_{1},S_{2}}\log\rho_{S_{1},S_{2}}\Big{)}=\exp(nH(Z))\leq\prod_{p\in\mathcal{P}_{W}}\exp(nH(Z_{p})).$

Note that $Z_{p}$ is supported on $(0,0),(1,0),(0,1)$ due to the constraint
$\rho_{S_{1},S_{2}}=0$ whenever we have $p\in S_{1}\cap S_{2}$. Additionally,
$\mathbb{P}[Z_{p}=(1,\cdot)]=\mathbb{P}[Z_{p}=(\cdot,1)]=\lfloor n/p\rfloor/n$
since $\sigma$ is a permutation and $|\\{p|\ell\\}\cap[n]|=\lfloor
n/p\rfloor$. These constraints determine the law of $Z_{p}$ and therefore
$H(Z_{p})$. In particular, it follows that
$\mathbb{P}[Z_{p}=(1,0)]=\mathbb{P}[Z_{p}=(0,1)]=1/p\pm 1/n$ and
$\mathbb{P}[Z_{p}=(0,0)]=1-2/p\pm 2/n$.

Putting everything together into 2.3, this gives an upper bound of

$\displaystyle
e^{7n\alpha^{1/3}}\cdot(n/e)^{n}\bigg{(}\prod_{p\in\mathcal{P}_{W}}(1-1/p)^{n(1-1/p)}(1/p)^{n/p}\bigg{)}^{2}\prod_{p\in\mathcal{P}_{W}}(1-2/p)^{-n(1-2/p)}(1/p)^{-2n/p}$
$\displaystyle\qquad=e^{7n\alpha^{1/3}}\cdot(n/e)^{n}\prod_{p\in\mathcal{P}_{W}}(p-1)^{2(1-1/p)n}(p-2)^{-n(1-2/p)}p^{-n}$
$\displaystyle\qquad\leq
e^{7n\alpha^{1/3}+n/W}\cdot(n/e)^{n}\prod_{p}(p-1)^{2(1-1/p)n}(p-2)^{-n(1-2/p)}p^{-n}$

where we have used that $(p-1)^{2(1-1/p)}(p-2)^{-(1-2/p)}p^{-1}=1+O(1/p^{2})$
which implies the inequality
$\prod_{p>W}(p-1)^{2(1-1/p)}(p-2)^{-(1-2/p)}p^{-1}=1+O(1/(W\log
W))\leq\exp(1/W)$. The desired result follows immediately using our
expressions for $\alpha,W$.

##  3\. Lower Bound

###  3.1. Setup and preliminaries

Recall that we have set $W=\exp(2^{-10}(\log n\log\log n)^{1/2})$ and
$k=2^{-5}(\log n/\log\log n)^{1/2}$, as well as $\alpha=\exp(-k\log k)$.

For the lower bound, note that any coprime permutation must map even numbers
to odd numbers, which are sets of approximately the same size. In order to
witness this phenomenon, we use the following reduction of Pomerance [10].
Define $[n]_{\mathrm{o}}=\\{1,3,\ldots,2n-1\\}$, the first $n$ odd positive
integers.

######  Lemma 3.1 ([10, Lemma 1]).

Let $C_{0}(n)$ be the number of bijective functions
$f\colon[n]_{\mathrm{o}}\to[n]$ with $\gcd(j,f(j))=1$ for all
$j\in\\{1,3,\ldots,2n-1\\}$. Then $C(2n)=C_{0}(n)^{2}$ and $C(2n+1)\geq
2C_{0}(n-1)^{2}$.

We will require slightly modified sets in order to capture the divisibility
structures of $[n]$ and $[n]_{\mathrm{o}}$. Let
$\mathcal{P}_{W}^{\prime}=\mathcal{P}_{W}\setminus\\{2\\}$. For
$S\subseteq\mathcal{P}_{W}^{\prime}$ define

$B_{S}=[n]\cap\bigcap_{p\in S}\\{m\colon
p|m\\}\cap\bigcap_{p\in\mathcal{P}_{W}^{\prime}\setminus S}\\{m\colon p\nmid
m\\},~{}B_{\geq
k}=[n]\setminus\bigcup_{S\in\binom{\mathcal{P}_{W}^{\prime}}{<k}}B_{S}$

and the analogous

$C_{S}=[n]_{\mathrm{o}}\cap\bigcap_{p\in S}\\{m\colon
p|m\\}\cap\bigcap_{p\in\mathcal{P}_{W}^{\prime}\setminus S}\\{m\colon p\nmid
m\\},~{}C_{\geq
k}=[n]_{\mathrm{o}}\setminus\bigcup_{S\in\binom{\mathcal{P}_{W}^{\prime}}{<k}}C_{S}.$

Similar to Lemma 2.1, we have the following control on the set sizes $B_{S}$,
$C_{S}$.

######  Lemma 3.2.

For all $S\subseteq\mathcal{P}_{W}^{\prime}$ with $|S|<k$ we have,

$|B_{S}|,|C_{S}|=(1\pm\alpha^{2})n\prod_{p\in
S}p^{-1}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus S}(1-p^{-1}).$

Furthermore,

$|B_{\geq k}|,|C_{\geq k}|\leq\alpha^{1/2}n.$

For the remainder of this section we focus on providing an asymptotic lower
bound on $C_{0}(n)$, which we treat as a matching problem between
$[n]_{\mathrm{o}}$ and $[n]$ with forbidden edges corresponding to non-coprime
pairs. Now, heuristically, the lower bound should match the bounding scheme
presented in Section 2 (ignoring the prime $p=2$). However, to actually ensure
we avoid coprimality with respect to larger primes $p$, we will require a
number of estimates. The crucial point for our lower bound is noting that the
divisibility structure of primes larger than $W$ removes a sparse graph (of
density bounded by say $W^{-1/2}$).

######  Fact 3.3.

For $\ell\in[2n]$ we have that

$\sum_{\begin{subarray}{c}p|\ell\\\ p>W\end{subarray}}\frac{1}{p}<W^{-1/2}.$

###### Proof.

This is immediate noting that at most $2\log n/(\log\log n)$ distinct primes
divide an integer $\ell\in[2n]$. ∎

We will also require that for a given “large class” only a small portion of
numbers are divisible by a fixed $p$.

######  Lemma 3.4.

For $W<p^{\prime}\leq n^{1/3}$ and $|S|<k$, we have that

$|B_{S}\cap\\{m\colon p^{\prime}|m\\}|\leq
2|A_{S}|/p^{\prime},~{}|C_{S}\cap\\{m\colon p|m\\}|\leq 2|B_{S}|/p^{\prime}.$

For $p^{\prime}>n^{1/3}$ we have that

$|B_{S}\cap\\{m\colon p^{\prime}|m\\}|\leq 2n^{2/3},~{}|C_{S}\cap\\{m\colon
p|m\\}|\leq 2n^{2/3}.$

###### Proof.

The second claim follows immediately as $|[2n]\cap\\{m:p^{\prime}|m\\}|\leq
2n^{2/3}$. For the first claim, we focus on $B_{S}$ as the analogous claim for
$C_{S}$ follows similarly. Let $Y$ and $Y_{p}$ be as in the proof of Lemma
2.1. Note that

$\displaystyle\mathbb{P}[Y\in A_{S}$ $\displaystyle\wedge
Y_{p^{\prime}}]=\mathbb{E}\bigg{[}Y_{p^{\prime}}\prod_{p\in
S}Y_{p}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus S}(1-Y_{p})\bigg{]}$
$\displaystyle=\sum_{0\leq j\leq
4k}(-1)^{j}\sum_{T\subseteq\binom{\mathcal{P}_{W}^{\prime}\setminus
S}{j}}\mathbb{E}\bigg{[}Y_{p^{\prime}}\prod_{p\in S}Y_{p}\prod_{p\in
T}Y_{p}\bigg{]}\pm\sum_{T\subseteq\binom{\mathcal{P}_{W}^{\prime}\setminus
S}{4k}}\mathbb{E}\bigg{[}Y_{p^{\prime}}\prod_{p\in S}Y_{p}\prod_{p\in
T}Y_{p}\bigg{]}$ $\displaystyle=\sum_{0\leq j\leq
4k}(-1)^{j}\sum_{T\subseteq\binom{\mathcal{P}_{W}^{\prime}\setminus
S}{j}}(p^{\prime})^{-1}\prod_{p\in S\cup
T}p^{-1}\pm\bigg{(}(p^{\prime})^{-1}\bigg{(}\prod_{p\in
S}p^{-1}\bigg{)}\exp(-3k\log k)+(2W)^{4k}/n\bigg{)}$
$\displaystyle=(1\pm\alpha^{2})(p^{\prime})^{-1}\prod_{p\in
S}p^{-1}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus S}(1-p^{-1})$

where we have used that $n^{1/3}(W)^{4k}\leq n^{1/4}$. The result follows
immediately. ∎

We will also require that in any bipartite graph with balanced vertex parts
and a very large minimum degree there are a large number of perfect matchings.
The argument is identical to that of Alon, Rödl, and Ruciński [1, Theorem 1],
but we reprove the result in order to quantify the dependencies.

######  Lemma 3.5.

Let $G=(A\cup B,E)$ be a bipartite graph such that $|A|=|B|=n$ and
$\Delta=\Delta(K_{A,B}\setminus G)$, the maximum degree of the bipartite
complement, satisfies $\Delta\leq n/3$. Then there are at least
$((n-2\Delta)/e)^{n}$ perfect matchings in $G$.

###### Proof.

By [1, Theorem 3], if for any $X\subseteq A$, $Y\subseteq B$,

$k(|X|+|Y|)+E(A\setminus X,B\setminus Y)\geq kn$

holds, then $G$ has a $k$-factor (and therefore a subgraph of $G$ which is
$k$-regular). Let $k=n-2\Delta$ and note that for $|X|+|Y|\geq n$ the above
inequality is immediate. If $|X|+|Y|\leq n$, suppose that $|X|\geq|Y|$ (the
other case being analogous) and note that $E(A\setminus X,B\setminus
Y)\geq(n-|X|)(n-|Y|)-(n-|X|)\Delta=(n-|X|)(n-|Y|-\Delta)$ and hence

$k(|X|+|Y|)+E(A\setminus X,B\setminus Y)\geq
k(|X|+|Y|)+(n-|X|)(n-|Y|-\Delta).$

Therefore it suffices to check that $k(|X|+|Y|)+(n-|X|)(n-|Y|-\Delta)\geq kn$
which is equivalent to

$k(n-|X|-|Y|)\leq(n-|X|)(n-|Y|-\Delta).$

Fix $|X|+|Y|=t$ and imagine varying $|X|$ with $t\geq|X|\geq|Y|=t-|X|$. By
concavity of the right side, it suffices to check when $|X|=t,|Y|=0$ and
$|X|=|Y|=t/2$. In the first case the inequality follows as $k\leq n-\Delta$
and in the second case the inequality is equivalent to

$k(n-t)\leq(n-t/2)(n-t/2-\Delta).$

Viewing the difference of the right and left hand sides as functions of $t$,
the derivative is seen to be $t/2+k-n+\Delta/2$ and thus the inequality is
closest when $t=2n-2k-\Delta=3\Delta$. It becomes

$(n-2\Delta)(n-3\Delta)\leq(n-3\Delta/2)(n-5\Delta/2)$

which is trivially true as $n\geq 3\Delta$.

Now let $G^{\prime}$ be the corresponding subgraph of $G$ which is exactly
$k$-regular. Then by the solution of the van der Waerden conjecture on
permanents due to Egorychev [4] and Falikman [5] it follows that the number of
perfect matchings in $G^{\prime}$ is at least $(k/n)^{n}n!\geq(k/e)^{n}$. ∎

Finally we will repeatedly require the Chernoff bound.

######  Lemma 3.6 (Chernoff bound).

Let $X$ be either:

  * •

a sum of independent random variables, each of which takes values in $[0,1]$

  * •

hypergeometrically distributed (with any parameters).

Then for any $\delta\geq 0$ we have

$\mathbb{P}[X\geq(1+\delta)\mathbb{E}X]\leq\exp(-\delta^{2}\mathbb{E}X/(2+\delta)).$

###  3.2. Proof of the lower bound

In order to complete the proof of Theorem 1.2, we devise a procedure that
produces a triplet $(\sigma,R_{B},R_{C})$ where $R_{B}\subseteq[n]$ and
$R_{C}\subseteq[n]_{\mathrm{o}}$ with $|R_{B}|,|R_{C}|\leq n\exp(n\alpha^{c})$
and such that $\sigma$ is coprime. The number of possible outcomes of such a
procedure is bounded above by $C(n)\exp(n\alpha^{c/2})$ due to the size bounds
on $R_{B},R_{C}$ and therefore it suffices to provide a lower bound on the
number possible outcomes.

Let $\mathcal{M}_{B}$ denote the set of labels
$\bigcup_{\begin{subarray}{c}S\subseteq\mathcal{P}_{W}^{\prime}\\\
|S|<k\end{subarray}}\\{S\\}\cup\\{\ast\\}$ and let $\mathcal{M}_{C}$ denote
the set of labels
$\bigcup_{\begin{subarray}{c}S\subseteq\mathcal{P}_{W}^{\prime}\\\
|S|<k\end{subarray}}\\{S\\}\cup\\{\ast\\}$. We first construct (via a
randomized procedure) a pair of maps $\phi_{1}\colon[n]\to\mathcal{M}_{B}$ and
$\phi_{2}\colon[n]\to\mathcal{M}_{C}$ which will serve as a template for the
procedure producing $(\sigma,R_{B},R_{C})$.

######  Lemma 3.7.

There exist assignment functions $\phi_{1}\colon[n]\to\mathcal{M}_{B}$ and
$\phi_{2}\colon[n]\to\mathcal{M}_{C}$ such that:

  * •

$\phi_{1}(m_{1})=\ast$ for $m_{1}\in B_{\geq k}$ and $\phi_{2}(m_{2})=\ast$
for $m_{2}\in C_{\geq k}$.

  * •

If $S\cap T\neq\emptyset$ then $x\in B_{S}$ implies $\phi_{1}(x)\neq T$, and
$y\in C_{T}$ implies $\phi_{2}(y)\neq S$.

  * •

For $|S_{1}|,|S_{2}|<k$, we have

$\beta_{S_{1},S_{2}}:=|\\{m_{1}\in
B_{S_{1}}\colon\phi_{1}(m_{1})=S_{2}\\}|=|\\{m_{2}\in
C_{S_{2}}\colon\phi_{2}(m_{2})=S_{1}\\}|.$

  * •

For all $|S_{1}|,|S_{2}|<k$,

$\beta_{S_{1},S_{2}}=(1\pm\exp(-\Omega((\log n\log\log n)^{1/2})))n\prod_{p\in
S_{1}}p^{-1}\prod_{p\in
S_{2}}p^{-1}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus(S_{1}\cup
S_{2})}(1-2p^{-1})$

  * •

$|\\{m_{1}\colon\phi_{1}(m_{1})=\ast\\}|\leq 2\alpha^{1/2}n$.

  * •

$|\\{m_{2}\colon\phi_{2}(m_{2})=\ast\\}|\leq 2\alpha^{1/2}n$.

###### Proof.

We construct initial random assignment functions which essentially only fail
the third bullet point, and then “throw out” a small fraction of values by
reassigning them to $\ast$ in order to balance the two sides. For $m_{1}\in
B_{\geq k}$ deterministically assign $\phi_{1}^{\prime}(m_{1})=\ast$ and for
$m_{2}\in C_{\geq k}$ deterministically assign
$\phi_{2}^{\prime}(m_{2})=\ast$. Given $m_{1}\in B_{S}$ with $|S|<k$, for each
$|T|<k$ we assign $\phi_{1}^{\prime}(m_{1})=T$ with probability

$\mathbbm{1}_{S\cap T=\emptyset}\cdot\prod_{p\in
T}(p-1)^{-1}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus(S\cup
T)}(1-(p-1)^{-1})$

and otherwise set $\phi_{1}^{\prime}(m_{1})=\ast$. Analogously, given
$m_{2}\in C_{S}$ with $|S|<k$, for $|T|<k$ we assign
$\phi_{2}^{\prime}(m_{2})=T$ with probability

$\mathbbm{1}_{S\cap T=\emptyset}\cdot\prod_{p\in
T}(p-1)^{-1}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus(S\cup
T)}(1-(p-1)^{-1})$

and else set $\phi_{2}^{\prime}(m_{2})=\ast$. Evidently the stated
probabilities add to at most $1$ when summing over $|T|<k$, so this is well-
defined.

By Lemma 3.2 we have for $|S_{1}|,|S_{2}|<k$ which are disjoint that

$\displaystyle\mathbb{E}$ $\displaystyle|\\{m_{1}\in
B_{S_{1}}\colon\phi_{1}^{\prime}(m_{1})=S_{2}\\}|$
$\displaystyle=(1\pm\alpha^{2})n\prod_{p\in
S_{1}}p^{-1}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus
S_{1}}(1-p^{-1})\prod_{p\in
S_{2}}(p-1)^{-1}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus(S_{1}\cup
S_{2})}(1-(p-1)^{-1})$ $\displaystyle=(1\pm\alpha^{2})n\prod_{p\in
S_{1}}p^{-1}\prod_{p\in
S_{2}}p^{-1}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus(S_{1}\cup
S_{2})}(1-2p^{-1})$

and analogously

$\mathbb{E}|\\{m_{2}\in
C_{S_{2}}\colon\phi_{2}^{\prime}(m_{2})=S_{1}\\}|=(1\pm\alpha^{2})n\prod_{p\in
S_{1}}p^{-1}\prod_{p\in
S_{2}}p^{-1}\prod_{p\in\mathcal{P}_{W}^{\prime}\setminus(S_{1}\cup
S_{2})}(1-2p^{-1}).$

Applying Chernoff (Lemma 3.6) it follows with probability at least $1/2$ that
if for $|S_{1}|,|S_{2}|<k$ disjoint we set

$\beta_{S_{1},S_{2}}=\min(|\\{m_{1}\in
B_{S_{1}}\colon\phi_{1}^{\prime}(m_{1})=S_{2}\\}|,|\\{m_{2}\in
C_{S_{2}}\colon\phi_{2}^{\prime}(m_{2})=S_{1}\\}|)$

then

$\beta_{S_{1},S_{2}}=(1\pm 2\alpha^{2})n\prod_{p\in S_{1}}p^{-1}\prod_{p\in
S_{2}}p^{-1}\prod_{p\in\mathcal{P}_{W}\setminus(S_{1}\cup S_{2})}(1-2p^{-1}).$

To construct $\phi_{1}$, for each pair of disjoint $|S_{1}|,|S_{2}|$ of size
less than $k$ assign $\beta_{S_{1},S_{2}}$ values of $S_{1}$ to map to the
label $S_{2}$. Send all remaining values in $S_{1}$, if any, to $\ast$ (and
set $B_{\geq k}=\ast$ as before). For $\phi_{2}$, assign $\beta_{S_{1},S_{2}}$
values of $S_{2}$ to map to the label $S_{1}$ and send all remaining values in
$S_{2}$, if any, to $\ast$. This is well-defined since clearly

$\sum_{|S_{2}|<k}\beta_{S_{1},S_{2}}\leq\sum_{|S_{2}|<k}|\\{m_{1}\in
B_{S_{1}}\colon\phi_{1}^{\prime}(m_{1})=S_{2}\\}|\leq|B_{S_{1}}|$

and similar for $\phi_{2}^{\prime}$.

The first four properties follow immediately. The final two properties follow
from noting that

$\Big{|}|\\{m_{1}\in
B_{S_{1}}\colon\phi_{1}^{\prime}(m_{1})=S_{2}\\}|-|\\{m_{2}\in
C_{S_{2}}\colon\phi_{2}^{\prime}(m_{2})=S_{1}\\}|\Big{|}\leq
8\alpha^{2}\beta_{S_{1},S_{2}},$

so the extra amount of values assigned the label $\ast$ is not significant,
and by citing the second part of Lemma 3.2. ∎

To construct $(\sigma,R_{B},R_{C})$, we define $\beta_{S_{1},S_{2}}$ as in
Lemma 3.7, $\beta_{\ast}^{1}=|\\{m_{1}\colon\phi_{1}(m_{1})=\ast\\}|$, and
$\beta_{\ast}^{2}=|\\{m_{2}\colon\phi_{2}(m_{2})=\ast\\}|$. Next we consider
the following procedure:

  * •

Fix $\phi_{1}$, $\phi_{2}$ satisfying Lemma 3.7 with these values of $\beta$.

  * •

Choose $\sigma_{1}$ uniformly randomly among all permutations of $[n]$ with
$\sigma_{1}(B_{S})=B_{S}$ for all $S$ and $\sigma_{2}$ uniformly randomly
among all permutations of $[n]_{\mathrm{o}}$ with $\sigma_{2}(C_{S})=C_{S}$
for all $S$.

  * •

Let $B_{S_{1},S_{2}}=\\{m_{1}\in
B_{S_{1}}\colon\phi_{1}(\sigma_{1}(m_{1}))=S_{2}\\}$ and
$C_{S_{1},S_{2}}=\\{m_{2}\in
C_{S_{2}}\colon\phi_{2}(\sigma_{2}(m_{2}))=S_{1}\\}$ for $|S_{1}|,|S_{2}|<k$.
Furthermore let $R_{B}=\\{m_{1}\colon\phi_{1}(\sigma_{1}(m_{1}))=\ast\\}$ and
$R_{C}=\\{m_{2}\colon\phi_{2}(\sigma_{2}(m_{2}))=\ast\\}$

  * •

For disjoint $S_{1},S_{2}$ of size at most $k$ with
$(S_{1},S_{2})\neq(\emptyset,\emptyset)$, let $\sigma$ form a uniformly random
matching between $B_{S_{1},S_{2}}$ and $C_{S_{1},S_{2}}$ not violating the
coprimality condition.

  * •

For $(S_{1},S_{2})=(\emptyset,\emptyset)$ let $\sigma$ be a uniformly random
matching between $B_{\emptyset,\emptyset}\cup R_{B}$ and
$C_{\emptyset,\emptyset}\cup R_{C}$.

  * •

Output $(\sigma,R_{B},R_{C})$.

Note that this fully specifies $\sigma$ since if $S_{1}\cap
S_{2}\neq\emptyset$ we have $x\in B_{S_{1}}$ implies $\phi_{1}(x)\neq S_{2}$
and similar for $\phi_{2}$. Furthermore, given the output
$(\sigma,R_{B},R_{C})$ we can mostly reconstruct the procedure. In particular,
for $|S_{1}|,|S_{2}|<k$, $B_{S_{1},S_{2}}$ is the set of values in $B_{S_{1}}$
mapped to $C_{S_{2}}$ by $\sigma$ for all
$(S_{1},S_{2})\neq(\emptyset,\emptyset)$, while $B_{\emptyset,\emptyset}$ is
the set of values in $B_{\emptyset}\setminus R_{B}$ mapped to
$C_{\emptyset}\cup R_{C}$. This determines $\sigma_{1}$ up to right-
multiplication by some $\xi_{1}$ which preserves the sets $B_{S_{1},S_{2}}$
for all $|S_{1}|,|S_{2}|<k$ as well as $R_{B}$, and similar for $\sigma_{2}$.

Therefore, if $X$ is the number of possible outputs $(\sigma,R_{B},R_{C})$
where $\sigma$ is coprime and $Y$ is the number of ways to run the above
procedure where the output is coprime, we have

$Y\leq
X\cdot\beta_{\ast}^{1}!\beta_{\ast}^{2}!\prod_{\begin{subarray}{c}|S_{1}|,|S_{2}|<k\\\
S_{1}\cap S_{2}=\emptyset\end{subarray}}\beta_{S_{1},S_{2}}!^{2}\leq
C_{0}(n)\cdot
n^{\beta_{\ast}^{1}+\beta_{\ast}^{2}}\prod_{\begin{subarray}{c}|S_{1}|,|S_{2}|<k\\\
S_{1}\cap S_{2}=\emptyset\end{subarray}}\beta_{S_{1},S_{2}}!^{2}.$ (3.1)

Now we provide a lower bound on $Y$.

Note that by Lemma 3.4 and the Chernoff bound (Lemma 3.6), with probability at
least $3/4$ we have that for all disjoint $|S_{1}|,|S_{2}|<k$ the coprimality
bipartite graph between parts $B_{S_{1},S_{2}}$ and $C_{S_{1},S_{2}}$ for
$(S_{1},S_{2})\neq(\emptyset,\emptyset)$ has at most
$4W^{-1/2}\beta_{S_{1},S_{2}}$ edges missing from each vertex. Furthermore
with probability $3/4$ the non-coprimality bipartite graph between
$B_{\emptyset,\emptyset}\cup R_{B}$ and $C_{\emptyset,\emptyset}\cup R_{C}$
has at most $4W^{-1/2}\beta_{\emptyset,\emptyset}+|R_{B}|$ edges missing from
each vertex. Using that $|\beta_{\emptyset,\emptyset}|\geq n\exp(-2\log\log
n)$ and $|R_{B}|=\beta_{\ast}^{1}\leq 2\alpha^{1/2}n$, we see that the maximum
degree of the complement bipartite graph is bounded by
$(4W^{-1/2}+\alpha^{1/3})\beta_{\emptyset,\emptyset}$, say.

Therefore applying Lemma 3.5 shows that there are

$Y\geq(3/4)^{2}\beta_{\ast}^{1}!\beta_{\ast}^{2}!\prod_{|S|<k}|B_{S}|!\prod_{|S|<k}|C_{S}|!\cdot\prod_{\begin{subarray}{c}|S_{1}|,|S_{2}|<k\\\
S_{1}\cap
S_{2}=\emptyset\end{subarray}}(\beta_{S_{1},S_{2}}/e)^{\beta_{S_{1},S_{2}}}e^{-4(4W^{-1/2}+\alpha^{1/3})\beta_{S_{1},S_{2}}}$
(3.2)

possible valid ways to run the procedure. Now, given 3.1 and 3.2, a nearly
identical calculation as in Section 2 demonstrates

$C_{0}(n)\geq\exp(-n\alpha^{1/4})n!\bigg{(}\prod_{p>2}\frac{(p-1)^{2(1-1/p)}}{p\cdot(p-2)^{(1-2/p)}}\bigg{)}^{n}.$

Finally, using Lemma 3.1 implies the result.

## References

  * [1] Noga Alon, Vojtech Rödl, and Andrzej Ruciński, _Perfect matchings in $\epsilon$-regular graphs_, Electron. J. Combin. 5 (1998), Research Paper 13, 4.
  * [2] Tom Bohman and Fei Peng, _Coprime mappings and lonely runners_ , arXiv:2109.09860.
  * [3] Carlo Bonferroni, _Teoria statistica delle classi e calcolo delle probabilita_ , Pubblicazioni del R Istituto Superiore di Scienze Economiche e Commericiali di Firenze 8 (1936), 3–62.
  * [4] G. P. Egorychev, _The solution of the van der Waerden problem for permanents_ , Dokl. Akad. Nauk SSSR 258 (1981), 1041–1044.
  * [5] D. I. Falikman, _Proof of the van der Waerden conjecture on the permanent of a doubly stochastic matrix_ , Mat. Zametki 29 (1981), 931–938, 957.
  * [6] D. M. Jackson, _The combinatorial interpretation of the Jacobi identity from Lie algebras_ , J. Combinatorial Theory Ser. A 23 (1977), 233–256.
  * [7] Peter Keevash, _Counting designs_ , J. Eur. Math. Soc. (JEMS) 20 (2018), 903–927.
  * [8] Nathan McNew, _private communication_.
  * [9] Richard Montgomery, _Spanning trees in random graphs_ , Adv. Math. 356 (2019), 106793, 92.
  * [10] Carl Pomerance, _Coprime permutations_ , arXiv:2203.03085.
  * [11] Carl Pomerance, _Coprime matchings_ , Integers 22 (2022), Paper No. A2, 9.
  * [12] Carl Pomerance and J. L. Selfridge, _Proof of D. J. Newman’s coprime mapping conjecture_ , Mathematika 27 (1980), 69–83.
  * [13] Michael Simkin, _The number of $n$-queens configurations_, arXiv:2107.13460.
  * [14] N. J. A. Sloane, _The On-Line Encyclopedia of Integer Sequences (OEIS)_.
  * [15] Endre Szemerédi, _Is laziness paying off? (“Absorbing” method)_ , Colloquium De Giorgi 2010–2012, Colloquia, vol. 4, Ed. Norm., Pisa, 2013, pp. 17–34.