# Theory of spin center sensing of diffusion

Denis R. Candido<EMAIL_ADDRESS>Department of Physics and
Astronomy, University of Iowa, Iowa City, Iowa 52242, USA  Michael E. Flatté
<EMAIL_ADDRESS>Department of Physics and Astronomy, University of
Iowa, Iowa City, Iowa 52242, USA  Department of Applied Physics, Eindhoven
University of Technology, P.O. Box 513, 5600 MB, Eindhoven, The Netherlands

###### Abstract

Surface electric (charge) noise influences spin defects due to fluctuation of
the surface charge density and also the electrostatic potential at the crystal
surface. Surprisingly, the two-point correlation function of both the charged
particles’ positions and the surface electrostatic potential strongly
influences the power of the polynomial decay of the electric noise spectral
density; this power is not determined solely by the character of the charge
fluctuators. Time-dependent crossover behavior near the correlation time of
the fluctuators, of the spin defect’s relaxation and decoherence, provide a
quantitative fingerprint of the diffusive behavior of charged particles at the
surface.

Introduction.— Spin defects in solids are promising solid-state qubits whose
long spin coherence times permit optical initialization and readout [1, 2].
The energy levels of spin-1 defects lacking inversion symmetry, such as NV
centers in diamond, shift in response to electric fields to first order,
motivating their application to quantum sensing and metrology of charge
dynamics [3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15], including for
nanoscale probes of many-body physics [16, 17, 18, 19]. Undesirable
fluctuating electric and magnetic fields also couple to the spin[5, 20, 21,
22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 16], and their effects increase near
the surface due to trapped charges and imperfections. The differing spatial
dependence of fluctuating electric fields produced by uncorrelated point-like
and dipole charges at the surface, $E\propto|\textbf{r}|^{-2}$ and
$E\propto|\textbf{r}|^{-3}$, can produce a surface noise spectral density
$S_{E}(\omega)\propto d^{-2}$ and $S_{E}(\omega)\propto d^{-4}$ [23, 24, 25,
31], where $d$ is the defect depth; this noise at a resonant frequency for
spin transitions leads directly to spin relaxation and decoherence and reduces
both the sensitivity to external fields and the optical contrast in readout.
The rapid increase in noise near a surface conflicts with the desire to place
the spin as close to the surface as possible to increase sensitivity to a
signal originating from outside the material and to improve the sensor’s
spatial resolution. A major focus, therefore, for quantum sensing, explores
techniques to reduce these noise sources[23, 20, 21, 32, 33]; similar goals
help advance defect-based hybrid quantum systems [34, 35, 36, 37, 38, 39, 40,
41, 42, 43, 44, 45, 46], as the coupling of a spin to another excitation (e.g.
magnon, phonon) also improves the closer the spin comes to the surface. These
efforts depend on an accurate understanding of the noise and its depth-
dependent scaling.

Figure 1: (a) Electric surface noise (b) Electric noise spectral density
$S_{E}(\omega)$ as a function of the defect depth, $d$, for different
correlation functions of the electronic density. For both exponential
correlation functions $e^{-|r|/\lambda}$ and $e^{-|r|^{2}/\lambda^{2}}$ we
obtain a constant depth dependence for $d\ll d_{c}$, followed by a $d^{-2}$
dependence for $d\gg d_{c}$ ($d_{c}=\lambda$). For a correlation function
obeying the diffusion equation, we obtain $-\ln d$ depth dependence for $d\ll
d_{c}$ and $d^{-4}$ for $d\gg d_{c}$ [$d_{c}=({\cal D}/\omega)^{1/2}$ where
${\cal D}$ is the diffusion coefficient]. Within the gray box inset, we obtain
$S_{E}(\omega)\propto d^{-\gamma}$ with $\gamma=1,2$ and $3$.

Here we show that very different depth dependencies for the noise spectrum
$S_{E}(\omega)$ emerge when spatial correlation is considered; our general
formalism relies on the two-point correlation functions of the fluctuating
point-like charges’ surface density,
$\left\langle\delta\rho(\textbf{r}^{\prime},t^{\prime})\delta\rho(\textbf{r},t)\right\rangle$,
and the fluctuating electrostatic potential at the crystal surface,
$\left\langle\delta\phi(\textbf{r}^{\prime},t^{\prime})\delta\phi(\textbf{r},t)\right\rangle$,
yielding $S_{E}^{\delta\rho}$ and $S_{E}^{\delta\phi}$, respectively. For
example, we find $S_{E}^{\delta\rho}\propto 1/d^{4}$ instead of $1/d^{2}$ for
fluctuation of the point-like particle’s density arising from diffusive
phenomena. We find $S_{E}^{\delta\rho}$ is independent of depth for
fluctuating charge density with an exponential spatial correlation function
(such as arises from crystal surface distortions and imperfections[47, 48])
and a spin depth much shorter than the correlation length. We further obtain
$S_{E}^{\delta\phi}\propto 1/d^{4}$ for exponential correlation of the
fluctuation of the surface potential. In general, the depth-dependent behavior
of $S_{E}^{\delta\rho}$ and $S_{E}^{\delta\phi}$ strongly depends on the
detailed physics of the surface noise phenomenon and the corresponding two-
point correlation functions of the fluctuating charge density and fluctuating
surface potential. $S_{E}^{\delta\rho}$ and $S_{E}^{\delta\phi}$ even exhibit
different depth dependence for the same correlation functions describing
fluctuating density and surface potential. These results suggest a new sensing
approach; correlating surface-dependent physical phenomena with depth-
dependent or frequency-dependent noise spectra in quantum sensing and
metrology. A specific example is the use of spin defects to sense diffusive
behavior, including charges within crystals, at the crystal surface, or even
in a liquid surrounding the crystal; diffusion constants and correlation times
can be extracted from the time-dependent spin coherence as a function of
frequency.

Fundamentals of electric noise at the diamond surface.— Hydrogen termination
generates an acceptor level; oxygen termination generates a donor level [49,
50, 51, 52]; similar effects emerge from imperfections of crystal termination
[53] along with spatial fluctuations of the surface potential [54, 55].
Interaction of these levels with trapped charges in the nitrogen doped
diamonds that host NV centers results in a creation of an effective surface 2D
hole or electron gas [56, 57, 58, 50, 59, 53]. These charges do not distribute
uniformly and move around due to the thermal fluctuations of the charge’s
position, collision between charges, capture and ionization by donors and
acceptors, and others processes. Moreover, as the read-out and initialization
of the sensing spin requires laser illumination, the light increases these
fluctuations. Thus the fluctuation of charged particles and the local surface
potential produce noise that is experienced by shallow spins, causing
decoherence of the spin state.

The temporal fluctuation of these charged particles relative to equilibrium,
and the surface potential, produce a temporally fluctuating electric field
$\textbf{E}(\textbf{r},t)=\left\langle\textbf{E}(\textbf{r})\right\rangle+\delta\textbf{E}(\textbf{r},t)$,
which causes decoherence [5, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31,
16], increases the photoluminescence linewidth [60, 32, 61, 33]. For Dirichlet
boundary conditions

${\textbf{E}\left(\textbf{r},t\right)=\frac{1}{\varepsilon}\int_{\mathcal{V}}d\textbf{r}^{\prime}\textbf{K}\left(\textbf{r}-\textbf{r}^{\prime}\right)\rho\left(\textbf{r}^{\prime},t\right)+\oint_{\mathcal{S}}da^{\prime}\phi\left(\textbf{r}^{\prime}_{S},t\right)\frac{\partial\textbf{K}\left(\textbf{r}-\textbf{r}^{\prime}_{S}\right)}{\partial
n^{\prime}}},$ (1)

and

$\displaystyle\left\langle\delta{E}_{\nu}\left(\textbf{r},t\right)\delta{E}_{\mu}\left({\textbf{r}},0\right)\right\rangle$
$\displaystyle=\left(\frac{1}{\varepsilon}\right)^{2}\int_{\mathcal{V}}\int_{\mathcal{V}}d\textbf{r}^{\prime}d\textbf{r}^{\prime\prime}{K}_{\mu}\left(\textbf{r}-\textbf{r}^{\prime}\right){K}_{\nu}\left({\textbf{r}}-\textbf{r}^{\prime\prime}\right)\left\langle\delta\rho\left(\textbf{r}^{\prime},t\right)\delta\rho\left(\textbf{r}^{\prime\prime},0\right)\right\rangle$
$\displaystyle+\int_{\mathcal{S}}\int_{\mathcal{S}}d\textbf{r}^{\prime}_{\parallel}d\textbf{r}^{\prime\prime}_{\parallel}\frac{\partial{K}_{\mu}(\textbf{r}-\textbf{r}^{\prime}_{\parallel})}{\partial
z}\frac{\partial{K}_{\nu}({\textbf{r}}-\textbf{r}^{\prime\prime}_{\parallel})}{\partial
z}\expectationvalue{\delta\phi(\textbf{r}^{\prime}_{\parallel},t)\delta\phi(\textbf{r}^{\prime\prime}_{\parallel},0)}$
(2)

for $\textbf{r}_{S}$ on ${\cal S}=\left\\{(x,y,z),\thinspace-L/2\leq x,y\leq
L/2,\thinspace\thinspace z=0\right\\}$, r on ${\cal
V}=\left\\{(x,y,z),\thinspace-L/2\leq x,y\leq L/2,\thinspace\thinspace
z>0\right\\}$, and $n^{\prime}$ the normal to this surface. The correlation of
the fluctuating electric field [Eq. (2)] is determined by the correlation of
the fluctuating surface electrostatic potential (as derived in Refs. 62 and
63) but also the fluctuating surface charge density.

To obtain the central equation for two-dimensional surface noise, the Kernel
and its derivative with respect to $z$ are expressed using a 2D Fourier
transformation for the spatial directions parallel to the surface (yielding
the Weyl representation of the Green’s function)

${\textbf{K}(\textbf{r}-\textbf{r}^{\prime}_{\parallel})=\left(2\pi\right)^{-2}\int{d\textbf{k}_{\parallel}}{\cal
F}_{2D}^{\textbf{k}_{\parallel}}[\textbf{K}(\textbf{r}-\textbf{r}^{\prime}_{\parallel})]e^{i\textbf{k}_{\parallel}\cdot\left(\textbf{r}{}_{\parallel}-\textbf{r}^{\prime}_{\parallel}\right)}},$
(3) ${{\cal
F}_{2D}^{\textbf{k}_{\parallel}}[K_{\mu}(\textbf{r})]=\frac{1}{2}\left[i\left(\delta_{\nu,x}+\delta_{\nu,y}\right)k_{\nu}/|\textbf{k}_{\parallel}|-\delta_{\nu,z}\right]e^{-|\textbf{k}_{\parallel}|\left|z\right|}},$
(4) ${\cal
F}_{2D}^{\textbf{k}_{\parallel}}\left[\partial_{z}K_{\mu}\left(\textbf{r}\right)\right]=\frac{1}{2}\left[i\left(\delta_{\nu,x}+\delta_{\nu,y}\right)k_{\nu}-\delta_{\nu,z}\left|\textbf{k}_{\parallel}\right|\right]e^{-\left|\textbf{k}_{\parallel}\right|\left|z\right|}.$
(5)

We restrict fluctuations of density and surface potential to occur only on
$\cal{S}$, thus
$\rho\left(\textbf{r},t\right)=Q\delta(z)n\left(\textbf{r}_{\parallel},t\right)$.
We assume transitional symmetry of the correlation function for the density
and surface potential along the surface,

$\left\langle\delta n(\textbf{r}^{\prime}_{\parallel},t)\delta
n(\textbf{r}_{\parallel},0)\right\rangle=\left\langle\delta
n(\textbf{r}^{\prime}_{\parallel}-\textbf{r}_{\parallel},t)\delta
n\left(0,0\right)\right\rangle$ (6)
${\left\langle\delta\phi(\textbf{r}^{\prime}_{\parallel},t)\delta\phi(\textbf{r}^{\prime\prime}_{\parallel},0)\right\rangle}=\left\langle\delta\phi(\textbf{r}^{\prime}_{\parallel}-\textbf{r}^{\prime\prime}_{\parallel},t)\delta\phi(0,0)\right\rangle$
(7)

Finally, ${S_{E_{\nu}}(\omega)=\int_{-\infty}^{\infty}dt\left\langle\delta
E_{\nu}\left(\textbf{r},t\right)\delta
E_{\nu}\left(\textbf{r},0\right)\right\rangle e^{-i\omega t}}$, yielding

$\displaystyle S_{E_{\nu}}^{\delta\rho}\left(\omega\right)$
$\displaystyle=\left(\frac{Q}{\varepsilon}\right)^{2}\int\frac{d\textbf{k}_{\parallel}}{\left(2\pi\right)^{2}}{\cal
F}_{2D}^{\textbf{k}_{\parallel}}\left[K_{\nu}\left(\textbf{r}\right)\right]{\cal
F}_{2D}^{-\textbf{k}_{\parallel}}\left[K_{\nu}\left(\textbf{r}\right)\right]{\cal
F}_{2D}^{\textbf{k}_{\parallel},\omega}[\left\langle\delta
n(\textbf{r}^{\prime}_{\parallel}-\textbf{r}_{\parallel},t)\delta
n\left(0,0\right)\right\rangle]$ (8) $\displaystyle
S_{E_{\nu}}^{\delta\phi}\left(\omega\right)$
$\displaystyle=\int\frac{d\textbf{k}_{\parallel}}{\left(2\pi\right)^{2}}{\cal
F}_{2D}^{\textbf{k}_{\parallel}}[\partial_{z}K_{\nu}(\textbf{r})]{\cal
F}_{2D}^{-{\textbf{k}}_{\parallel}}[\partial_{z}K_{\nu}(\textbf{r})]{\cal
F}_{2D}^{\textbf{k}_{\parallel},\omega}[\left\langle\delta\phi(\textbf{r}^{\prime}_{\parallel}-\textbf{r}^{\prime\prime}_{\parallel},t)\delta\phi(0,0)\right\rangle]$
(9)

Equations (3)-(9) establish the relation between electric spectral noise and
the time and space Fourier transform of the correlation of the surface density
and potential at different times and positions.

Fluctuating surface density from diffusion — We assume $\left\langle
n\left(\textbf{r}_{\parallel},t\right)\right\rangle$ satisfies

$\left(\frac{\partial}{\partial
t}+\frac{1}{\tau}+\mu\textbf{E}_{\parallel}\cdot\nabla_{\parallel}-{\cal
D}\nabla^{2}_{\parallel}\right)\left\langle
n\left(\textbf{r}_{\parallel},t\right)\right\rangle=0,$ (10)

which includes diffusion, drift and a finite carrier lifetime. Accordingly,
the Green’s function satisfies

${(\frac{\partial}{\partial
t}+\frac{1}{\tau}+\mu\textbf{E}_{\parallel}\cdot\nabla_{\parallel}-{\cal
D}\nabla^{2}_{\parallel}){\cal
G}\left(\textbf{r}_{\parallel},t\right)=\delta(\textbf{r}_{\parallel})\delta(t)}$
(11)

and ${\left\langle\delta n(\textbf{r}_{\parallel},t)\right\rangle=\int
d\textbf{r}^{\prime}{\cal
G}(\textbf{r}_{\parallel}-\textbf{r}^{\prime}_{\parallel},t)\left\langle\delta
n(\textbf{r}^{\prime}_{\parallel},0)\right\rangle}$, so

$\left\langle\delta
n(\textbf{r}^{\prime}_{\parallel}-\textbf{r}_{\parallel},t)\delta
n(0,0)\right\rangle=\int d\textbf{r}^{\prime\prime}_{\parallel}{\cal
G}(\textbf{r}^{\prime\prime}_{\parallel}-\textbf{r}^{\prime}_{\parallel},t)\left\langle\delta
n(\textbf{r}^{\prime\prime}_{\parallel}-\textbf{r}_{\parallel},0)\delta
n(0,0)\right\rangle$ (12)

With this, Eqs. (8) show that the spectral density noise is determined as
function of Fourier transform of both the electrostatic Green’s function and
the Green’s function related to the differential equation describing the
dynamics of our charged particles. Using ${\left\langle\delta
n(\textbf{r}^{\prime}_{\parallel},0)\delta
n(\textbf{r}_{\parallel},0)\right\rangle=\delta(\textbf{r}^{\prime}_{\parallel}-\textbf{r}_{\parallel})\left\langle\delta
n(\textbf{r}_{\parallel},0)\right\rangle}$, we obtain $\left\langle\delta
n(\textbf{r}^{\prime}_{\parallel},t)\delta
n\left(\textbf{r}_{\parallel},0\right)\right\rangle={\cal
G}(\textbf{r}^{\prime}_{\parallel}-\textbf{r}_{\parallel},t)\left\langle
n\left(\textbf{r}_{\parallel},0\right)\right\rangle$ and ${\cal
F}_{2D}^{\textbf{k}_{\parallel},\omega}\left[\left\langle\delta
n\left(\textbf{r}{}_{\parallel},t\right)\delta
n\left(0,0\right)\right\rangle\right]={\cal
F}_{2D}^{\textbf{k}_{\parallel},\omega}\left[{\cal
G}\left(\textbf{r}{}_{\parallel},t\right)\right]\left\langle
n\left(0,0\right)\right\rangle$ with

${\cal F}_{2D}^{\textbf{k}_{\parallel},\omega}[{\cal
G}(\textbf{r}_{\parallel},t)]=\frac{1}{{i\omega+1/\tau+i\mu\textbf{E}_{\parallel}\cdot\textbf{k}_{\parallel}+{\cal
D}\textbf{k}_{\parallel}^{2}}},$ (13)

following from Eq. (10).

Figure 2: (a) $S_{E}^{\delta\rho}(\omega)$ [Eq. (14)] as a function of
$d/d_{c}$ ($d_{c}=\sqrt{{\cal D}/\omega}$). The dashed and dotted lines are
guides for the eyes for $1/d^{\gamma}$ dependencies with $\gamma=1,2,3$ and
$4$, and the green circles gives the depths in which $\gamma$-dependence
changes. (b) Same as (a) for different ${\cal D}=1,10,100$ (nm2 GHz) and $d=1$
nm. Inset shows the $S_{E}^{\delta\rho}(\omega)$ [Eq. (14)] dependence on
${\cal D}$ for $d=1$ nm and $\omega=1$ GHz. (c) $S_{E}^{\delta\rho}(\omega)$
[Eq. (14)] as a function $\omega/\omega_{c}$ ($\omega_{c}={\cal D}/d^{2}$),
where dotted and dashed correspond to the asymptotic analytical forms within
Eq. (15). (d) $S_{E}^{\delta\rho}(\omega)$ [Eq. (18)] as a function of $d$ for
$\lambda=0.1$ and $1$ nm, with dotted, dot-dashed and dashed lines
corresponding to the asymptotic analytical forms [Eq. (20)]. (e) Same as (d)
for $S_{E}^{\delta\phi}(\omega)$, Eq. (22), with inset showing the $d$
dependencies for different $\Lambda$ values. (f) Same as (e) for
$S_{E}^{\delta\phi}(\omega)$, Eq. (23).

We now use this equation within Eq. (9) assuming first
$\textbf{E}_{\parallel}=0$ and $1/\tau\rightarrow 0$, i.e., purely diffusive
regime, which for $\mu=x,y$ yields

$\displaystyle
S_{E_{\mu}}^{\delta\rho}\left(\omega\right)=\left(\frac{Q}{4\pi\varepsilon}\right)^{2}{n_{{\cal
S}}\pi}\times\frac{2{\cal
D}}{\omega^{2}d^{4}}\int_{0}^{\infty}d\epsilon\thinspace\frac{\epsilon^{3}e^{-2\epsilon}}{1+\frac{{\cal
D}^{2}}{\omega^{2}d^{4}}\epsilon^{4}},$ (14)

where $\epsilon$ is a dimensionless variable and ${S_{E_{z}}=2S_{E_{x}}\equiv
S_{E}}$. It is important to notice that we cannot decouple the temporal part
of the spectral noise from the spatial part. However, we see there are two
limit behaviors, dictated by the value of $d^{2}\omega/{\cal D}$, in which the
Eq. (14) can be solved analytically, yielding,

$\displaystyle
S_{E_{x}}^{\delta\rho}\left(\omega\right)\approx\begin{cases}{B_{0}}\frac{3{\cal
D}}{4d^{4}\omega^{2}}&\thinspace d^{2}\omega/{\cal D}\gg 1\\\
B_{0}\frac{1}{{\cal D}}\left[\ln\left(\frac{{\cal D}}{4\omega
d^{2}}\right)+2\Gamma\right]&d^{2}\omega/{\cal D}\ll 1\end{cases}$ (15)

with $B_{0}=\left({Q}/{4\pi\varepsilon}\right)^{2}n_{S}\pi$ and the Euler-
Gamma number $\Gamma=0.577$ . In Fig. 1 (a) we plot $S_{E}^{\delta\rho}$ [Eq.
(14)] as a function of $d/d_{c}$ with $d_{c}=({{\cal D}/\omega})^{1/2}$. We
see that while $S_{E}^{\delta\rho}$ behaves as $d^{-4}$ for $d\gg d_{c}$, we
obtain a $-\ln{d}$ dependence for $d\ll d_{c}$, with the intermediate region
$d\approx d_{c}$ yielding $d^{-\gamma}$ (with $\gamma=1,2$ and $3$) shown in
Fig. 2(a).

We emphasize that a $d^{-4}$ depth dependence on defect coherence time
reported in the literature was attributed to the dipole character of the
fluctuators [29, 30, 31, 64]. Our results show that this depth dependence can
also be obtained from point-like charges that diffuse. Diffusive motion of
charges occurs on crystal surfaces in which crystal imperfections create traps
in which the charge can be trapped and released. Inside the gray region of
Fig. 1(a), i.e., for $0.1\lesssim d/d_{c}\lesssim 10$ we see a $d^{-\gamma}$
dependence with $\gamma=1,2$ and $3$, as shown in Fig. 2(a). In Fig. 2(b) we
plot $S_{E}$ as a function of $d$ for different diffusion coefficient. We note
that for depth $d\ll 1$ nm, faster diffusion (larger $\cal{D}$) suppresses
noise, similar to the effect of motional narrowing [65, 66]. However for $d\gg
1$ nm. we obtain the opposite behavior, and a plot of $S_{E}^{\delta\rho}$ as
a function of ${\cal D}$ shows a maximum at ${\cal D}\approx 1$ nm2 GHz for
$\omega=1$ GHz and $d=1$ nm [Inset Fig. 2(b)]. We plot in Fig. 2(c)
$S_{E}^{\delta\rho}$ as a function of $\omega$ where we see the crossover
between $1/\omega^{2}$ and $-\ln{\omega}$ taking place at $\omega_{c}={\cal
D}/d^{2}$.

Fluctuating density with exponential spatial correlation function.— Here we
assume that the correlation function of the fluctuating charge density is
given by an exponential, i.e., $\left\langle\delta
n(\textbf{r}^{\prime}_{\parallel},t)\delta
n(\textbf{r}_{\parallel},0)\right\rangle=\chi\left(t\right)(n_{S}/A)e^{-{\left|\textbf{r}^{\prime}_{\parallel}-\textbf{r}_{\parallel}\right|^{n}}/\lambda^{n}}$
with $n>0\in{\mathbb{Z}}$ and $A$ the surface area. This assumption is very
common in crystals with inhomogeneous surfaces, and it can be obtained
theoretically assuming random crystal surface distortions and imperfections
[47, 48]. This behavior has also being experimentally realized through X-ray
scattering intensity of the crystal surfaces [47, 48, 67, 68], which directly
measures the degree of correlation between the fluctuation of the electronic
density at two different points as a function of its distance [47, 48, 67,
68]. Accordingly,

${{\cal F}_{2D}^{\textbf{k}_{\parallel}}\left[\left\langle\delta
n\left(\textbf{r}{}_{\parallel},t\right)\delta
n\left(0,0\right)\right\rangle\right]=\chi\left(t\right)(n_{S}/A){2\pi\lambda^{2}}/{(1+\lambda^{2}\textbf{k}_{\parallel}^{2})^{3/2}}}$
(16) ${{\cal F}_{2D}^{\textbf{k}_{\parallel}}\left[\left\langle\delta
n\left(\textbf{r}{}_{\parallel},t\right)\delta
n\left(0,0\right)\right\rangle\right]=\chi\left(t\right)(n_{S}/A)\pi\lambda^{2}e^{-{\lambda^{2}\textbf{k}_{\parallel}^{2}}/{4}}}$
(17)

for $n=1$ and $n=2$, respectively. This yields through Eq. (8)

$\displaystyle{S_{E_{x}}^{\delta\phi,n=1}\left(\omega\right)}{}$
$\displaystyle=\left(\frac{Q}{\varepsilon}\right)^{2}\frac{n_{S}}{8A}\frac{{\cal
F}_{2D}^{\omega}\left[\chi\left(t\right)\right]}{\left(d/\lambda\right)^{2}}\int_{0}^{\infty}d\epsilon\frac{\epsilon
e^{-2\epsilon}}{[1+(\lambda/d)^{2}\epsilon^{2}]^{3/2}},$ (18)
$\displaystyle{S_{E_{x}}^{\delta\phi,n=2}\left(\omega\right)}{}$
$\displaystyle=\left(\frac{Q}{\varepsilon}\right)^{2}\frac{n_{S}}{8A}\frac{{\cal
F}_{2D}^{\omega}\left[\chi\left(t\right)\right]}{\left(d/\lambda\right)^{2}}\int_{0}^{\infty}d\epsilon\thinspace\epsilon
e^{-2\epsilon-\frac{\epsilon^{2}}{4(d/\lambda)^{2}}}.$ (19)

Similarly to the diffusion case, here both exponential correlation functions
contain a critical length defined by $d_{c}=\lambda$, and the corresponding
asymptotic behaviors follows

$\displaystyle
S_{E_{x}}^{\delta\phi,n=1}\left(\omega\right)=\begin{cases}B_{1}{\cal
F}_{2D}^{\omega}\left[\chi\left(t\right)\right]&d\ll d_{c},\\\ B_{1}{\cal
F}_{2D}^{\omega}\left[\chi\left(t\right)\right]\frac{1}{4\left(d/d_{c}\right)^{2}}&\thinspace
d\gg d_{c},\\\ \end{cases}$ (20) $\displaystyle
S_{E_{x}}^{\delta\phi,n=2}\left(\omega\right)=\begin{cases}B_{1}{\cal
F}_{2D}^{\omega}\left[\chi\left(t\right)\right]&d\ll d_{c},\\\ (1/2)B_{1}{\cal
F}_{2D}^{\omega}\left[\chi\left(t\right)\right]\frac{1}{4\left(d/d_{c}\right)^{2}}&\thinspace
d\gg d_{c},\end{cases}$ (21)

with $B_{1}=\left({Q}/{\varepsilon}\right)^{2}n_{S}/(8A)$. We see that both
exponential correlation functions produce similar asymptotic behaviors, with a
corresponding $d^{-2}$ depth dependence for $d\gg d_{c}$. Compared to Eq.
(15), Eqs. (20) and (21) show that different classes of density spatial
correlation function yields different depth dependence. In Fig. 2(d) we plot
Eq. (18) as a function of $d$ for $\lambda=0.1$ and $1$ nm, together with its
asymptotic functions given by Eq. (20). Surprisingly, there is a plateau of
$S_{E}^{\delta\rho}$ for depth much shorter than $d_{c}=\lambda$, with nominal
value independent on the power of the exponential correlation [See Fig. 1(b)].
Additionally, we also see that for intermediate depth values, $d\approx
d_{c}$, $S_{E}^{\delta\phi}$ scales with $1/d$.

Effect of the fluctuating surface potential.— Differently from the charged
particles, the fluctuation of the electrostatic potential cannot be described
by the diffusive model Eq. (10). Alternatively, we assume our surface composed
of various plaquettes [55, 62] with corresponding local uniform potential,
$\phi_{i}\left(t\right)$ and $\left\langle\phi_{i}(t)\right\rangle=\phi_{eq}$.
Moreover, we assume
$\delta\phi\left(\textbf{r}_{\parallel},t\right)=\sum_{{i=1}}^{N}\delta\phi_{i}\left(t\right)\Theta_{i}(\textbf{r}_{\parallel})$
where $\Theta_{i}(\textbf{r}_{\parallel})$ is the square function associated
to the plaquette $i$, i.e., it assumes value $1$ $(0)$ inside (outside)
plaquette $i$, centered at $\textbf{r}_{\parallel}$ [55, 62], and with
$\int_{\cal
S}d\textbf{r}_{\parallel}\sum{}_{i}\Theta_{i}(\textbf{r}_{\parallel})=A$.
Using
${\left\langle\delta\phi_{i}\left(t\right)\delta\phi_{j}\left(0\right)\right\rangle=\delta_{ij}f_{\phi}\left(t\right)}$,
we can write
$\sum{}_{i}\left\langle\Theta_{i}(\textbf{r}_{\parallel})\Theta_{j}(\textbf{r}^{\prime}_{\parallel})\right\rangle=N\left\langle\Theta_{i}(\textbf{r}_{\parallel})\Theta_{i}(\textbf{r}^{\prime}_{\parallel})\right\rangle$
due to the randomness of $\Theta$, where $N$ is the number of plaquettes.
Assuming the spatial transitional symmetry of the correlation function of
plaquettes’ positions, we model
$\left\langle\Theta_{i}(\textbf{r}_{\parallel})\Theta_{i}(\textbf{r}^{\prime}_{\parallel})\right\rangle=e^{-|\textbf{r}_{\parallel}-\textbf{r}^{\prime}_{\parallel}|^{n}/\Lambda^{n}}$
[62, 63]. For $n=1$ and $n=2$, we obtain

$\displaystyle S_{E_{x}}^{\delta\phi,n=1}\left(\omega\right)$
$\displaystyle={\cal
F}_{2D}^{\omega}\left[f_{\phi}\left(t\right)\right]N\frac{\Lambda^{2}}{8d^{4}}\int_{0}^{\infty}d\epsilon\frac{\epsilon^{3}e^{-2\epsilon}}{[1+(\Lambda/d)^{2}\epsilon^{2}]^{3/2}},$
(22) $\displaystyle S_{E_{x}}^{n=2}\left(\omega\right)$ $\displaystyle={\cal
F}_{2D}^{\omega}\left[f_{\phi}\left(t\right)\right]N\frac{\Lambda^{2}}{4d^{4}}\int_{0}^{\infty}d\epsilon\epsilon^{3}e^{-2\epsilon}e^{-\left(\Lambda/d\right)^{2}\epsilon^{2}/4},$
(23)

where Eq. (22) has been already reported previously [62]. Similarly to the
previous cases, the above integrals exhibits asymptotic analytical behavior
given by

$\displaystyle S_{E_{x}}^{\delta\phi,n=1}\left(\omega\right)$
$\displaystyle=\begin{cases}{\cal
F}_{2D}^{\omega}\left[f_{\phi}\left(t\right)\right]\frac{N}{8}\frac{1}{2d\Lambda}&d\ll\Lambda,\\\
{\cal
F}_{2D}^{\omega}\left[f_{\phi}\left(t\right)\right]\frac{N}{8}\frac{3\Lambda^{2}}{8d^{4}}&d\gg\Lambda,\end{cases}$
(24) $\displaystyle\vspace{5mm}S_{E_{x}}^{n=2}\left(\omega\right)$
$\displaystyle=\begin{cases}{\cal
F}_{2D}^{\omega}\left[f_{\phi}\left(t\right)\right]\frac{N}{8}\frac{8}{\Lambda^{2}}&d\ll\Lambda,\\\
{\cal
F}_{2D}^{\omega}\left[f_{\phi}\left(t\right)\right]\frac{N}{8}\frac{3\Lambda^{2}}{8d^{4}}&d\gg\Lambda,\end{cases}$
(25)

from where we obtain a crossover between $d^{-1}$ ($d^{0}$) and $d^{-4}$ depth
dependencies for $n=1$ ($n=2$) as can be seen in Fig. 2(e) [Fig. 2(f)]. This
result also corroborate the fact that the dependence $d^{-4}$ cannot be
attributed solely to the electric noise arising from fluctuating dipole
charges. Additionally, this also shows that the $d$-dependence for short
characteristic correlation length ($\Lambda$) depends strongly on the
correlation form, i.e., $e^{-|r|/\Lambda}$ or $e^{-|r|^{2}/\Lambda^{2}}$. We
also see in Figs. 2(e) and (f) a $d^{-1}$, $d^{-2}$ and $d^{-3}$ dependencies
for intermediate values of $d$. Finally, for $n=1$ and $n=2$ we see that the
bigger (smaller) the $\Lambda$ the stronger (weaker) the suppression of the
noise for $d\ll\Lambda$ ($d\gg\Lambda$).

Figure 3: Plot of $\left\langle\phi^{2}\right\rangle=-\ln{\cal C}(t)$ [Eq.
(26)] as a function $t/\tau_{c}$. The solid line represent the exact numeric
evaluation of Eq. (26), while the dashed and dot-dashed represent its
asymptotic behaviors shown by Eq. (27).

Sensing the diffusion with spin defects.— Our theory allows for the probe and
sensing of the diffusion regime using spin defects, together with the
experimental extraction of the corresponding ${\cal D}$. We envision two
different possible ways to do that. First, we can have different defects
located at different depths. If there is a diffusion regime happening, the
crossover from $d^{-4}\rightarrow d^{-3}\rightarrow d^{-2}\rightarrow
d^{-1}\rightarrow-\ln d$ should be seen through either $T_{1}^{-1}$ or
$T_{2}^{-1}$, with experimental results showing a trend similar to Figs. 1 and
2(a) as $T_{1,2}^{-1}\propto S_{E}(\omega)$. By identifying the critical depth
in which the crossover happens ($d_{c}$), we can extract the diffusion
coefficient through ${\cal D}=d_{c}^{2}{\omega_{c}}$. Secondly, we could also
verify the diffusion regime with only one defect located at $d$. For this
case, we have to plot either $1/T_{1}$ or $1/T_{2}$ as a function of $\omega$.
Here, the experimental results should show similar trend to the one in Fig.
2(c), with a crossover located at $\omega_{c}$. Similarly, the extraction of
the diffusion coefficient can be assessed through ${\cal
D}=d^{2}{\omega_{c}}$.

Similarly to the second case, and in analogy with recent works on the use of
spin defects to sense many-body interactions [18, 19], here we plot the
average probe coherence, ${\cal C}(t)$, which characterizes the time evolution
of both the dephasing or relaxation processes due to the electric noise.
Accordingly, we have for Gaussian noise ${\cal
C}(t)=e^{-(\delta^{2}/2)\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}\left\langle
E_{\mu}(t^{\prime\prime})E_{\mu}(t^{\prime})\right\rangle
e^{i\Omega(t^{\prime\prime}-t^{\prime})}}$ [33, 31] where $\delta$ is
proportional to the spin-defect electric dipole moment. $\Omega$ characterizes
either spin relaxation, with its value given by the corresponding frequency
difference of the levels; or spin decoherence for $\Omega=0$ [33, 64, 31]. By
rewriting this equation using the electric noise spectral density, we obtain

${\cal C}(t)={\rm
exp}\left\\{-\frac{\delta^{2}}{2}\int_{-\infty}^{\infty}\frac{d\omega}{2\pi}S_{E}\left(\omega\right)\frac{\sin^{2}\left[\left(\omega+\Omega\right)t/2\right]}{\left[\left(\omega+\Omega\right)/2\right]^{2}}\right\\},$
(26)

with corresponding limits ${\cal
C}(t)\approx\text{exp}\left[-\frac{\delta^{2}}{2}S_{E}(t=0)t^{2}\right]$ with
$S_{E}(t=0)=\int_{-\infty}^{\infty}d\omega/(2\pi)S_{{E}}(\omega)$ for
$t\ll\tau_{c}$, and ${\cal
C}(t)\approx\text{exp}\left[-\frac{\delta^{2}}{2}S_{{E}}(\omega=-\Omega)t\right]$
for $t\gg\tau_{c}$. Accordingly, for the diffusion case, we obtain the
following corresponding asymptotic behaviors 111As $\lim_{\omega\rightarrow
0}S_{E}^{\delta\rho}(\omega)\rightarrow\infty$, this asymptotic behavior only
holds for $\Omega\neq 0$. For $\Omega=0$ we still find asymptotic behaviors,
but they are obtained by first performing the integral of Eq. (26) and then
taking the limit $\Omega\rightarrow 0$.

${\cal C}(t)=\begin{cases}{\rm
exp}\left[{-\frac{\delta^{2}}{2}\frac{B_{0}}{4d^{2}}t^{2}}\right]&\quad
t\ll\tau_{c}\\\ \vspace{-3mm}\\\ {\rm
exp}\left[-\frac{\delta^{2}}{2}S_{E}^{\delta\rho}(-\Omega)t\right]&\quad
t\gg\tau_{c}\end{cases},$ (27)

which indicates a crossover between two different exponential decays, namely,
$e^{-At^{2}}$ and $e^{-Bt}$. We therefore claim that a diffusion leading
decoherence and decay can be demonstrated and observed experimentally. This
makes it possible to determine the correlation time of diffuse mediums through
ordinary defect sensing techniques as $\tau_{c}=\tau_{c}({\cal D})$.
Additionally, we obtain a depth dependence $\propto d^{-2}$ only for the
condition $t\ll\tau_{c}$. For $t\gg\tau_{c}$, we obtain a depth condition
dictated by the characteristic of our modeling, determined by Eq. (15).

Our framework predicts a non-trivial temporal profile for the average
coherence of probe spins, which exhibits a crossover between two distinct
stretched exponential decays [Fig. 1] [10–16, 32]. Crucially, we demonstrate
theoretically, with experimental support, that the associated stretch powers
contain a wealth of information about both the static and dynamical properties
of the many-body spin system. We focus on three particular properties. First,
the stretch power can distinguish between different forms of spectral
diffusion, shedding light on the nature of local spin fluctuations. Second,
the crossover in time between different stretch powers allows one to extract
the many-body system’s correlation time. Finally, the stretch power also
contains direct signatures about the dimensionality and disorder intrinsic to
the system.

Sensing shot noise with spin defects.— Our approach also permits sensing shot
noise phenomena with spin defects. Since the noise of the current is also
produce by the fluctuation of the charged particles’ density when applied an
electric field, i.e., $\left\langle\delta j_{\nu}(t^{\prime})\delta
j_{\nu}(t)\right\rangle\propto e^{2}\left\langle\delta n(t^{\prime})\delta
n(t)\right\rangle$, a clear relation between the shot noise and $T_{1,2}$ can
be established, although this is beyond the scope of this work.

###### Acknowledgements.

We thank D. D. Awschalom, P. C. Maurer, M. Fukami, and J. C. Karsch for useful
discussions. The work was supported as part of the Center for Molecular
Quantum Transduction, an Energy Frontier Research Center funded by the U.S.
Department of Energy, Office of Science, Basic Energy Sciences, under Award
Number DE-SC0021314.

## References

  * Awschalom  _et al._ [2002] D. D. Awschalom, D. Loss, and N. Samarth, eds.,  _Semiconductor spintronics and quantum computation_ (Springer Science & Business Media, 2002).
  * Anderson  _et al._ [2021] C. P. Anderson, E. O. Glen, C. Zeledon, A. Bourassa, Y. Jin, Y. Zhu, C. Vorwerk, A. L. Crook, H. Abe, J. Ul-Hassan, _et al._ , arXiv preprint arXiv:2110.01590 (2021).
  * Taylor  _et al._ [2008] J. M. Taylor, P. Cappellaro, L. Childress, L. Jiang, D. Budker, P. R. Hemmer, A. Yacoby, R. Walsworth, and M. D. Lukin, Nat. Phys. 4, 810 (2008).
  * Koenraad and Flatté [2011] P. M. Koenraad and M. E. Flatté, Nature Materials 10, 91 (2011).
  * Dolde  _et al._ [2011] F. Dolde _et al._ , Nature Physics 7, 459 (2011).
  * Dolde  _et al._ [2014] F. Dolde, M. W. Doherty, J. Michl, I. Jakobi, B. Naydenov, S. Pezzagna, J. Meijer, P. Neumann, F. Jelezko, N. B. Manson, and J. Wrachtrup, Phys. Rev. Lett. 112, 097603 (2014).
  * Schirhagl  _et al._ [2014] R. Schirhagl, K. Chang, M. Loretz, and C. L. Degen, Annual review of physical chemistry 65, 83 (2014).
  * van der Sar  _et al._ [2015] T. van der Sar, F. Casola, R. Walsworth, and A. Yacoby, Nat. Commun. 6, 7886 (2015).
  * Degen  _et al._ [2017] C. L. Degen, F. Reinhard, and P. Cappellaro, Rev. Mod. Phys. 89, 035002 (2017).
  * Flebus and Tserkovnyak [2018] B. Flebus and Y. Tserkovnyak, Phys. Rev. Lett. 121, 187204 (2018).
  * Casola  _et al._ [2018] F. Casola, T. van der Sar, and A. Yacoby, Nat. Rev. Mater. 3, 17088 (2018).
  * Mittiga  _et al._ [2018a] T. Mittiga, S. Hsieh, C. Zu, B. Kobrin, F. Machado, P. Bhattacharyya, N. Z. Rui, A. Jarmola, S. Choi, D. Budker, _et al._ , Phys. Rev. Lett. 121, 246402 (2018a).
  * Zhou  _et al._ [2020] B. B. Zhou, P. C. Jerger, K.-H. Lee, M. Fukami, F. Mujid, J. Park, and D. D. Awschalom, Phys. Rev. X 10, 011003 (2020).
  * Lee-Wong  _et al._ [2020] E. Lee-Wong, R. Xue, F. Ye, A. Kreisel, T. van der Sar, A. Yacoby, and C. R. Du, Nano Letters 20, 3284 (2020).
  * Rustagi  _et al._ [2020] A. Rustagi, I. Bertelli, T. van der Sar, and P. Upadhyaya, Phys. Rev. B 102, 220403(R) (2020).
  * Choi  _et al._ [2017] J. Choi, S. Choi, G. Kucsko, P. C. Maurer, B. J. Shields, H. Sumiya, S. Onoda, J. Isoya, E. Demler, F. Jelezko, N. Y. Yao, and M. D. Lukin, Phys. Rev. Lett. 118, 093601 (2017).
  * Kucsko  _et al._ [2018] G. Kucsko, S. Choi, J. Choi, P. C. Maurer, H. Zhou, R. Landig, H. Sumiya, S. Onoda, J. Isoya, F. Jelezko, E. Demler, N. Y. Yao, and M. D. Lukin, Phys. Rev. Lett. 121, 023601 (2018).
  * Davis  _et al._ [2021] E. J. Davis, B. Ye, F. Machado, S. A. Meynell, T. Mittiga, W. Schenken, M. Joos, B. Kobrin, Y. Lyu, D. Bluvstein, _et al._ , arXiv preprint arXiv:2103.12742 (2021).
  * Dwyer  _et al._ [2021] B. L. Dwyer, L. V. Rodgers, E. K. Urbach, D. Bluvstein, S. Sangtawesin, H. Zhou, Y. Nassab, M. Fitzpatrick, Z. Yuan, K. De Greve, _et al._ , arXiv preprint arXiv:2103.12757 (2021).
  * Jamonneau  _et al._ [2016] P. Jamonneau, M. Lesik, J. P. Tetienne, I. Alvizu, L. Mayer, A. Dréau, S. Kosen, J.-F. Roch, S. Pezzagna, J. Meijer, T. Teraji, Y. Kubo, P. Bertet, J. R. Maze, and V. Jacques, Phys. Rev. B 93, 024305 (2016).
  * Sangtawesin  _et al._ [2019] S. Sangtawesin, B. L. Dwyer, S. Srinivasan, J. J. Allred, L. V. H. Rodgers, K. De Greve, A. Stacey, N. Dontschuk, K. M. O’Donnell, D. Hu, D. A. Evans, C. Jaye, D. A. Fischer, M. L. Markham, D. J. Twitchen, H. Park, M. D. Lukin, and N. P. de Leon, Phys. Rev. X 9, 031052 (2019).
  * Shin  _et al._ [2013] C. S. Shin, C. E. Avalos, M. C. Butler, H.-J. Wang, S. J. Seltzer, R.-B. Liu, A. Pines, and V. S. Bajaj, Phys. Rev. B 88, 161412 (2013).
  * Kim  _et al._ [2015] M. Kim, H. J. Mamin, M. H. Sherwood, K. Ohno, D. D. Awschalom, and D. Rugar, Phys. Rev. Lett. 115, 087602 (2015).
  * Chrostoski  _et al._ [2018] P. Chrostoski, H. R. Sadeghpour, and D. H. Santamore, Phys. Rev. Applied 10, 064056 (2018).
  * Mittiga  _et al._ [2018b] T. Mittiga, S. Hsieh, C. Zu, B. Kobrin, F. Machado, P. Bhattacharyya, N. Z. Rui, A. Jarmola, S. Choi, D. Budker, and N. Y. Yao, Phys. Rev. Lett. 121, 246402 (2018b).
  * Balasubramanian  _et al._ [2009] G. Balasubramanian, P. Neumann, D. Twitchen, M. Markham, R. Kolesov, N. Mizuochi, J. Isoya, J. Achard, J. Beck, J. Tissler, V. Jacques, P. R. Hemmer, F. Jelezko, and J. Wrachtrup, Nature Materials 8, 383 (2009).
  * Meriles  _et al._ [2010] C. A. Meriles, L. Jiang, G. Goldstein, J. S. Hodges, J. Maze, M. D. Lukin, and P. Cappellaro, The Journal of Chemical Physics 133, 124105 (2010).
  * Rosskopf  _et al._ [2014] T. Rosskopf, A. Dussaux, K. Ohashi, M. Loretz, R. Schirhagl, H. Watanabe, S. Shikata, K. M. Itoh, and C. L. Degen, Phys. Rev. Lett. 112, 147602 (2014).
  * Myers  _et al._ [2014] B. A. Myers, A. Das, M. C. Dartiailh, K. Ohno, D. D. Awschalom, and A. C. Bleszynski Jayich, Phys. Rev. Lett. 113, 027602 (2014).
  * Romach  _et al._ [2015] Y. Romach, C. Müller, T. Unden, L. J. Rogers, T. Isoda, K. M. Itoh, M. Markham, A. Stacey, J. Meijer, S. Pezzagna, B. Naydenov, L. P. McGuinness, N. Bar-Gill, and F. Jelezko, Phys. Rev. Lett. 114, 017601 (2015).
  * Myers  _et al._ [2017] B. A. Myers, A. Ariyaratne, and A. C. B. Jayich, Phys. Rev. Lett. 118, 197201 (2017).
  * Anderson  _et al._ [2019] C. P. Anderson, A. Bourassa, K. C. Miao, G. Wolfowicz, P. J. Mintun, A. L. Crook, H. Abe, J. Ul Hassan, N. T. Son, T. Ohshima, and D. D. Awschalom, Science 366, 1225 (2019).
  * Candido and Flatté [2021] D. R. Candido and M. E. Flatté, PRX Quantum 2, 040310 (2021).
  * Trifunovic  _et al._ [2013] L. Trifunovic, F. L. Pedrocchi, and D. Loss, Phys. Rev. X 3, 041023 (2013).
  * Li  _et al._ [2015] P.-B. Li, Y.-C. Liu, S.-Y. Gao, Z.-L. Xiang, P. Rabl, Y.-F. Xiao, and F.-L. Li, Phys. Rev. Applied 4, 044003 (2015).
  * Li  _et al._ [2016] P.-B. Li, Z.-L. Xiang, P. Rabl, and F. Nori, Phys. Rev. Lett. 117, 015502 (2016).
  * Andrich  _et al._ [2017] P. Andrich, C. F. de las Casas, X. Liu, H. L. Bretscher, J. R. Berman, F. J. Heremans, P. F. Nealey, and D. D. Awschalom, npj Quantum Inf 3, 28 (2017).
  * Lemonde  _et al._ [2018] M.-A. Lemonde, S. Meesala, A. Sipahigil, M. J. A. Schuetz, M. D. Lukin, M. Loncar, and P. Rabl, Phys. Rev. Lett. 120, 213603 (2018).
  * Flebus and Tserkovnyak [2019] B. Flebus and Y. Tserkovnyak, Phys. Rev. B 99, 140403 (2019).
  * Mühlherr  _et al._ [2019] C. Mühlherr, V. O. Shkolnikov, and G. Burkard, Phys. Rev. B 99, 195413 (2019).
  * Zou  _et al._ [2020] J. Zou, S. K. Kim, and Y. Tserkovnyak, Phys. Rev. B 101, 014416 (2020).
  * Candido  _et al._ [2021] D. R. Candido, G. D. Fuchs, E. Johnston-Halperin, and M. E. Flatté, Mat. Quantum Technol. 1, 011001 (2021).
  * Neuman  _et al._ [2020] T. Neuman, D. S. Wang, and P. Narang, Phys. Rev. Lett. 125, 247702 (2020).
  * Wang  _et al._ [2021] D. S. Wang, T. Neuman, and P. Narang, The Journal of Physical Chemistry C 125, 6222 (2021).
  * Solanki  _et al._ [2020] A. B. Solanki, S. I. Bogdanov, A. Rustagi, N. R. Dilley, T. Shen, M. M. Rahman, W. Tong, P. Debashis, Z. Chen, J. Appenzeller, _et al._ , arXiv:2012.01497 (2020).
  * Fukami  _et al._ [2021] M. Fukami, D. R. Candido, D. D. Awschalom, and M. E. Flatté, PRX Quantum 2, 040314 (2021).
  * Debye and Bueche [1949] P. Debye and A. M. Bueche, Journal of Applied Physics 20, 518 (1949).
  * Debye  _et al._ [1957] P. Debye, H. R. Anderson, and H. Brumberger, Journal of Applied Physics 28, 679 (1957).
  * Sque  _et al._ [2006] S. J. Sque, R. Jones, and P. R. Briddon, Phys. Rev. B 73, 085313 (2006).
  * Sussmann [2009] R. S. Sussmann,  _CVD diamond for electronic devices and sensors_ , Vol. 26 (John Wiley & Sons, 2009).
  * Kaviani  _et al._ [2014] M. Kaviani, P. Deák, B. Aradi, T. Frauenheim, J.-P. Chou, and A. Gali, Nano Letters 14, 4772 (2014).
  * Chou and Gali [2017] J.-P. Chou and A. Gali, MRS Communications 7, 551–562 (2017).
  * Stacey  _et al._ [2019] A. Stacey, N. Dontschuk, J.-P. Chou, D. A. Broadway, A. K. Schenk, M. J. Sear, J.-P. Tetienne, A. Hoffman, S. Prawer, C. I. Pakes, A. Tadich, N. P. de Leon, A. Gali, and L. C. L. Hollenberg, Advanced Materials Interfaces 6, 1801449 (2019).
  * Morrison [1956] S. R. Morrison, Phys. Rev. 104, 619 (1956).
  * Turchette  _et al._ [2000] Q. A. Turchette, Kielpinski, B. E. King, D. Leibfried, D. M. Meekhof, C. J. Myatt, M. A. Rowe, C. A. Sackett, C. S. Wood, W. M. Itano, C. Monroe, and D. J. Wineland, Phys. Rev. A 61, 063418 (2000).
  * Kawarada [1996] H. Kawarada, Surface Science Reports 26, 205 (1996).
  * Maier  _et al._ [2000] F. Maier, M. Riedel, B. Mantel, J. Ristein, and L. Ley, Phys. Rev. Lett. 85, 3472 (2000).
  * Takeuchi  _et al._ [2003] D. Takeuchi, M. Riedel, J. Ristein, and L. Ley, Phys. Rev. B 68, 041304 (2003).
  * Crawford  _et al._ [2021] K. G. Crawford, I. Maini, D. A. Macdonald, and D. A. Moran, Progress in Surface Science 96, 100613 (2021).
  * Tamarat  _et al._ [2006] P. Tamarat, T. Gaebel, J. R. Rabeau, M. Khan, A. D. Greentree, H. Wilson, L. C. L. Hollenberg, S. Prawer, P. Hemmer, F. Jelezko, and J. Wrachtrup, Phys. Rev. Lett. 97, 083002 (2006).
  * de las Casas  _et al._ [2017] C. F. de las Casas, D. J. Christle, J. Ul Hassan, T. Ohshima, N. T. Son, and D. D. Awschalom, Applied Physics Letters 111, 262403 (2017).
  * Dubessy  _et al._ [2009] R. Dubessy, T. Coudreau, and L. Guidoni, Phys. Rev. A 80, 031402 (2009).
  * Low  _et al._ [2011] G. H. Low, P. F. Herskind, and I. L. Chuang, Phys. Rev. A 84, 053425 (2011).
  * [64] D. R. Candido and M. E. Flatté, in preparation .
  * Anderson [1954] P. W. Anderson, Journal of the Physical Society of Japan 9, 316 (1954).
  * Berthelot  _et al._ [2006] A. Berthelot, I. Favero, G. Cassabois, C. Voisin, C. Delalande, P. Roussignol, R. Ferreira, and J. M. Gérard, Nature Physics 2, 759 (2006).
  * Sinha  _et al._ [2014] S. K. Sinha, Z. Jiang, and L. B. Lurio, Advanced Materials 26, 7764 (2014).
  * Marshall  _et al._ [2021] M. C. Marshall, D. F. Phillips, M. J. Turner, M. J. H. Ku, T. Zhou, N. Delegan, F. J. Heremans, M. V. Holt, and R. L. Walsworth, Phys. Rev. Applied 16, 054032 (2021).
  * Note [1] As $\mathop{lim}\displaylimits_{\omega\rightarrow 0}S_{E}^{\delta\rho}(\omega)\rightarrow\infty$, this asymptotic behavior only holds for $\Omega\neq 0$. For $\Omega=0$ we still find asymptotic behaviors, but they are obtained by first performing the integral of Eq. (26) and then taking the limit $\Omega\rightarrow 0$.


proposition@alttheorem lemma@alttheorem corollary@alttheorem
question@alttheorem definition@alttheorem remark@alttheorem
construction@alttheorem

# On an infinite family of generalized PBIBDs

Daniel Kalmanovich  Einstein Institute of Mathematics
The Hebrew University of Jerusalem<EMAIL_ADDRESS>
(Date: August 13, 2015)

###### Abstract.

We consider a generalization of the notion of partially balanced incomplete
block designs (PBIBDs), by relaxing the requirement that the underlying
association scheme be commutative. An infinite family of such generalizations
is constructed, one for each prime power $q$ congruent to $1$ modulo $4$.

##  1\. Introduction and Preliminaries

Combinatorial designs are fascinating objects which have been studied for
centuries, and have origins in diverse areas (cf. the introduction of [5] for
a brief historical account). This paper focuses on a particular kind of
combinatorial design, originating in the theory of the design of experiments,
called partially balanced incomplete block design (PBIBD). Formally introduced
by Bose and Nair in [2], a PBIBD is a certain incidence structure together
with an underlying symmetric association scheme111There are longstanding
terminology disagreements in this area. In this paper association scheme means
homogeneous coherent configuration. which provides many regularity properties
to the whole structure. We refer to Bailey’s book [1] which gives a
statistician-friendly account of association schemes and their use in the
design of experiments.

Association schemes arose naturally in algebraic combinatorics, independent of
any applications to statistics, in works pertaining to permutation groups and,
independently, to the graph isomorphism problem. We refer to Cameron’s survey
[4] for a brief overview of the permutation groups aspect, and the celebrated
monograph [3] for the many graph theoretic aspects of association schemes.

Natural generalizations of the notion of a PBIBD were suggested and studied
over the years. Among these, notions that take a more general underlying
structure by relaxing the various requirements in the definition of an
association scheme were introduced and studied in [9] and in [10]. In [7],
Nair suggested a generalization of PBIBDs to incidence structures with
underlying association schemes which are not necessarily commutative.

In this paper we present a construction of an infinite family of such
generalizations of PBIBDs, one for each prime power $q$ congruent to $1$
modulo $4$. This construction stemmed from the analysis of designs constructed
by Nevo in [8].

We assume that the reader is familiar with basics of association scheme
theory, but for notational purposes we recall the definitions of an
association scheme and a $\mathrm{PBIBD}$.

######  Definition .

An association scheme with $m$ associate classes is a finite set $X$ together
with $m+1$ binary relations $R_{i}\subseteq X^{2}$, $0\leq i\leq m$, such
that:

  * •

$R_{0}=\left\\{(x,x)\ \middle|\ x\in X\right\\}$,

  * •

for each relation $R_{i}$, the inverse relation $R_{i}^{-1}=\left\\{(y,x)\
\middle|\ (x,y)\in R_{i}\right\\}$ is also an associate class,

  * •

$\bigcup\limits_{i=0}^{d}R_{i}=X^{2}$, and $R_{i}\cap R_{j}=\emptyset$
whenever $i\neq j$,

  * •

for any pair $(x,z)\in R_{k}$, the number of elements $y\in X$ such that
$(x,y)\in R_{i}$ and $(y,z)\in R_{j}$ depends only on $i,j,k$. This number is
denoted by $p_{ij}^{k}$.

The numbers in the last item above are called the intersection numbers of the
association scheme. We will also use the term $i$-associates when referring to
a pair $(x,y)\in R_{i}$.

A symmetric association scheme is one in which all the binary relations are
symmetric.

######  Definition .

A partially balanced incomplete block design with $m$ associate classes
denoted $\mathrm{PBIBD}(m)$ is a block design $D=(P,B)$ with $|P|=v$ points,
$|B|=b$ blocks, where each block has size $k$, each point lies on $r$ blocks,
and such that there is a symmetric association scheme defined on $P$ where, if
$(x,y)$ are $j$-th associates then they occur together in precisely
$\lambda_{j}$ blocks for $1\leq j\leq m$. The numbers
$v,b,r,k,\lambda_{j}(1\leq j\leq m)$ are called the parameters of the
$\mathrm{PBIBD}(m)$.

Our definition of a generalized PBIBD is verbatim section 1 dropping the word
symmetric. It follows that if $R_{j}$ is the inverse relation of $R_{i}$ then
$\lambda_{j}=\lambda_{i}$.

A special feature of our family of $\mathrm{GPBIBD}$s is that the number of
non-vanishing values of the $\lambda_{j}$s is constant: we have just the value
$5$ when $q=5^{\alpha}$, and the values $4$ and $1$ when $q$ is not a power of
$5$. Moreover, there is exactly one class for each non-zero value.

##  2\. The construction

Let $q=p^{\alpha}$ be a prime power congruent to $1$ modulo $4$, and consider
the vector space $V={{\mathbbm{F}}_{q}}^{2}$, where ${\mathbbm{F}}_{q}$
denotes the finite field with $q$ elements. Let $\omega$ be a generator of
${{\mathbbm{F}}_{q}}^{\times}$, and denote $i=\omega^{\frac{q-1}{4}}$. Our set
of points is
$P={\raisebox{1.99997pt}{$V\setminus\\{0\\}$}\left/\raisebox{-1.99997pt}{$\langle
i\rangle$}\right.}$, that is, we identify the four vectors $\pi$, $-\pi$,
$i\pi$ and $-i\pi$, for all non-zero vectors $\pi\in V$, to obtain the set $P$
of $\frac{q^{2}-1}{4}$ points. Let $G=\operatorname{PSL}(2,q)$, the quotient
group of the matrix group
$\operatorname{SL}(2,q)=\left\\{\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\ \middle|\ ad-bc=1\right\\}$ by its center
$Z(\operatorname{SL}(2,q))=\left\\{\scriptsize\begin{pmatrix}1&0\\\
0&1\end{pmatrix},\scriptsize\begin{pmatrix}-1&0\\\
0&-1\end{pmatrix}\right\\}$. Each element of $G$ is a set of $2$ matrices, and
consider the natural action of $G$ on $P$ via multiplication (of
representatives):

$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\cdot\scriptsize\begin{pmatrix}x\\\
y\end{pmatrix}=\scriptsize\begin{pmatrix}ax+by\\\ cx+dy\end{pmatrix},$

here $\scriptsize\begin{pmatrix}a&b\\\ c&d\end{pmatrix}$ is a representative
of the element $\left\\{\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix},\scriptsize\begin{pmatrix}-a&-b\\\
-c&-d\end{pmatrix}\right\\}\in G$, and $\scriptsize\begin{pmatrix}x\\\
y\end{pmatrix}$ is a representative of the element
$\left\\{\scriptsize\begin{pmatrix}x\\\
y\end{pmatrix},\scriptsize\begin{pmatrix}-x\\\
-y\end{pmatrix},\scriptsize\begin{pmatrix}ix\\\
iy\end{pmatrix},\scriptsize\begin{pmatrix}-ix\\\ -iy\end{pmatrix}\right\\}\in
P$. Define

$T=\left\\{\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix},\scriptsize\begin{pmatrix}0\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}1\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}1+i\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}i\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}1\\\ 1-i\end{pmatrix}\right\\},$

and denote $B=T^{G}$, the orbit of $T$ under the natural action of $G$ (the
induced action on $6$-sets). We have thus defined an incidence structure
$D=(P,B)$ on which $G$ acts transitively both on the points and on the blocks.
In fact, if we analyse the action of $G$ on pairs of points $P$ we can give a
geometric interpretation to the basic block $T$, and thus to all blocks of
$B$. Consider the points of $T$ as the vertices of a regular octahedron, with
edges as depicted in the following figure:

$\binom{0}{1}$$\binom{1}{i-1}$$\binom{i}{1}$$\binom{1}{1}$$\binom{1}{0}$$\binom{i+1}{1}$
Figure 1. The basic octahedron $T$

With this description at hand we will freely switch between the terms block
and octahedron when we refer to elements of $B$. Also, we will use the terms
edges and diagonals for pairs of points of $D$ which form an edge in an
octahedron, and pairs of points of $D$ which form a diagonal in an octahedron
(a pair of antipodal points), respectively. Of course, any other pair of
points do not lie in a block together.

######  Lemma .

  1. (1)

$G$ acts transitively on the edges of $D$.

  2. (2)

$G$ acts transitively on the diagonals of $D$.

Let us denote by ${\mathcal{S}}$ the Schurian association scheme of the action
of $G$ on $P$. Then we have

######  Corollary .

The block design $D=(P,B)$ is a $\mathrm{GPBIBD}$ with underlying association
scheme ${\mathcal{S}}$.

######  Lemma .

The stabilizer $G_{\pi}$ of a point $\pi\in P$ in $G$ is isomorphic to
${C_{p}}^{\alpha}\rtimes C_{2}$, where $q=p^{\alpha}$.

###### Proof.

Since $G$ acts transitively on $P$ it suffices to consider the stabilizer
$G_{\binom{1}{0}}$ of $\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}$. So let
$\scriptsize\begin{pmatrix}a&b\\\ c&d\end{pmatrix}\in G_{\binom{1}{0}}$, that
is

$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\cdot\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix}=\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}.$

It follows that $a=1$ or $i$, and $c=0$, thus

$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}1&b\\\
0&d\end{pmatrix}\quad\text{or}\quad\scriptsize\begin{pmatrix}i&b\\\
0&d\end{pmatrix}$

the determinant $=1$ condition implies that $d=1$ in the first case, and
$d=-i$ in the second case, so we have

$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}1&b\\\
0&1\end{pmatrix}\quad\text{or}\quad\scriptsize\begin{pmatrix}i&b\\\
0&-i\end{pmatrix}.$

Therefore

$G_{\binom{1}{0}}=\left\\{\scriptsize\begin{pmatrix}u&x\\\
0&u^{-1}\end{pmatrix}\ \middle|\
u\in\\{1,i\\},x\in{\mathbbm{F}}_{q}\right\\}.$

∎

The cases $q=5^{\alpha}$ and $q=p^{\alpha}$, $p\neq 5$ are different in an
essential way, so we analyse them separately. The reason for the difference
between the case of characteristic $5$ and other characteristics is the fact
that in the case of characteristic $5$, the subgroup $\\{1,-1,i,-i\\}$ forms
the multiplicative group of the prime subfield ${\mathbbm{F}}_{5}$ (and not
just a subgroup of the multiplicative group of the field). A manifestation of
this fact that will be useful for our purposes later is formulated in the
following lemma.

######  Lemma .

$1+i=-i$ if and only if $q=5^{\alpha}$.

###### Proof.

$1+i=-i\qquad\Longleftrightarrow\qquad 2i=-1\qquad\Longleftrightarrow\qquad
2i=i^{2}\qquad\Longleftrightarrow\qquad i=2$

so we have $-1=i^{2}=2^{2}=4$, that is, ${\mathbbm{F}}_{q}$ is a field of
characteristic $5$. ∎

## The case $p\neq 5$

In this case the stabilizer $G_{T}$ of the basic block $T$ in $G$ is a group
of rotations of the regular octahedron. It is in fact half of the full group
of rotations of the regular octahedron. There are two types of rotations of a
regular octahedron:

  1. (a)

rotation about the axis defined by a pair of antipodal points, and

  2. (b)

rotation about the axis formed by connecting the centers of a pair of opposite
faces.

In the full group of rotations of the octahedron, the rotations of type (a)
are of order $4$, and the rotations of type (b) are of order $3$ thus
producing a group of order $24$. In our group $G_{T}$ we take the square of
rotations of type (a) and all rotations of type (b), obtaining a group of
order $12$. In this description $G_{T}=\langle g,f\rangle$ where

  * •

$g=\scriptsize\begin{pmatrix}-i&-1+i\\\ 0&i\end{pmatrix}$ is the $180^{\circ}$
rotation around the axis defined by the antipodal points
$\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}$ and
$\scriptsize\begin{pmatrix}1\\\ 1-i\end{pmatrix}$, and

  * •

$f=\scriptsize\begin{pmatrix}i&1\\\ i&1-i\end{pmatrix}$ is the rotation around
the axis formed by connecting the centers of the opposite faces
$\left\\{\scriptsize\begin{pmatrix}1\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}1+i\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}\right\\}$ and
$\left\\{\scriptsize\begin{pmatrix}0\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}1\\\
1-i\end{pmatrix},\scriptsize\begin{pmatrix}i\\\ 1\end{pmatrix}\right\\}$.

######  Lemma .

The stabilizer $G_{\beta}$ of a block $\beta\in B$ in $G$ is isomorphic to the
alternating group $A_{4}$.

###### Proof.

First, since $G$ acts transitively on the block set $B$ it is enough to
consider $G_{T}$, the stabilizer of the basic block $T$. In fact, we will
explicitly write representatives for the elements of $G_{T}$. We will first
show that $|G_{T}|=12$, and then one can verify that

$\left\\{{\scriptsize\begin{aligned} &\begin{pmatrix}1&0\\\ 0&1\\\
\end{pmatrix},\begin{pmatrix}0&-i\\\ -i&-1\\\
\end{pmatrix},\begin{pmatrix}-1&i\\\ i&0\\\
\end{pmatrix},\begin{pmatrix}-i&0\\\ -1-i&i\\\
\end{pmatrix},\begin{pmatrix}-1+i&1\\\ i&-i\\\
\end{pmatrix},\begin{pmatrix}1&-1\\\ 1&0\\\ \end{pmatrix},\\\
&\begin{pmatrix}1&-1-i\\\ 1-i&-1\\\ \end{pmatrix},\begin{pmatrix}-1-i&i\\\
-1&i\\\ \end{pmatrix},\begin{pmatrix}i&1\\\ i&1-i\\\
\end{pmatrix},\begin{pmatrix}-i&-1+i\\\ 0&i\\\
\end{pmatrix},\begin{pmatrix}0&1\\\ -1&1\\\
\end{pmatrix},\begin{pmatrix}i&-i\\\ 1&-1-i\\\
\end{pmatrix}\end{aligned}}\right\\}$

is a set of $12$ representatives of distinct elements of $G_{T}$. We use the
orbit-stabilizer theorem:

$\left|G_{T}\right|=\left|{G_{T}}_{\binom{1}{0}}\right|\cdot\left|\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix}^{G_{T}}\right|.$

One can check directly that the elements in the above list leave $T$ fixed,
and that there are elements in that list that take
$\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}$ to each of the points in $T$,
so $\left|\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}^{G_{T}}\right|=6$. To
compute ${G_{T}}_{\binom{1}{0}}=G_{T}\cap G_{\binom{1}{0}}$ we recall that by
section 2

$G_{\binom{1}{0}}=\left\\{\scriptsize\begin{pmatrix}u&x\\\
0&u^{-1}\end{pmatrix}\ \middle|\
u\in\\{1,i\\},x\in{\mathbbm{F}}_{q}\right\\},$

now it is enough to show that the only elements in $G_{\binom{1}{0}}$ that fix
$T$ are $\scriptsize\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}$ and
$\scriptsize\begin{pmatrix}i&1-i\\\ 0&-i\end{pmatrix}$. So let
$\scriptsize\begin{pmatrix}u&x\\\ 0&u^{-1}\end{pmatrix}$ be an element that
stabilizes $T$, since $p\neq 5$ it follows from section 2 that
$\scriptsize\begin{pmatrix}u&x\\\
0&u^{-1}\end{pmatrix}\cdot\scriptsize\begin{pmatrix}0\\\
1\end{pmatrix}=\scriptsize\begin{pmatrix}x\\\ u^{-1}\end{pmatrix}$ is either
$\scriptsize\begin{pmatrix}0\\\ 1\end{pmatrix}$,
$\scriptsize\begin{pmatrix}i\\\ 1\end{pmatrix}$,
$\scriptsize\begin{pmatrix}1\\\ 1\end{pmatrix}$ or
$\scriptsize\begin{pmatrix}1+i\\\ 1\end{pmatrix}$. Each of the above four
cases gives $2$ options for $\scriptsize\begin{pmatrix}u&x\\\
0&u^{-1}\end{pmatrix}$, and out of these eight options the only elements that
indeed fix $T$ are $\scriptsize\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}$ and
$\scriptsize\begin{pmatrix}i&1-i\\\ 0&-i\end{pmatrix}$, so
$\left|{G_{T}}_{\binom{1}{0}}\right|=2$ and $\left|G_{T}\right|=2\cdot 6=12$.
The isomorphism $G_{T}\cong A_{4}$ can now be deduced by noticing that $G_{T}$
has no element of order $6$, and the only non-abelian group of order $12$
which has no elements of order $6$ is the alternating group $A_{4}$. ∎

######  Lemma .

  1. (a)

The group $G$ acts transitively on the edges of $D$.

  2. (b)

Every edge of $D$ is contained in $4$ ocathedra.

###### Proof.

  1. (a)

Since $G$ acts transitively on the octahedra, it suffices to show that $G_{T}$
acts transitively on the $12$ edges of $T$. This will follow from the orbit-
stabilizer theorem once we show that stabilizer ${G_{T}}_{e}$ of the edge
$e=\left\\{\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix},\scriptsize\begin{pmatrix}0\\\ 1\end{pmatrix}\right\\}$ (in
$G_{T}$) is trivial. So let $g\in{G_{T}}_{e}$, first we will show that $g$
must fix both vertices $\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}$ and
$\scriptsize\begin{pmatrix}0\\\ 1\end{pmatrix}$. Write
$g=\scriptsize\begin{pmatrix}a&b\\\ c&d\end{pmatrix}$, and assume that

$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\cdot\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix}=\scriptsize\begin{pmatrix}0\\\
1\end{pmatrix}\quad\text{and}\quad\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\cdot\scriptsize\begin{pmatrix}0\\\
1\end{pmatrix}=\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}$

then from the first equality it follows that $a=0$ and $c=1,i,-1$ or $-i$,
from the second equality if follows that $d=0$ and $b=1,i,-1$ or $-i$, and by
the determinant $=1$ condition we have

$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}0&-1\\\
1&0\end{pmatrix}\quad\text{or}\quad\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}0&i\\\ i&0\end{pmatrix}.$

Both elements above are not in $G_{T}$, thus $g$ must fix each of the vertices
of $e$. It follows that $g\in{G_{T}}_{\binom{1}{0}}$, thus
$g=\scriptsize\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}$ or
$g=\scriptsize\begin{pmatrix}-i&-1+i\\\ 0&i\end{pmatrix}$, but
$\scriptsize\begin{pmatrix}-i&-1+i\\\ 0&i\end{pmatrix}$ does not fix
$\scriptsize\begin{pmatrix}0\\\ 1\end{pmatrix}$ so we must have
$g=\scriptsize\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}$ and so ${G_{T}}_{e}$ is
trivial.

  2. (b)

From the calculation in the proof of $(a)$ it follows that

$G_{e}=\left\langle\scriptsize\begin{pmatrix}0&-1\\\
1&0\end{pmatrix},\scriptsize\begin{pmatrix}0&i\\\
i&0\end{pmatrix}\right\rangle$

thus $\left|G_{e}\right|=4$. Using the orbit-stabilizer theorem we can now
compute the total number of edges of $D$, that is, the size of the orbit of
$e$ under the action of $G$:

$\left|e^{G}\right|=\frac{\left|G\right|}{\left|G_{e}\right|}=\frac{\frac{q(q^{2}-1)}{2}}{4}=\frac{q(q^{2}-1)}{8}.$

Let $E$ denote the set of edges of $D$. Now we apply a standard double-
counting argument, we count the size of the set

$S=\left\\{(\epsilon,\beta)\ \middle|\ \epsilon\in E,\beta\in B,\epsilon\text{
is an edge of }\beta\right\\}$

in two ways:

     * •

$\left|S\right|=12\cdot\left|B\right|$ because every block has $12$ edges, and

     * •

$\left|S\right|=\left|E\right|\cdot x$, where $x$ is the number of blocks that
contain a given edge.

Therefore, we have

$x=\frac{12\cdot\left|B\right|}{\left|E\right|}=\frac{12\cdot\frac{q(q^{2}-1)}{24}}{\frac{q(q^{2}-1)}{8}}=4.$

∎

######  Lemma .

  1. (a)

The group $G$ acts transitively on the diagonals of $D$.

  2. (b)

Every diagonal of $D$ is contained in $1$ ocathedron.

###### Proof.

  1. (a)

As in the previous proofs, since $G$ acts transitively on the octahedra, it
suffices to show that $G_{T}$ acts transitively on the $3$ diagonals of $T$.
To see that, it is enough to consider just the rotation $f$:

$\left\\{\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix},\scriptsize\begin{pmatrix}1\\\
1-i\end{pmatrix}\right\\}\overset{f}{\mapsto}\left\\{\scriptsize\begin{pmatrix}1\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}i\\\
1\end{pmatrix}\right\\}\overset{f}{\mapsto}\left\\{\scriptsize\begin{pmatrix}1+i\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}0\\\ 1\end{pmatrix}\right\\}.$

  2. (b)

The following arguments are similar to those given in the proof of The case
$p\neq 5$, we will compute the stabilizer $G_{d}$ of the diagonal
$d=\left\\{\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix},\scriptsize\begin{pmatrix}1\\\ 1-i\end{pmatrix}\right\\}$ in
order to obtain the total number of diagonals in $D$ via the orbit-stabilizer
theorem, and then use a double-counting argument to compute the number of
octahedra that contain a given diagonal. So let
$g=\scriptsize\begin{pmatrix}a&b\\\ c&d\end{pmatrix}\in G_{d}$, then we have
either

$(I)\quad\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\cdot\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix}=\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix}\quad\text{and}\quad\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\cdot\scriptsize\begin{pmatrix}1\\\
1-i\end{pmatrix}=\scriptsize\begin{pmatrix}1\\\ 1-i\end{pmatrix}$

or

$(II)\quad\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\cdot\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix}=\scriptsize\begin{pmatrix}1\\\
1-i\end{pmatrix}\quad\text{and}\quad\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}\cdot\scriptsize\begin{pmatrix}1\\\
1-i\end{pmatrix}=\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}.$

In case $(I)$, the first equality implies $a=1$ or $i$, and $c=0$, thus
$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}1&b\\\ 0&1\end{pmatrix}$ or
$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}i&b\\\ 0&-i\end{pmatrix}$, plugging
this into the second equation we get $b=0$ in the first case and $b=1-i$ in
the second case, so

$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}1&0\\\ 0&1\end{pmatrix}\text{ or
}\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}i&1-i\\\ 0&-i\end{pmatrix}.$

In case $(II)$, the first equality implies $a=1$ or $i$, and $c=1-i$ or $1+i$
(respectively), thus $\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}1&b\\\ 1-i&d\end{pmatrix}$ or
$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}i&b\\\ 1+i&d\end{pmatrix}$,
plugging this into the second equation we get $d=-1$ in the first case, and
$d=-i$ in the second, so $\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}1&b\\\ 1-i&-1\end{pmatrix}$ or
$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}i&b\\\ 1+i&-i\end{pmatrix}$, now
the determinant $=1$ condition determines $b$ and we get

$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}1&-1-i\\\ 1-i&-1\end{pmatrix}\text{
or }\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}=\scriptsize\begin{pmatrix}i&0\\\ 1+i&-i\end{pmatrix}.$

It is now easy to check that the four representative matrices we found form a
group, thus $\left|G_{d}\right|=4$. Again, the orbit-stabilizer theorem allows
us to compute the total number of diagonals of $D$, that is, the size of the
orbit of $d$ under the action of $G$:

$\left|d^{G}\right|=\frac{\left|G\right|}{\left|G_{d}\right|}=\frac{\frac{q(q^{2}-1)}{2}}{4}=\frac{q(q^{2}-1)}{8}.$

Let $\Delta$ denote the set of diagonals of $D$. Yet again we apply a double-
counting argument, we count the size of the set

$S=\left\\{(\delta,\beta)\ \middle|\ \delta\in\Delta,\beta\in B,\delta\text{
is a diagonal of }\beta\right\\}$

in two ways:

     * •

$\left|S\right|=3\cdot\left|B\right|$ because every block has $3$ diagonals,
and

     * •

$\left|S\right|=\left|\Delta\right|\cdot x$, where $x$ is the number of blocks
that contain a given diagonal.

Therefore, we have

$x=\frac{3\cdot\left|B\right|}{\left|\Delta\right|}=\frac{3\cdot\frac{q(q^{2}-1)}{24}}{\frac{q(q^{2}-1)}{8}}=1.$

∎

######  Theorem 1.

The block design $D$ is a $\mathrm{GPBIBD}(m)$ with parameters

$v=\frac{q^{2}-1}{4},\quad b=\frac{q(q^{2}-1)}{24},\quad k=6,\quad
r=q,\quad\lambda_{1}=4,\quad\lambda_{2}=1,\quad\lambda_{3}=\dots=\lambda_{m}=0.$

The underlying association scheme is the Schurian scheme ${\mathcal{S}}(G,P)$
of $(G,P)$, and the number of associate classes is

$m=\frac{q-3}{2}.$

###### Proof.

The number of points is clear from the construction. The number of blocks can
be calculated by using the orbit-stabilizer theorem and the fact that, by The
case $p\neq 5$, $\left|G_{T}\right|=12$:

$b=\left|B\right|=\left|T^{G}\right|=\frac{\left|G\right|}{\left|G_{T}\right|}=\frac{\frac{q(q^{2}-1)}{2}}{12}=\frac{q(q^{2}-1)}{24}.$

The fact that each block has size $k=6$ follows from the fact that the basic
block $T$ has size $6$, and that the set of blocks is the orbit $T^{G}$ of $T$
under the action of $G$. To prove that each point lies on $r=q$ blocks it
suffices to show that $\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}$ lies on
$q$ blocks, because $G$ acts transitively on $P$. For that purpose we will
show that $\left|T^{G_{\binom{1}{0}}}\right|=q$. Again, the orbit-stabilizer
theorem applied to $G_{\binom{1}{0}}$ will do the job:

$2q=\left|G_{\binom{1}{0}}\right|=\left|{G_{\binom{1}{0}}}_{T}\right|\cdot\left|T^{G_{\binom{1}{0}}}\right|=2\cdot\left|T^{G_{\binom{1}{0}}}\right|,$

thus $\left|T^{G_{\binom{1}{0}}}\right|=q$. A pair of points of $P$ is either
an edge of $D$, a diagonal of $D$ or neither (that is a pair that does not lie
in a block together), it follows from The case $p\neq 5$ that an edge lies on
$\lambda_{1}=4$ blocks, from The case $p\neq 5$ that a diagonal lies on
$\lambda_{2}=1$ block, and any other pair lies on $\lambda_{j}=0$ ($3\leq
j\leq m$) blocks. Next, the fact that the Schurian scheme of the action
$(G,P)$ is the underlying association scheme follows from the fact that $G$
acts transitively on the edges of $D$ (proved in The case $p\neq 5$), and on
the diagonals of $D$ (proved in The case $p\neq 5$).

Finally, we compute the number of associate classes, that is, the number of
classes of ${\mathcal{S}}(G,P)$, or in other words, the rank of the transitive
permutation group $(G,P)$ minus $1$. The rank of a transitive permutation
group is equal to number of orbits of the stabilizer of a point, so we compute
the number of orbits of $\left(G_{\binom{1}{0}},P\right)$. Recall that

$G_{\binom{1}{0}}=\left\\{\scriptsize\begin{pmatrix}u&x\\\
0&u^{-1}\end{pmatrix}\ \middle|\
u\in\\{1,i\\},x\in{\mathbbm{F}}_{q}\right\\},$

so for each point of the form $\scriptsize\begin{pmatrix}a\\\
0\end{pmatrix}\in P$ we have $\scriptsize\begin{pmatrix}a\\\
0\end{pmatrix}^{G_{\binom{1}{0}}}=\left\\{\scriptsize\begin{pmatrix}a\\\
0\end{pmatrix}\right\\}$, these are $\frac{q-1}{4}$ orbits. For points of the
form $\scriptsize\begin{pmatrix}a\\\ b\end{pmatrix}\in P$, with $b\neq 0$, we
have

$\scriptsize\begin{pmatrix}a\\\
b\end{pmatrix}^{G_{\binom{1}{0}}}=\left\\{\scriptsize\begin{pmatrix}y\\\
b\end{pmatrix}\ \middle|\ y\in{\mathbbm{F}}_{q}\right\\},$

because $\scriptsize\begin{pmatrix}1&\frac{y-a}{b}\\\ 0&1\end{pmatrix}$ takes
$\scriptsize\begin{pmatrix}a\\\ b\end{pmatrix}$ to
$\scriptsize\begin{pmatrix}y\\\ b\end{pmatrix}$, so we get $\frac{q-1}{4}$
more orbits. Altogether we counted the $\frac{q-1}{2}$ orbits of
$G_{\binom{1}{0}}$, thus the number of associate classes is $m=\frac{q-3}{2}$.
∎

## The case $p=5$

Let $q=5^{\alpha}$. In this case the basic block $T$ is in fact just the
projective line over ${\mathbbm{F}}_{5}$, and we have

######  Lemma .

The stabilizer $G_{\beta}$ of a block $\beta\in B$ in $G$ is isomorphic to
$\operatorname{PSL}(2,5)$.

###### Proof.

As before, by the transitive action of $G$ on $B$ it is enough to consider
$G_{T}$. By section 2, since $q=5^{\alpha}$, we have that $1+i=-i$. It follows
that $1-i=-1$ and thus we have

$T=\left\\{\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix},\scriptsize\begin{pmatrix}0\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}1\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}-i\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}i\\\
1\end{pmatrix},\scriptsize\begin{pmatrix}-1\\\ 1\end{pmatrix}\right\\},$

that is, $T$ is the projective line over ${\mathbbm{F}}_{5}$, and it follows
that subgroup of $\operatorname{PSL}(2,q)$ that fixes the projective line over
${\mathbbm{F}}_{5}$ is $\operatorname{PSL}(2,5)$. ∎

######  Corollary .

$G$ acts transitively on pairs of points that lie together in a block.

We will refer to pairs of points that lie together on a block adjacent points,
and correspondingly to pairs of points that do not lie together on a block as
non-adjacent points.

######  Theorem 2.

The block design $D$ is a $\mathrm{GPBIBD}(m)$ with parameters

$v=\frac{q^{2}-1}{4},\quad b=\frac{q(q^{2}-1)}{120},\quad k=6,\quad
r=\frac{q}{5}=5^{\alpha-1},\quad\lambda_{1}=1,\quad\lambda_{2}=\dots=\lambda_{m}=0.$

The underlying association scheme is the Schurian scheme ${\mathcal{S}}(G,P)$
of $(G,P)$, and the number of associate classes is

$m=\frac{q-3}{2}.$

###### Proof.

The number of points is clear from the construction. The number of blocks can
be calculated by using the orbit-stabilizer theorem and the fact that, by The
case $p=5$, $\left|G_{T}\right|=60$:

$b=\left|B\right|=\left|T^{G}\right|=\frac{\left|G\right|}{\left|G_{T}\right|}=\frac{\frac{q(q^{2}-1)}{2}}{60}=\frac{q(q^{2}-1)}{120}.$

The fact that each block has size $k=6$ follows from the fact that the basic
block $T$ has size $6$, and that the set of blocks is the orbit $T^{G}$ of $T$
under the action of $G$. To prove that each point lies on
$r=\frac{q}{5}=5^{\alpha-1}$ blocks it suffices to show that
$\scriptsize\begin{pmatrix}1\\\ 0\end{pmatrix}$ lies on $5^{\alpha-1}$ blocks,
because $G$ acts transitively on $P$. For that purpose we will show that
$\left|T^{G_{\binom{1}{0}}}\right|=5^{\alpha-1}$. Again, the orbit-stabilizer
theorem applied to $G_{\binom{1}{0}}$ will do the job:

$2q=\left|G_{\binom{1}{0}}\right|=\left|{G_{\binom{1}{0}}}_{T}\right|\cdot\left|T^{G_{\binom{1}{0}}}\right|.$

By The case $p=5$ we have

${G_{\binom{1}{0}}}_{T}=G_{T}\cap G_{\binom{1}{0}}=\operatorname{PSL}(2,5)\cap
G_{\binom{1}{0}},$

that is, ${G_{\binom{1}{0}}}_{T}=\left\\{\scriptsize\begin{pmatrix}u&x\\\
0&u^{-1}\end{pmatrix}\ \middle|\ u\in\\{1,i\\},x\in{\mathbbm{F}}_{5}\right\\}$
and $\left|{G_{\binom{1}{0}}}_{T}\right|=10$, thus
$2q=10\cdot\left|T^{G_{\binom{1}{0}}}\right|$ and therefore
$\left|T^{G_{\binom{1}{0}}}\right|=\frac{2q}{10}=\frac{2\cdot
5^{\alpha}}{2\cdot 5}=5^{\alpha-1}$. The number of pairs of adjacent points is

$\binom{6}{2}\cdot\left|B\right|=15\cdot\frac{q(q^{2}-1)}{120}=\frac{q(q^{2}-1)}{8}.$

To show that each pair of adjacent points determines a unique block it is
enough to show that the stabilizer of a pair of adjacent points stabilizes
that whole block. So by The case $p=5$ we may assume that
$e=\left\\{\scriptsize\begin{pmatrix}1\\\
0\end{pmatrix},\scriptsize\begin{pmatrix}0\\\ 1\end{pmatrix}\right\\}$ is our
pair of adjacent points on the block $T$. We need to show that $G_{e}\leq
G_{T}$. Recall that the calculation in The case $p\neq 5$ showed that

$G_{e}=\left\langle\scriptsize\begin{pmatrix}0&-1\\\
1&0\end{pmatrix},\scriptsize\begin{pmatrix}0&i\\\
i&0\end{pmatrix}\right\rangle,$

in this case, clearly $G_{e}\leq\operatorname{PSL}(2,5)=G_{T}$. Finally, the
fact that the Schurian scheme of the action $(G,P)$ is the underlying
association scheme follows from the fact that $G$ acts transitively on the
pairs of adjacent points $D$ (stated in The case $p=5$), the computation of
the number of associate classes is identical to that given in Theorem 1. ∎

At this point it seems interesting to ask:

######  Question .

Can we find an underlying association scheme for $D$ with less vanishing
classes?

The following two sections deal with this question. In section 3 we find the
largest (in a certain natural sense) group acting on $D$ producing an
underlying Schurian association scheme with less vanishing classes. In section
4 we compute the underlying association scheme with the smallest possible
number of vanishing classes for all $\mathrm{GPBIBD}$s with $q\leq 169$, and
propose a general conjecture about the minimal association schemes.

##  3\. A smaller rank underlying association scheme

In this section we prove that in fact we have a larger group acting on $D$.
Let $\phi$ be the Frobenius automorphism of ${\mathbbm{F}}_{q}$, that is, the
map that takes each element $a\in{\mathbbm{F}}_{q}$ to its $p$-th power
$a^{p}\in{\mathbbm{F}}_{q}$. We denote by $F=\langle\phi\rangle$ the
automorphism group $\operatorname{Aut}({\mathbbm{F}}_{q})$. The projective
special semilinear group, commonly denoted by $\operatorname{P\Sigma L}(2,q)$,
is the semidirect product $\operatorname{PSL}(2,q)\rtimes F$, where the cyclic
group $F$ is acting on $\operatorname{PSL}(2,q)$ in the natural way:

$\scriptsize\begin{pmatrix}a&b\\\
c&d\end{pmatrix}^{\phi}=\scriptsize\begin{pmatrix}a^{p}&b^{p}\\\
c^{p}&d^{p}\end{pmatrix}.$

######  Theorem 3.

The largest automorphism group ${\mathcal{G}}$ of $D$ that preserves the
octahedral structure is isomorphic to $\operatorname{P\Sigma L}(2,q)\times
C_{2}$.

###### Proof.

Clearly $\operatorname{P\Sigma L}(2,q)$ acts on $D$. The additional involution
$C_{2}=\langle\sigma\rangle$ is the missing $90^{\circ}$ rotation about the
axis formed by two antipodal vertices of the octahedron. ∎

### Computing the number of orbits of $F$ on $P$

We use Burnside’s lemma and standard Galois cohomology tools to compute the
number of orbits
${\raisebox{1.99997pt}{$P$}\left/\raisebox{-1.99997pt}{$F$}\right.}$ of $F$ on
$P={\raisebox{1.99997pt}{${{\mathbbm{F}}_{q}}^{\times}$}\left/\raisebox{-1.99997pt}{$\mu_{4}$}\right.}$
and $q=p^{n}$. As usual we denote by $P^{g}$ the fixed points of $P$ under the
action of $g$. Burnside’s lemma establishes that

$\left|{\raisebox{1.99997pt}{$P$}\left/\raisebox{-1.99997pt}{$F$}\right.}\right|=\frac{1}{|F|}\sum_{g\in
F}\left|P^{g}\right|.$

For each divisor $d$ of $n$ we have a subfield ${\mathbbm{F}}_{p^{d}}$ of
${\mathbbm{F}}_{p^{n}}$, the elements of ${\mathbbm{F}}_{p^{d}}$ are precisely
the fixed points of $\phi^{d}$. Moreover, if $k$ is a multiple of $d$ that
does not divide $n$ then the fixed points of $\phi^{k}$ are again the elements
of ${\mathbbm{F}}_{p^{d}}$, thus we can rewrite the above sum as follows:

$\left|{\raisebox{1.99997pt}{$P$}\left/\raisebox{-1.99997pt}{$F$}\right.}\right|=\frac{1}{|F|}\sum_{g\in
F}\left|P^{g}\right|=\frac{1}{|F|}\sum_{d|n}\varphi\left(\frac{n}{d}\right)\left|P^{\phi^{d}}\right|$

where $\varphi$ is Euler totient function. To compute
$\left|P^{\phi^{d}}\right|$ we use Galois cohomology. For clearer notation we
denote $K:=\phi^{d}$. The short exact sequence

$1\longrightarrow\mu_{4}\longrightarrow{{\mathbbm{F}}_{q}}^{\times}\longrightarrow{\raisebox{1.99997pt}{${{\mathbbm{F}}_{q}}^{\times}$}\left/\raisebox{-1.99997pt}{$\mu_{4}$}\right.}\longrightarrow
1$

provides a long exact sequence in group cohomology

$1\longrightarrow{\mu_{4}}^{K}\longrightarrow\left({{\mathbbm{F}}_{q}}^{\times}\right)^{K}\longrightarrow\left({\raisebox{1.99997pt}{${{\mathbbm{F}}_{q}}^{\times}$}\left/\raisebox{-1.99997pt}{$\mu_{4}$}\right.}\right)^{K}\longrightarrow
H^{1}(K,\mu_{4})\longrightarrow
H^{1}(K,{{\mathbbm{F}}_{q}}^{\times})\longrightarrow\dots$

By Hilbert’s Theorem 90 we have

$H^{1}(K,{{\mathbbm{F}}_{q}}^{\times})=1$

so we get the following exact sequence

$1\longrightarrow{\mu_{4}}^{K}\longrightarrow\left({{\mathbbm{F}}_{q}}^{\times}\right)^{K}\longrightarrow\left({\raisebox{1.99997pt}{${{\mathbbm{F}}_{q}}^{\times}$}\left/\raisebox{-1.99997pt}{$\mu_{4}$}\right.}\right)^{K}\longrightarrow
H^{1}(K,\mu_{4})\longrightarrow 1$

and we are interested in

$\left|\left({\raisebox{1.99997pt}{${{\mathbbm{F}}_{q}}^{\times}$}\left/\raisebox{-1.99997pt}{$\mu_{4}$}\right.}\right)^{K}\right|=\frac{\left|H^{1}(K,\mu_{4})\right|\cdot\left|\left({{\mathbbm{F}}_{q}}^{\times}\right)^{K}\right|}{\left|{\mu_{4}}^{K}\right|}.$

Clearly, we have
$\left|\left({{\mathbbm{F}}_{q}}^{\times}\right)^{K}\right|=p^{d}-1$. We deal
with cases $p\equiv 1\pmod{4}$ and $p\equiv 3\pmod{4}$ separately.

### The case $p\equiv 1\pmod{4}$

In this case $\mu_{4}$ is already contained in ${\mathbbm{F}}_{p}$ and thus is
fixed by $K$ so $\left|{\mu_{4}}^{K}\right|=4$ and we have

$\left|\left({\raisebox{1.99997pt}{${{\mathbbm{F}}_{q}}^{\times}$}\left/\raisebox{-1.99997pt}{$\mu_{4}$}\right.}\right)^{K}\right|=\left|H^{1}(K,\mu_{4})\right|\cdot\frac{p^{d}-1}{4}.$

Now, since $\mu_{4}$ is a trivial $K$-module we have

$H^{1}(K,\mu_{4})=\operatorname{Hom}(K,\mu_{4})\cong\begin{cases}C_{4}&|K|\equiv
0\pmod{4}\\\ C_{2}&|K|\equiv 2\pmod{4}\\\ 1&\text{otherwise}\end{cases}$

so if we denote $h(d):=|H^{1}(K,\mu_{4})|$ (this number depends on the residue
of $\frac{n}{d}$ modulo $4$) we have

$\left|{\raisebox{1.99997pt}{$P$}\left/\raisebox{-1.99997pt}{$F$}\right.}\right|=\frac{1}{\left|F\right|}\sum_{d|n}\varphi\left(\frac{n}{d}\right)\cdot
h(d)\cdot\frac{p^{d}-1}{4}.$

### The case $p\equiv 3\pmod{4}$

In this case $\mu_{4}$ is contained in ${\mathbbm{F}}_{p^{2}}$ but not in
${\mathbbm{F}}_{P}$ so for every odd power of $\phi$ there are just $2$ fixed
points, and for every even power of $\phi$ there are $4$ fixed points, that
is:

$\left|{\mu_{4}}^{K}\right|=\begin{cases}2&d\text{ odd}\\\ 4&d\text{
even}\end{cases}.$

The computation of $\left|H^{1}(K,\mu_{4})\right|$ in this case is a bit more
involved. If $\frac{n}{d}$ is even then $\mu_{4}$ is again a trivial
$K$-module and we have

$H^{1}(K,\mu_{4})=\operatorname{Hom}(K,\mu_{4})\cong\begin{cases}C_{4}&|K|\equiv
0\pmod{4}\\\ C_{2}&|K|\equiv 2\pmod{4}\\\ \end{cases}.$

For $d$ odd, $\mu_{4}$ is no longer a trivial $K$-module. There are two cyclic
groups here, and an action of one of them on the other, considered as a
module, so we will use additive notation for $\mu_{4}$ (that is, the group of
integers modulo $4$), and write $a\cdot m$ to denote the action of $a$ on $m$,
where $a\in K$ and $m\in\mu_{4}$, the operation in $K$ will be denoted simply
$ab$ for $a,b\in K$. Fix a generator $\alpha$ of $K$, the the action of $K$ on
$\mu_{4}$ is

$\alpha\cdot 0=0,\quad\alpha\cdot 1=3,\quad\alpha\cdot 2=2,\quad\alpha\cdot
3=1.$

The $1$-cocycles $Z^{1}(K,\mu_{4})$ are thus the functions
$f:K\longrightarrow\mu_{4}$ satisfying

$f(ab)=f(a)+a\cdot f(b)$

and the $1$-coboundaries $B^{1}(K,\mu_{4})$ are those satisfying

$f(a)=a\cdot m-m\qquad\text{for some}\quad m\in\mu_{4}.$

A $1$-cocycle is determined by its value on $\alpha$. So let $f\in
B^{1}(K,\mu_{4})$, one can verify that out of the $4$ options (one option for
each $m\in\mu_{4}$) there are just $2$ elements in $B^{1}(K,\mu_{4})$,
explicitly, these are the functions

$f(\alpha)=\alpha\cdot 0-0\quad\text{and}\quad f(\alpha)=\alpha\cdot 1-1.$

Consider now $Z^{1}(K,\mu_{4})$. Let $f\in Z^{1}(K,\mu_{4})$ then out of the
$4$ options

$f(\alpha)=0,\quad f(\alpha)=1,\quad f(\alpha)=2,\quad f(\alpha)=3,$

it is easy to verify that $f(\alpha)=0$ and $f(\alpha)=2$ are both in
$B^{1}(K,\mu_{4})$. Also, $f(\alpha)=1$ and $f(\alpha)=3$ are cohomologous
(their difference is in $B^{1}(K,\mu_{4})$), so altogether we get

$H^{1}(K,\mu_{4})={\raisebox{1.99997pt}{$Z^{1}(K,\mu_{4})$}\left/\raisebox{-1.99997pt}{$B^{1}(K,\mu_{4})$}\right.}\cong
C_{2}.$

Denote
$h(d):=\frac{\left|H^{1}(K,\mu_{4})\right|}{\left|{\mu_{4}}^{K}\right|}$
summarizing the above we have

$h(d)=\begin{cases}1&\text{$d$ odd; or $d$ even and $\frac{n}{d}\equiv
0\pmod{4}$}\\\ \frac{1}{2}&\text{$d$ even and $\frac{n}{d}\equiv
2\pmod{4}$}\\\ \frac{1}{4}&\text{$d$ even and $\frac{n}{d}$ odd}\end{cases}$

and

$\left|{\raisebox{1.99997pt}{$P$}\left/\raisebox{-1.99997pt}{$F$}\right.}\right|=\frac{1}{\left|F\right|}\sum_{d|n}\varphi\left(\frac{n}{d}\right)\cdot
h(d)\cdot(p^{d}-1).$

######  Corollary .

$D$ is a $\mathrm{GPBIBD}(m)$ with
$m=2\cdot\left|{\raisebox{1.99997pt}{$P$}\left/\raisebox{-1.99997pt}{$F$}\right.}\right|-1$
associate classes.

##  4\. The smallest rank underlying association scheme

We have constructed all the $\mathrm{GPBIBD}$s with $q\leq 169$ on GAP [6], we
used the program stabil2 to perform WL-stabilization on the partition of
$P\times P$ into $\lambda$-classes (a $\lambda$-class consists of all pairs in
$P\times P$ that lie on $\lambda$ blocks) in order to obtain an underlying
association scheme with the smallest number of classes possible for each
incidence structure. It turns out that in a considerable number of cases we in
fact obtain an underlying scheme of smaller rank than what is provided by
section 3. It follows from Theorem 3 that in those cases the underlying scheme
is non-Schurian.

$q$ |  Class Number in section 3 |  Smallest Class Number |  Parameters $(v,b,r,k;\lambda_{1},\dots,\lambda_{m})$ |  Remarks
---|---|---|---|---
$9$ | $3$ | $3$ | $(20,30,6,9;4,1,0)$ |
$13$ | $5$ | $5$ | $(42,91,6,13;4,1,0,0,0)$ |
$17$ | $7$ | $7$ | $(72,204,6,17;4,1,0,\dots,0)$ |
$25$ | $7$ | $3$ | $(156,130,6,5;1,0,0)$ |  Non-Schurian. This is an antipodal distance regular graph of diameter $3$, a $6$-fold cover of $K_{4}$.
$29$ | $13$ | $13$ | $(210,1015,6,29;4,1,0,\dots,0)$ |
$37$ | $17$ | $17$ | $(342,2109,6,37;4,1,0,\dots,0)$ |
$41$ | $19$ | $11$ | $(420,2870,6,41;4,1,0,\dots,0)$ |  Non-Schurian. Non-commutative.
$49$ | $17$ | $13$ | $(600,4900,6,49;4,1,0,\dots,0)$ |  Non-Schurian. Non-commutative.
$53$ | $25$ | $25$ | $(702,6201,6,53;4,1,0,\dots,0)$ |
$61$ | $29$ | $29$ | $(930,9455,6,61;4,1,0,\dots,0)$ |
$73$ | $35$ | $35$ | $(1332,16206,6,73;4,1,0,\dots,0)$ |
$81$ | $13$ | $5$ | $(1640,22140,6,81;4,1,0,0,0)$ |  non-Schurian.
$89$ | $43$ | $43$ | $(1980,29370,6,89;4,1,0,\dots,0)$ |
$97$ | $47$ | $47$ | $(2352,38024,6,97;4,1,0,\dots,0)$ |
$101$ | $49$ | $49$ | $(2550,42925,6,101;4,1,0,\dots,0)$ |
$109$ | $53$ | $19$ | $(2970,53955,6,109;4,1,0,\dots,0)$ |  Non-Schurian. Non-commutative.
$113$ | $55$ | $15$ | $(3192,60116,6,113;4,1,0,\dots,0)$ |  Non-Schurian. Non-commutative.
$121$ | $39$ | $21$ | $(3660,73810,6,121;4,1,0,\dots,0)$ |  non-Schurian.
$125$ |  | $3$ | $(3906,16275,6,25;1,0,\dots,0)$ |  non-Schurian. This is an antipodal distance regular graph of diameter $3$, a $31$-fold cover of $K_{4}$.
$137$ |  |  |  |
$149$ |  |  |  |
$157$ |  |  |  |
$169$ | $47$ | $7$ | $(7140,201110,6,169;4,1,0,\dots,0)$ |  non-Schurian.

### Acknowledgements

I thank Rosemary Bailey and Eran Nevo for helpful comments on earlier versions
of this text.

## References

  * [1] R. A. Bailey.  Association schemes, volume 84 of Cambridge Studies in Advanced Mathematics.  Cambridge University Press, Cambridge, 2004.  Designed experiments, algebra and combinatorics.
  * [2] R. C. Bose and K. R. Nair.  Partially balanced incomplete block designs.  Sankhyā: The Indian Journal of Statistics (1933-1960), 4(3):337–372, 1939.
  * [3] A. E. Brouwer, A. M. Cohen, and A. Neumaier.  Distance-regular graphs, volume 18 of Ergebnisse der Mathematik und ihrer Grenzgebiete (3) [Results in Mathematics and Related Areas (3)].  Springer-Verlag, Berlin, 1989.
  * [4] Peter J. Cameron.  Coherent configurations, association schemes and permutation groups.  In Groups, combinatorics & geometry (Durham, 2001), pages 55–71. World Sci. Publ., River Edge, NJ, 2003.
  * [5] Charles J. Colbourn and Jeffrey H. Dinitz, editors.  Handbook of combinatorial designs.  Discrete Mathematics and its Applications (Boca Raton). Chapman & Hall/CRC, Boca Raton, FL, second edition, 2007.
  * [6] The GAP Group.  GAP – Groups, Algorithms, and Programming, Version 4.11.1, 2021\.
  * [7] C. Ramankutty Nair.  A new class of designs.  J. Amer. Statist. Assoc., 59:817–833, 1964.
  * [8] Eran Nevo.  Polyhedrons and PBIBDs from hyperbolic manifolds.  Geometriae Dedicata, pages 1–8, 2015.
  * [9] B. V. Shah.  A generalisation of partially balanced incomplete block designs.  Ann. Math. Statist., 30:1041–1050, 1959.
  * [10] Bikas Kumar Sinha.  Some aspects of simplicity in the analysis of block designs.  J. Statist. Plann. Inference, 6(2):165–172, 1982.