k},n_{1}-1,n_{2}-1}\right)\,,$ (A.19) $U_{\mathrm{C},\alpha,{\bf
k}}=\sum_{n_{1},n_{2}}\sqrt{n_{1}(n_{2}+1)}\left(c^{0}_{n_{1}-1,n_{2}+1}u_{\alpha,{\bf
k},n_{1},n_{2}}+c^{0}_{n_{1},n_{2}}\,v_{\alpha,{\bf
k},n_{1}-1,n_{2}+1}\right)\,,$ (A.20) $V_{\mathrm{P},\alpha,{\bf
k}}=\sum_{n_{1},n_{2}}\sqrt{n_{1}\,n_{2}}\left(c^{0}_{n_{1},n_{2}}\,u_{\alpha,{\bf
k},n_{1}-1,n_{2}-1}+c^{0}_{n_{1}-1,n_{2}-2}\,v_{\alpha,{\bf
k},n_{1},n_{2}}\right)\,,$ (A.21) $V_{\mathrm{C},\alpha,{\bf
k}}=\sum_{n_{1},n_{2}}\sqrt{n_{1}(n_{2}+1)}\left(c^{0}_{n_{1},n_{2}}\,u_{\alpha,{\bf
k},n_{1}-1,n_{2}+1}+c^{0}_{n_{1}-1,n_{2}+2}\,v_{\alpha,{\bf
k},n_{1},n_{2}}\right)\,,$ (A.22)

in complete analogy with Eqs. (18).

Recalling the discussion on the completeness relation concerning the one-body
Bose fields (A.11), a natural question that arises from the decompositions in
Eqs. (A.16)-(A.17) is whether they reproduce the canonical commutators

$\left[\hat{a}_{1,{\bf r}}\,\hat{a}_{2,{\bf r}},\hat{a}^{\dagger}_{1,{\bf
s}}\,\hat{a}^{\dagger}_{2,{\bf s}}\right]=\delta_{\bf
r,s}\left[1+\hat{n}_{d}\right]\,,$ (A.23) $\left[\hat{a}_{1,{\bf
r}}\,\hat{a}^{\dagger}_{2,{\bf r}},\hat{a}^{\dagger}_{1,{\bf
s}}\,\hat{a}_{2,{\bf s}}\right]=-\delta_{\bf r,s}\,\hat{n}_{s}\,,$ (A.24)

where $\hat{n}_{d}$ and $\hat{n}_{s}$ were defined in Sec. 2.1. Inserting the
QGW expression of the pairing fields (A.18), the first-order approximation of
the commutation rules (A.23)-(A.24) reads

$\left[\delta_{1}\hat{\psi}_{\mathrm{P/C}}({\bf
r}),\delta_{1}\hat{\psi}^{\dagger}_{\mathrm{P/C}}({\bf
s})\right]=\frac{1}{I}\sum_{\bf k}e^{i\,{\bf k}\cdot({\bf r}-{\bf
s})}\sum_{\alpha}\left(\left|U_{\alpha,\mathrm{P/C},{\bf
k}}\right|^{2}-\left|V_{\alpha,\mathrm{P/C},{\bf k}}\right|^{2}\right)\,.$
(A.25)

Notably, we find numerically that, due to the only approximate description of
the excitation spectrum of the CFSF and PSF phases, in particular the flat
bands describing pairing and antipairing excitations in Figs. 3 and 6
respectively, the above canonical relations are not satisfied in general.
Likewise, we expect a similar result when the higher-order contributions are
included. For these reasons, only intraspecies coherence functions have been
considered in Sec. 4.1.

###  A.3 Second-order fluctuations

Having illustrated how first-order quantized excitations build on top of the
mean-field observables, we now estimate the corrections to observables due to
second-order quantum fluctuations.

Figure 17:  (a) Control parameter $F$ as defined in Eq. (A.26) plotted in two
dimensions as a function of $2\,d\,J/U$ for fixed $n_{0,d}$ corresponding to
the transitions shown in the right panel. The point located on each curve
indicates the location of the respective phase transition. (b)-(d) Quantum
corrections to the order parameters and total density for the transitions
considered in the (a). (b) First-order MI-to-SF transition for $U_{12}/U=0.9$
and $n_{d}=2$. (c) Second-order MI-to-SF transition for $U_{12}/U=0.5$ and
$n_{d}=2$. (d) CFSF-to-SF transition for $U_{12}/U=0.9$ and $n_{d}=1$. (e)
PSF-to-SF transition for $U_{12}/U=-0.7$ and $n_{d}\approx 1.47$. (f) Quantum
corrections to the CFSF order parameter for fixed total density $n_{d}=1$ at
$U_{12}/U=0.9$ and along the $\mu/U=0.5$ line within the CFSF phase. (g)
Quantum corrections to the PSF order parameter for fixed total density
$n_{d}\approx 1.47$ at $U_{12}/U=-0.7$ and along the $(\mu/U)_{c}=-0.35$
critical line within the PSF phase.

(Control parameter of the theory) – Following Ref. [35], we start our analysis
by inspecting the magnitude of quantum fluctuations around the Gutzwiller
mean-field state, identified by the control parameter

$F=1-\langle\hat{A}^{2}({\bf
r})\rangle=\sum_{n_{1},n_{2}}\langle\delta\hat{c}^{\dagger}_{n_{1},n_{2}}({\bf
r})\,\delta\hat{c}_{n_{1},n_{2}}({\bf
r})\rangle=\frac{1}{I}\sum_{\alpha}\sum_{\bf
k}\sum_{n_{1},n_{2}}\left|v_{\alpha,{\bf k},n_{1},n_{2}}\right|^{2}\,.$ (A.26)

Within the QGW theory of the one-component Bose Hubbard model, it was found in
Ref. [35] that, for $d=3$ and fixed $\langle\hat{n}_{d}\rangle$, $F$ is peaked
at the MI-to-SF transition and approaches zero in the limiting regimes away
from the transition. For the present two-component case, we show results for
$d=2$ in Fig. 17(a) for fixed $n_{d}$ across the first and second-order MI-to-
SF transitions [blue and orange lines, respectively], the PSF-to-SF transition
[purple line], and the CFSF-to-SF transition [yellow line].

In fixing the total filling and approaching the deep SF regime, one travels
along the line $\mu\sim-2\,d\,J$ for increasing $J$, such that the limit
$n_{d}\,U\to 0$ is reached. In three dimensions, this limit coincides with the
dilute, weakly-interacting regime, as the two-body scattering amplitude is
density-independent. On the contrary, from renormalization group arguments it
is known that in two dimensions the dilute limit is more subtle, since the
two-body scattering amplitude in the weakly-interacting regime depends
logarithmically on the diluteness parameter $n_{d}\,a^{2}_{\mathrm{2D}}$,
where $a_{\mathrm{2D}}$ is the two-dimensional scattering length [78, 79, 80,
33]. The bare interaction strength $U$ used throughout our analysis cannot
capture this behavior. Consequently, we do not extend our study of quantum
fluctuations for $d=2$ beyond $2\,d\,J/U\sim 1$.

The pathologies that arise when the dilute limit is approached in our $d=2$
QGW theory can be seen in Fig. 17(a). Here, it is shown how the control
parameter increases monotonically for large $2\,d\,J/U$, which strongly
contrasts with the monotonic decay seen in the $d=3$ QGW study of the one-
component BH model [35]. Although the mean-field description provides the
usual one-body coherent state in this limit [29], the growth of $F$ indicates
that this is not the exact ground state of the dilute gas. Nonetheless, for
smaller values of $2\,d\,J/U$, including the critical regions and the deep MI,
CFSF and PSF regimes, we observe in Fig. 17(b)-(e) that $F$ remains small,
which supports the overall reliability of our approach.

(Quantum corrections) – Now that we identified the range of validity of the
QGW method, we turn to study the quantum corrections to the local order
parameters and densities. Expanding the corresponding Gutzwiller operators to
second order in the fluctuations, these quantities are modified as

$\langle\hat{\mathcal{N}}_{i}({\bf
r})\rangle=n_{0,i}+\langle\delta_{2}\hat{\mathcal{N}}_{i}({\bf r})\rangle\,,$
(A.27) $\langle\hat{\mathcal{D}}_{i}({\bf
r})\rangle=d_{0,i}+\langle\delta_{2}\hat{\mathcal{D}}_{i}({\bf r})\rangle\,,$
(A.28) $\langle\hat{\psi}_{i}({\bf
r})\rangle=\psi_{0,i}+\langle\delta_{2}\hat{\psi}_{i}({\bf r})\rangle\,,$
(A.29) $\langle\hat{\psi}_{\mathrm{P}}({\bf
r})\rangle=\psi_{0,\mathrm{P}}+\langle\delta_{2}\hat{\psi}_{\mathrm{P}}({\bf
r})\rangle\,,$ (A.30) $\langle\hat{\psi}_{\mathrm{C}}({\bf
r})\rangle=\psi_{0,\mathrm{C}}+\langle\delta_{2}\hat{\psi}_{\mathrm{C}}({\bf
r})\rangle\,.$ (A.31)

where we note that the first-order fluctuations analysed in App. A.2 make a
vanishing contribution and hence do not correct the mean-field predictions.
The second-order terms appearing on the right-hand side of Eqs. (A.27)-(A.31)
are explicitly given by

$\langle\delta_{2}\hat{\mathcal{N}}_{i}({\bf
r})\rangle=-F\,n_{0,i}+\sum_{n_{1},n_{2}}\left(n_{1}\,\delta_{i,1}+n_{2}\,\delta_{i,2}\right)\langle\delta\hat{c}^{\dagger}_{n_{1},n_{2}}({\bf
r})\,\delta\hat{c}_{n_{1},n_{2}}({\bf r})\rangle\,,$ (A.32)
$\langle\delta_{2}\hat{\mathcal{D}}_{i}({\bf
r})\rangle=-F\,d_{0,i}+\sum_{n_{1},n_{2}}\left(n^{2}_{1}\,\delta_{i,1}+n^{2}_{2}\,\delta_{i,2}\right)\langle\delta\hat{c}^{\dagger}_{n_{1},n_{2}}({\bf
r})\,\delta\hat{c}_{n_{1},n_{2}}({\bf r})\rangle\,,$ (A.33)
$\langle\delta_{2}\hat{\psi}_{1}({\bf
r})\rangle=-F\,\psi_{0,1}+\sum_{n_{1},n_{2}}\sqrt{n_{1}}\,\langle\delta\hat{c}^{\dagger}_{n_{1}-1,n_{2}}({\bf
r})\,\delta\hat{c}_{n_{1},n_{2}}({\bf r})\rangle\,,$ (A.34)
$\langle\delta_{2}\hat{\psi}_{2}({\bf
r})\rangle=-F\,\psi_{0,2}+\sum_{n_{1},n_{2}}\sqrt{n_{2}}\,\langle\delta\hat{c}^{\dagger}_{n_{1},n_{2}-1}({\bf
r})\,\delta\hat{c}_{n_{1},n_{2}}({\bf r})\rangle\,,$ (A.35)
$\displaystyle\langle\delta_{2}\hat{\psi}_{\mathrm{P}}({\bf r})\rangle_{c}=$
$\displaystyle-F\,\psi_{0,\mathrm{P}}+\sum_{n_{1},n_{2}}\sqrt{n_{1}\,n_{2}}\,\langle\delta\hat{c}^{\dagger}_{n_{1}-1,n_{2}-1}({\bf
r})\,\delta\hat{c}_{n_{1},n_{2}}({\bf r})\rangle$ (A.36)
$\displaystyle-\left[\psi_{0,2}\,\langle\delta_{2}\hat{\psi}_{1}({\bf
r})\rangle+\psi_{0,1}\,\langle\delta_{2}\hat{\psi}_{2}({\bf
r})\rangle\right]\,,$
$\displaystyle\langle\delta_{2}\hat{\psi}_{\mathrm{C}}({\bf r})\rangle_{c}=$
$\displaystyle-F\,\psi_{0,\mathrm{C}}+\sum_{n_{1},n_{2}}\sqrt{n_{1}(n_{2}+1)}\,\langle\delta\hat{c}^{\dagger}_{n_{1}-1,n_{2}+1}({\bf
r})\,\delta\hat{c}_{n_{1},n_{2}}({\bf r})\rangle$ (A.37)
$\displaystyle-\left[\psi^{*}_{0,2}\,\langle\delta_{2}\hat{\psi}_{1}({\bf
r})\rangle+\psi_{0,1}\,\langle\delta_{2}\hat{\psi}^{\dagger}_{2}({\bf
r})\rangle\right]\,,$

where the symbol $\langle\cdot\rangle_{c}$ on the left hand-side of Eqs.
(A.36)-(A.37) emphasizes that one-body correlations are subtracted on the
right hand side so that only genuine pairing/antipairing quantum fluctuations
are retained, in accordance with Eqs. (5). At zero temperature, the second-
order expectation values appearing in Eqs. (A.32)-(A.37) can be evaluated
straightforwardly by a generalization of the following two examples,

$\langle\delta_{2}\hat{\mathcal{N}}_{i}({\bf
r})\rangle=-F\,n_{0,i}+\frac{1}{I}\sum_{\alpha}\sum_{\bf
k}\sum_{n_{1},n_{2}}\left(n_{1}\,\delta_{i,1}+n_{2}\,\delta_{i,2}\right)\left|v_{\alpha,{\bf
k},n_{1},n_{2}}\right|^{2}\,,$ (A.38) $\langle\delta_{2}\hat{\psi}_{i}({\bf
r})\rangle=-F\,\psi_{0,i}+\frac{1}{I}\sum_{\alpha}\sum_{\bf
k}\sum_{n_{1},n_{2}}\left(\delta_{i,1}\,\sqrt{n_{1}}\,v_{\alpha,{\bf
k},n_{1}-1,n_{2}}+\delta_{i,2}\,\sqrt{n_{2}}\,v_{\alpha,{\bf
k},n_{1},n_{2}-1}\right)v_{\alpha,{\bf k},n_{1},n_{2}}\,.$ (A.39)

from which we obtain self-contained expressions for the one-species filling
and the one-body order parameter corrected by quantum fluctuations,

$\langle\hat{\mathcal{N}}_{i}({\bf
r})\rangle=\left(1-F\right)n_{0,i}+\frac{1}{I}\sum_{\alpha}\sum_{\bf
k}\sum_{n_{1},n_{2}}\left(n_{1}\,\delta_{i,1}+n_{2}\,\delta_{i,2}\right)\left|v_{\alpha,{\bf
k},n_{1},n_{2}}\right|^{2}\,,$ (A.40) $\langle\hat{\psi}_{1}({\bf
r})\rangle=\left(1-F\right)\psi_{0,1}+\frac{1}{I}\sum_{\alpha}\sum_{\bf
k}\sum_{n_{1},n_{2}}\sqrt{n_{1}}\,v_{\alpha,{\bf
k},n_{1}-1,n_{2}}\,v_{\alpha,{\bf k},n_{1},n_{2}}\,.$ (A.41)

Accordingly, the expectation value of the total density is approximated by
$n_{d}\approx n_{0,d}+\sum_{i}\langle\delta_{2}\hat{\mathcal{N}}_{i}({\bf
r})\rangle$, as taken into account in Secs. 3 and 4 when discussing the
results of the QGW approach.

It is interesting to observe that, within the QGW formalism, quantum-corrected
local observables are always given by the sum of two distinct terms, one given
by quantum fluctuations only and the other, proportional to the mean-field
average, deriving exclusively from the normalization operator via the
expectation value $\big{\langle}\hat{A}^{2}{\left({\bf
r}\right)}\big{\rangle}$ defining the control parameter $F$. This result makes
more explicit the physical role of $\hat{A}{\left({\bf r}\right)}$, which
accounts for the feedback of quantum fluctuations onto the Gutzwiller ground
state. Along these lines, we also remark that this contribution is of key
importance in giving accurate predictions for the local density and spin
correlations presented in Sec. 4.2.

The relative quantum corrections to the total density and the one-body order
parameters are shown in Fig. 17(b)-(e) for fixed $n_{d}$. In the deep CFSF,
PSF and MI phases, quantum fluctuations are always small. This remains true
also near the transition points, whereas the corrections grow as either the
collapse point or the limit $n_{d}\,U\to 0$ is approached. These increasingly
large corrections mirror the pathologies of this limit in two dimensions. It
is worth noticing that, in Fig. 17(b)-(e), the corrections to the one-body
order parameters appear to diverge on the brink of the phase transitions. We
expect that the self-consistent inclusion of these effects into the phase
diagram description would shift the CFSF, PSF and MI boundaries towards larger
hopping energies in agreement with the results of quantum Monte Carlo
simulations [35, 51].

The quantum corrections to the CFSF and PSF order parameters are shown in Fig.
17(f)-(g) for fixed $n_{d}$ and across the respective phase transitions to the
SF regime. In the deep CFSF and PSF phases, these corrections are small, which
indicates that the effect of fluctuations in these phases within the QGW
picture is minimal. The corrections remain of the order of $10^{-1}$ across
the transitions and then monotonically increase as the pathological limit
$n_{d}\,U\to 0$ is reached. However, we note that the amplitude of relative
quantum corrections is amplified by the fact that both $\psi_{0,\mathrm{C}}$
and $\psi_{0,\mathrm{P}}$ decrease rapidly for weak interactions – see panels
(c) of Fig. 1 and Fig. 4.

##  Appendix B Linear response formalism within the QGW theory

This appendix is dedicated to a detailed review of the linear response
formalism (c.f. Ref. [54]) applied to the QGW quantum theory. For this
purpose, we derive the relevant expressions for the density, spin and current
response functions; in particular, the latter objects are the basic quantities
from which the superfluid matrix is calculated in Sec. 3.2.

###  B.1 Linear response formalism

Suppose that a Hermitian, time-dependent perturbation is applied to the
Hamiltonian of the system at time $t=0$,

$\hat{H}_{\mathrm{p}}(t)=\hat{H}+\theta(t)\,\hat{G}(t)\,.$ (B.42)

Let us define $|\psi^{N}_{n}\rangle$ as the $n^{\text{th}}$ eigenstate of
$\hat{H}$ with a total number of particles equal to $N$. For $t<0$, the system
is in the true many-body ground state $|\psi^{N}_{0}\rangle$ and the time
evolution of the system is given by
$|\psi(t)\rangle=\hat{U}(t)|\psi^{N}_{0}\rangle$ where $\hat{U}(t)$ is the
evolution operator with boundary condition $\hat{U}(0)=1$. Assuming that the
perturbation $\hat{G}(t)$ is sufficiently small with respect to all the energy
scales of $\hat{H}$, $\hat{U}(t)$ can be approximated by the Born
approximation of the solution to the equation
$i\,\hbar\,\partial_{t}\hat{U}(t)=\hat{H}_{\mathrm{p}}(t)\,\hat{U}(t)$, namely

$\hat{U}(t)\approx
e^{-i\,\hat{H}\,t/\hbar}-\frac{i}{\hbar}\int_{0}^{t}d\tau\,e^{-i\,\hat{H}(t-\tau)/\hbar}\,\hat{G}(\tau)\,e^{-i\,\hat{H}\,\tau/\hbar}\,.$
(B.43)

Now, we wish to calculate the time evolution of some observable
$\langle\hat{F}(t)\rangle$, which we determine in the Heisenberg picture via
the relation $\hat{F}(t)=e^{i\,\hat{H}\,t}\,\hat{F}\,e^{-i\,\hat{H}\,t}$.
Evaluating $\langle\hat{F}(t)\rangle$ using Eq. (B.43), we find

$\left\langle\hat{F}(t)\right\rangle=\left\langle\psi^{N}_{0}\left|\hat{U}^{\dagger}(t)\,\hat{F}\,\hat{U}(t)\right|\psi^{N}_{0}\right\rangle\approx\left\langle\hat{F}\right\rangle_{\mathrm{eq}}-\frac{i}{\hbar}\int_{0}^{t}d\tau\left\langle\psi^{N}_{0}\left|\left[\hat{G}(\tau),\hat{F}(t-\tau)\right]\right|\psi^{N}_{0}\right\rangle\,,$
(B.44)

where the “eq” subscript indicates the equilibrium expectation value. The
kernel of the above integral is known as the time-domain response function

$\chi_{\hat{F},\hat{G}}(t)=-\frac{i}{\hbar}\theta(t)\left\langle\psi^{N}_{0}\left|\left[\hat{F}(t),\hat{G}(0)\right]\right|\psi^{N}_{0}\right\rangle\,.$
(B.45)

In frequency space, we have

$\chi_{\hat{F},\hat{G}}(\omega)=\lim_{\varepsilon\to
0^{+}}\int_{-\infty}^{\infty}dt\,e^{-\varepsilon\,t}\,e^{i\,\omega\,t}\,\chi_{\hat{F},\hat{G}}(t)\,,$
(B.46)

where the infinitesimal regularization parameter $\varepsilon\to 0^{+}$
ensures that at $t=-\infty$ the system is governed by the unperturbed
Hamiltonian $\hat{H}$. Written in terms of the perturbation operator and the
probed observable, we find then the general expression

$\chi_{\hat{F},\hat{G}}(\omega)=-\frac{i}{\hbar}\lim_{\varepsilon\to
0^{+}}\int_{0}^{\infty}dt\,e^{-\varepsilon\,t}\,e^{i\,\omega\,t}\left\langle\left[\hat{F}(t),\hat{G}(0)\right]\right\rangle\,.$
(B.47)

In the following subsections, we proceed to specialize Eq. (B.47) to the case
of the density and current response functions of the two-component BH model
under the QGW quantum theory. Although the linear fluctuations around the
mean-field Gutzwiller state would provide a similar result at lowest order,
the QGW formalism makes the calculation straightforward and particularly
transparent. Moreover, the higher-order correlations accessed by the QGW
quantum theory are shown to introduce the effect of multi-mode fluctuations
into the response functions: in particular, this last feature is at the roots
of our accurate predictions for the current response functions and, in turn,
of the superfluid components of the system.

###  B.2 Density response function

As a standard yet relevant case of study, we are interested in the linear
response of the ${\bf q}$-component of the density operator

$\hat{\rho}_{i,{\bf q}}=\sum_{\bf k}\hat{a}^{\dagger}_{i,{\bf
k-q}}\,\hat{a}_{i,{\bf k}}$ (B.48)

for the $i^{\text{th}}$ species under the same type of density perturbation.
Therefore, we set $\hat{F}=\hat{G}=\delta\,\hat{\rho}_{i,{\bf q}}$ with
$\delta\hat{\rho}_{i,{\bf q}}=\hat{\rho}_{i,{\bf
q}}-\left\langle\hat{\rho}_{i,{\bf q}}\right\rangle_{\mathrm{eq}}$, where the
equilibrium contribution vanishes in a uniform system unless at ${\bf q=0}$.
The density-density response function is given by rephrasing Eq. (B.47) into

$\chi_{\hat{n}_{i}}({\bf q},\omega)=-\frac{i}{\hbar}\lim_{\varepsilon\to
0^{+}}\int_{0}^{\infty}dt\,e^{-\varepsilon\,t}\,e^{i\,\omega\,t}\,\left\langle\left[\delta\hat{\rho}_{i,{\bf
q}}(t),\delta\hat{\rho}_{i,{\bf-q}}(0)\right]\right\rangle\,.$ (B.49)

In order to evaluate this response function within the QGW formalism, we apply
the usual quantization procedure to the expectation value in Eq. (B.49), which
at the lowest order in the fluctuations is mapped into

$\langle\left[\delta\hat{\rho}_{i,{\bf
q}}(t),\delta\hat{\rho}_{i,{\bf-q}}(0)\right]\rangle\longrightarrow\sum_{\bf
r}e^{-i\,{\bf q}\cdot{\bf
r}}\left\langle\left[\delta_{1}\hat{\mathcal{N}}_{i}({\bf
r},t),\delta_{1}\hat{\mathcal{N}}_{i}({\bf 0},0)\right]\right\rangle\,,$
(B.50)

where $\delta_{1}\hat{\mathcal{N}}_{i}({\bf r},t)$ is the first-order
expansion of the QGW density operator defined in Eq. (A.3), whose time-
dependence is controlled by the interaction picture of the QGW Hamiltonian
(42). Therefore, evaluating the commutator in Eq. (B.50), we readily obtain

$\chi_{\hat{n}_{i}}({\bf
q},\omega)=\frac{2}{\hbar}\,\sum_{\alpha}\frac{N^{2}_{i,\alpha,{\bf
q}}\,\omega_{\alpha,{\bf
q}}}{\left(\omega+i\,0^{+}\right)^{2}-\omega^{2}_{\alpha,{\bf q}}}\,.$ (B.51)

As anticipated, our result for $\chi_{\hat{n}_{i}}({\bf q},\omega)$ agrees
exactly with the linear-response result in Ref. [32] – see in particular Eq.
(29) and App. A therein – based on the time-dependent Gutzwiller
approximation, although our quantum formalism offers a more transparent and
straightforward tool for the estimation of response functions, as well as the
unexplored possibility of including non-Gaussian fluctuations in the density
(spin) channel.

###  B.3 Current response functions

####  B.3.1 Current operator

As a starting ground, we determine the species-resolved local current
operators by resorting to the continuity equation, which is given by

$\nabla_{\bf r}\,\hat{j}_{i,{\bf r}}=-\frac{\partial\hat{n}_{i,{\bf
r}}}{\partial t}=\frac{i}{\hbar}\left[\hat{n}_{i,{\bf r}},\hat{H}\right]\,.$
(B.52)

where $\nabla_{\bf r}$ has to be intended according to its definition on a
lattice. Inserting the BH Hamiltonian (2) on the right-hand side of Eq.
(B.52), one finds

$\hat{j}_{i,{\bf
r}}=\frac{i\,J\,a}{\hbar}\sum_{j=1}^{d}\left(\hat{a}^{\dagger}_{i,{\bf
r}+\bf{e}_{j}}\,\hat{a}_{i,{\bf r}}-\text{h.c.}\right)\bf{e}_{j}\,.$ (B.53)

For later convenience, we introduce the current operator acting along a
specific direction of the square lattice, for instance the $x$-directed links,
for which we have

$\hat{j}_{i,{\bf r}}|_{x}=\frac{i\,J\,a}{\hbar}\left(\hat{a}^{\dagger}_{i,{\bf
r}+\bf{e}_{x}}\,\hat{a}_{i,{\bf r}}-\text{h.c.}\right)\,.$ (B.54)

In momentum space, the unidirectional current operator Eq. (B.54) transforms
into

$\hat{j}_{i,{\bf q}}|_{x}=\frac{i\,J\,a}{\hbar}\sum_{\bf
k}\left[e^{-i\,k_{x}\,a}-e^{i(k_{x}+q_{x})a}\right]\hat{a}^{\dagger}_{i,{\bf
k}}\,\hat{a}_{i,{\bf k+q}}\,,$ (B.55)

which in the uniform limit ${\bf q\to 0}$ becomes

$\hat{j}_{i,{\bf q\to 0}}|_{x}=\frac{\hbar}{m\,a}\sum_{\bf
k}\sin(k_{x}\,a)\,\hat{a}^{\dagger}_{i,{\bf k}}\,\hat{a}_{i,{\bf k}}\,.$
(B.56)

The uniform current operator derived above is the fundamental object by which
we estimate the relevant current response functions in the following
subsections.

####  B.3.2 Intraspecies response

First, we consider the case where both the probe and the response operators
correspond to the same current species $\hat{j}_{i,{\bf q}}$. Let us take the
momentum vector ${\bf q}$ to lie along the $z$ axis, then the longitudinal and
transverse components of the current operator are given by $\hat{j}_{i,{\bf
q}}|_{z}$ and $\hat{j}_{i,{\bf q}}|_{x}$, respectively. In the following, we
consider only the static ($\omega=0$) and uniform (${\bf q\to 0}$) limits of
the transverse response function [81, 82], which represents the paramagnetic
contribution to the superfluid matrix discussed in Sec. 3.2. The relevant
response function is

$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}({\bf q\to
0},\omega=0)=-\frac{i}{\hbar}\lim_{\varepsilon\to
0^{+}}\int_{0}^{\infty}dt\,e^{-\varepsilon\,t}\,\left\langle\left[\hat{j}_{i}|_{x}(t),\hat{j}_{i}|_{x}(0)\right]\right\rangle\,.$
(B.57)

Starting from Eq. (B.56) and making use of the mapping (A.10), the QGW
quantization of the uniform current operator $\hat{j}_{i,{\bf q\to
0}}(t)|_{x}$ is found to be

$\displaystyle\hat{\mathcal{J}}_{i}|_{x}=\frac{\hbar}{m\,a}\sum_{\alpha,\beta}\sum_{\bf
k}\left[U^{*}_{i,\alpha,{\bf k}}\,\hat{b}^{\dagger}_{\alpha,{\bf
k}}+V_{i,\alpha,{\bf k}}\,\hat{b}_{\alpha,-{\bf
k}}\right]\left[U_{i,\beta,{\bf k}}\,\hat{b}_{\beta,{\bf
k}}+V^{*}_{i,\beta,{\bf k}}\,\hat{b}^{\dagger}_{\beta,{-\bf
k}}\right]\sin(k_{x}\,a)\,,$ (B.58)

in which the presence of the condensate plays no role as $\sin(k_{x}\,a)$
vanishes at ${\bf k=0}$. Next, on the right-hand side of Eq. (B.57) we insert
a complete basis of excited states as a resolution of the identity operator to
obtain

$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}({\bf
0},0)=-\frac{i}{\hbar}\lim_{\varepsilon\to
0^{+}}\int_{0}^{\infty}dt\,e^{-\varepsilon\,t}\sum_{N,n>0}\left[\left\langle\psi^{N}_{0}\left|\hat{\mathcal{J}}_{i}|_{x}(t)\right|\psi^{N}_{n}\right\rangle\left\langle\psi^{N}_{n}\left|\hat{\mathcal{J}}_{i}|_{x}(0)\right|\psi^{N}_{0}\right\rangle-\mathrm{c.c.}\right]\,.$
(B.59)

At zero temperature, we obtain that the first-order term of the QGW current
operator (B.58) produces a vanishing contribution to the transverse current
response. Therefore, we immediately recognize that only expectation values of
the kind $\langle\psi^{N}_{0}|\,\hat{b}_{\alpha,{\bf k}}\,\hat{b}_{\beta,{\bf
p}}\,|\psi^{N+2}_{n}\rangle$ and $\langle\psi^{N}_{0}|\,\hat{b}_{\alpha,{\bf
k}}\,\hat{b}_{\beta,{\bf p}}\,|\psi^{N+2}_{n}\rangle$ provide a finite result.
More precisely, the second-order of $\hat{\mathcal{J}}_{i}(t)|_{x}$ generating
such contributions can be written in a symmetric form as

$\displaystyle\sum_{\alpha,\beta}\sum_{\bf k}U_{i,\alpha,{\bf
k}}\,V_{i,\beta,{\bf k}}\,\sin(k_{x}\,a)\,\hat{b}_{\alpha,{\bf
k}}\,\hat{b}_{\beta,{\bf k}}$ (B.60)
$\displaystyle=\frac{1}{2}\sum_{\alpha,\beta}\sum_{\bf
k}\left[U_{i,\alpha,{\bf k}}\,V_{i,\alpha,{\bf k}}-U_{i,\beta,{\bf
k}}\,V_{i,\alpha,{\bf k}}\right]\sin(k_{x}\,a)\,\hat{b}_{\alpha,{\bf
k}}\,\hat{b}_{\beta,{\bf k}}\,,$

where we have used the fact that the particle (hole) amplitudes
$U_{i,\alpha,{\bf k}}$ $\left(V_{i,\alpha,{\bf k}}\right)$ are even functions
of momentum in our system. Plugging the right-hand side of Eq. (B.60) back
into the linear response function (B.59) and computing the time integral, we
find

$\chi^{T}_{\hat{j}_{i},\hat{j}_{i}}({\bf
0},0)=-\frac{\hbar}{\left(m\,a\right)^{2}}\frac{2^{2}}{4}\sum_{\alpha,\beta}\sum_{\bf
k}\frac{\left(U_{i,\alpha,{\bf k}}\,V_{i,\beta,{\bf k}}-U_{i,\beta,{\bf
k}}\,V_{i,\alpha,{\bf k}}\right)^{2}}{\omega_{\alpha,{\bf
k}}+\omega_{\beta,{\bf k}}}\sin^{2}(k_{x}\,a)\,,$ (B.61)

where we have neglected the conjugation of the particle (hole) amplitudes
$U_{i,\alpha,{\bf k}}$ $\left(V_{i,\alpha,{\bf k}}\right)$ for simplicity, as
they turn out to be always real in the present case. In the above equation,
although the numerical factors cancel each other, we have chosen to write them
explicitly in order to emphasize their origins. The factor of $1/4$ arises
from the symmetrization of the current operator in Eq. (B.60); one factor of
$2$ descends from the inner product $\langle\psi^{N}_{0}|\hat{b}_{\alpha,{\bf
k}}\,\hat{b}_{\beta,{\bf p}}|\psi^{N+2}_{n}\rangle$ (appearing twice), while
the other factor of $2$ comes from the complex conjugate term in Eq. (B.59).

####  B.3.3 Interspecies response

Now, we consider the case in which the probe is proportional to the current of
one species while we measure the linear response of the current of the other
species, in the transverse direction and in the static-uniform limit as
before. This corresponds to the off-diagonal response function

$\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}({\bf q\to
0},\omega=0)=-\frac{i}{\hbar}\lim_{\varepsilon\to
0^{+}}\int_{0}^{\infty}dt\,e^{-\varepsilon\,t}\,\left\langle\left[\hat{j}_{1}|_{x}(t),\hat{j}_{2}|_{x}(0)\right]\right\rangle\,.$
(B.62)

The calculation of $\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}({\bf q\to 0},\omega=0)$
within the QGW framework can be performed analogously to the intraspecies case
by considering the current operator (B.58) for both the bosonic flavours. The
final result is

$\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}({\bf
0},0)=-\frac{\hbar}{\left(m\,a\right)^{2}}\frac{2^{2}}{4}\sum_{\alpha,\beta}\sum_{\bf
k}\frac{\prod_{i=1}^{2}\left(U_{i,\alpha,{\bf k}}\,V_{i,\beta,{\bf
k}}-U_{i,\beta,{\bf k}}\,V_{i,\alpha,{\bf k}}\right)}{\omega_{\alpha,{\bf
k}}+\omega_{\beta,{\bf k}}}\sin^{2}(k_{x}\,a)\,,$ (B.63)

where the numerical factors have the same interpretations as for Eq. (B.61).

It is instructive to briefly comment on how Eq. (B.63) fundamentally differs
from the equivalent result of Ref. [27] within the Bogoliubov approximation.
In the first place, the expression of $\chi^{T}_{\hat{j}_{1},\hat{j}_{2}}({\bf
q\to 0},\omega=0)$ in that work contains an additional factor of $1/2$. This
is due to the fact that the $U$’s and $V$’s for the spin and density Goldstone
modes are each normalized to one in Ref. [27], whereas in the present case the
normalization condition for the $U$’s and $V$’s is given by the sum over all
the excitation branches (see the relative discussion in App. A.2), which
typically takes more than two branches to saturate numerically. Furthermore,
the two results differ by a minus sign. When only the first two excitation
branches are considered in the QGW method, one finds that the $U$’s and $V$’s
for the density Goldstone mode differ by a minus sign between the two species,
whereas they are identical for the spin mode. This results in an overall minus
sign when only these two branches are considered in Eq. (B.63). However, we
notice that the normalization condition is far from being saturated if limited
to the Goldstone modes, so a good agreement with the Bogoliubov result is
expected only in the weakly-interacting regime. These remarks resolve the
apparent discrepancies between the Bogoliubov and QGW predictions.

###  B.4 Average kinetic energy

Albeit not being a proper response function, we conclude by calculating the
mean value of the average kinetic energy along $x$-directed links of the
lattice, $K_{i}|_{x}$, which is exactly proportional to the density $n_{i}$ in
the well-known expression for the superfluid matrix on the continuum [27]. The
local kinetic energy operator for the $i^{\text{th}}$ species along the $x$
direction is defined as

$\hat{K}_{i}|_{x}({\bf r})=-J\left(\hat{a}^{\dagger}_{{\bf
r}+\bf{e}_{x},i}\,\hat{a}_{{\bf r},i}+\mathrm{h.c.}\right)\,.$ (B.64)

Within the mean-field Gutzwiller approximation, one obtains $K_{0,i}|_{x}({\bf
r})=-2\,J\left|\psi_{0,i}\right|^{2}$, which is trivially proportional to the
condensate fraction. To improve this physical description of the kinetic
energy, we evaluate the average kinetic energy by rewriting the operator
(B.64) in terms of the QGW Bose fields (A.7)-(A.8). It follows that, up to the
lowest-order in the quantum fluctuations, the QGW local kinetic energy
operator reads

$\displaystyle\hat{\mathcal{K}}_{i}|_{x}({\bf r})\approx$
$\displaystyle-2\,J\left|\psi_{0,i}({\bf r})\right|^{2}-J\left[\psi_{0,i}({\bf
r})\,\delta_{1}\hat{\psi}^{\dagger}_{i}({\bf r}+{\bf
e}_{x})+\psi^{*}_{0,i}({\bf r}+{\bf e}_{x})\,\delta_{1}\hat{\psi}_{i}({\bf
r})+\mathrm{h.c.}\right]$ (B.65)
$\displaystyle-J\left[\delta_{1}\hat{\psi}^{\dagger}_{i}({\bf r}+{\bf
e}_{x})\,\delta_{1}\hat{\psi}_{i}({\bf r})+\mathrm{h.c.}\right]\,.$

For a uniform lattice, the average value of $\hat{\mathcal{K}}_{i}|_{x}({\bf
r})$ is given by

$K_{i}|_{x}\approx-2\,J\left|\psi_{0,i}\right|^{2}-\frac{2\,J}{I}\sum_{\alpha}\sum_{\bf
k}\left|V_{i,\alpha,{\bf k}}\right|^{2}\,\cos(k_{x}\,a)\,,$ (B.66)

where we have made use of Eq. (A.9) and neglected finite-temperature
contributions. We notice that only the first term on the right-hand side of
Eq. (B.66) agrees with the prediction of mean-field theory, while the second
term is due to zero-point fluctuations properly accounted by our quantum
theory. In this regard, the expression of $K_{i}|_{x}$ recalls the average
kinetic energy of the one-component system derived in Ref. [35].

## References

  * [1] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner and P. Zoller,  _Cold bosonic atoms in optical lattices_ ,  Phys. Rev. Lett. 81, 3108 (1998),  10.1103/PhysRevLett.81.3108.
  * [2] M. P. A. Fisher, P. B. Weichman, G. Grinstein and D. S. Fisher,  _Boson localization and the superfluid-insulator transition_ ,  Phys. Rev. B 40, 546 (1989),  10.1103/PhysRevB.40.546.
  * [3] A. Altland and B. D. Simons,  _Condensed matter field theory_ ,  Cambridge university press (2010).
  * [4] M. Greiner, O. Mandel, T. Esslinger, T. W. Hänsch and I. Bloch,  _Quantum phase transition from a superfluid to a Mott insulator in a gas of ultracold atoms_ ,  Nature 415(6867), 39 (2002),  10.1038/415039a.
  * [5] E. Altman, W. Hofstetter, E. Demler and M. D. Lukin,  _Phase diagram of two-component bosons on an optical lattice_ ,  New Journal of Physics 5, 113 (2003),  10.1088/1367-2630/5/1/113.
  * [6] A. Kuklov, N. Prokof’ev and B. Svistunov,  _Commensurate two-component bosons in an optical lattice: ground state phase diagram_ ,  Phys. Rev. Lett. 92, 050402 (2004),  10.1103/PhysRevLett.92.050402.
  * [7] A. B. Kuklov and B. V. Svistunov,  _Counterflow superfluidity of two-species ultracold atoms in a commensurate optical lattice_ ,  Phys. Rev. Lett. 90, 100401 (2003),  10.1103/PhysRevLett.90.100401.
  * [8] A. Isacsson, M.-C. Cha, K. Sengupta and S. M. Girvin,  _Superfluid-insulator transitions of two-species bosons in an optical lattice_ ,  Phys. Rev. B 72, 184507 (2005),  10.1103/PhysRevB.72.184507.
  * [9] A. Hubener, M. Snoek and W. Hofstetter,  _Magnetic phases of two-component ultracold bosons in an optical lattice_ ,  Phys. Rev. B 80, 245109 (2009),  10.1103/PhysRevB.80.245109.
  * [10] A. Hu, L. Mathey, I. Danshita, E. Tiesinga, C. J. Williams and C. W. Clark,  _Counterflow and paired superfluidity in one-dimensional Bose mixtures in optical lattices_ ,  Phys. Rev. A 80, 023619 (2009),  10.1103/PhysRevA.80.023619.
  * [11] A. Andreev and E. Bashkin,  _Three-velocity hydrodynamics of superfluid solutions_ ,  JETP 42, 164 (1976).
  * [12] M. A. Alpar, S. A. Langer and J. A. Sauls,  _Rapid postglitch spin-up of the superfluid core in pulsars._ ,  The Astrophysical Journal 282, 533 (1984),  10.1086/162232.
  * [13] O. Sjöberg,  _On the landau effective mass in asymmetric nuclear matter_ ,  Nuclear Physics A 265(3), 511 (1976),  https://doi.org/10.1016/0375-9474(76)90558-3.
  * [14] M. Borumand, R. Joynt and W. Kluźniak,  _Superfluid densities in neutron-star matter_ ,  Phys. Rev. C 54, 2745 (1996),  10.1103/PhysRevC.54.2745.
  * [15] M. Baldo, J. Cugnon, A. Lejeune and U. Lombardo,  _Proton and neutron superfluidity in neutron star matter_ ,  Nuclear Physics A 536(2), 349 (1992),  https://doi.org/10.1016/0375-9474(92)90387-Y.
  * [16] N. Chamel,  _Two-fluid models of superfluid neutron star cores_ ,  Monthly Notices of the Royal Astronomical Society 388(2), 737 (2008),  10.1111/j.1365-2966.2008.13426.x,  https://academic.oup.com/mnras/article-pdf/388/2/737/3006118/mnras0388-0737.pdf.
  * [17] S. B. Chung, H. Bluhm and E.-A. Kim,  _Stability of half-quantum vortices in ${p}_{x}+i{p}_{y}$ superconductors_,  Phys. Rev. Lett. 99, 197002 (2007),  10.1103/PhysRevLett.99.197002.
  * [18] S. B. Chung, D. F. Agterberg and E.-A. Kim,  _Fractional vortex lattice structures in spin-triplet superconductors_ ,  New Journal of Physics 11(8), 085004 (2009),  10.1088/1367-2630/11/8/085004.
  * [19] J. Garaud, K. A. H. Sellin, J. Jäykkä and E. Babaev,  _Skyrmions induced by dissipationless drag in U(1) $\times{}$U(1) superconductors_,  Phys. Rev. B 89, 104508 (2014),  10.1103/PhysRevB.89.104508.
  * [20] D. V. Fil and S. I. Shevchenko,  _Nondissipative drag of superflow in a two-component Bose gas_ ,  Phys. Rev. A 72, 013616 (2005),  10.1103/PhysRevA.72.013616.
  * [21] J. Linder and A. Sudbø,  _Calculation of drag and superfluid velocity from the microscopic parameters and excitation energies of a two-component Bose-Einstein condensate in an optical lattice_ ,  Phys. Rev. A 79, 063610 (2009),  10.1103/PhysRevA.79.063610.
  * [22] J. Nespolo, G. E. Astrakharchik and A. Recati,  _Andreev–Bashkin effect in superfluid cold gases mixtures_ ,  New Journal of Physics 19(12), 125005 (2017),  10.1088/1367-2630/aa93a0.
  * [23] L. Parisi, G. E. Astrakharchik and S. Giorgini,  _Spin dynamics and Andreev-Bashkin effect in mixtures of one-dimensional Bose gases_ ,  Phys. Rev. Lett. 121, 025302 (2018),  10.1103/PhysRevLett.121.025302.
  * [24] V. Karle, N. Defenu and T. Enss,  _Coupled superfluidity of binary Bose mixtures in two dimensions_ ,  Phys. Rev. A 99, 063627 (2019),  10.1103/PhysRevA.99.063627.
  * [25] K. Sellin and E. Babaev,  _Superfluid drag in the two-component Bose-Hubbard model_ ,  Phys. Rev. B 97, 094517 (2018),  10.1103/PhysRevB.97.094517.
  * [26] D. Contessi, D. Romito, M. Rizzi and A. Recati,  _Collisionless drag for a one-dimensional two-component Bose-Hubbard model_ ,  Phys. Rev. Research 3, L022017 (2021),  10.1103/PhysRevResearch.3.L022017.
  * [27] D. Romito, C. Lobo and A. Recati,  _Linear response study of collisionless spin drag_ ,  Phys. Rev. Research 3, 023196 (2021),  10.1103/PhysRevResearch.3.023196.
  * [28] E. Blomquist, A. Syrwid and E. Babaev,  _Borromean supercounterfluidity_ (2021), 2108.11935.
  * [29] K. V. Krutitsky and P. Navez,  _Excitation dynamics in a lattice Bose gas within the time-dependent Gutzwiller mean-field approach_ ,  Phys. Rev. A 84, 033602 (2011),  10.1103/PhysRevA.84.033602.
  * [30] W. Krauth, M. Caffarel and J.-P. Bouchaud,  _Gutzwiller wave function for a model of strongly interacting bosons_ ,  Phys. Rev. B 45, 3137 (1992),  10.1103/PhysRevB.45.3137.
  * [31] K. V. Krutitsky,  _Ultracold bosons with short-range interaction in regular optical lattices_ ,  Physics Reports 607, 1 (2016),  https://doi.org/10.1016/j.physrep.2015.10.004,  Ultracold bosons with short-range interaction in regular optical lattices.
  * [32] T. Ozaki, I. Danshita and T. Nikuni,  _Bose-Bose mixtures in an optical lattice: first-order superfluid-insulator transition and elementary excitations_ (2012), 1210.1370.
  * [33] L. Pitaevskii and S. Stringari,  _Bose-Einstein condensation and superfluidity_ , vol. 164,  Oxford University Press (2016).
  * [34] A. M. Rey, K. Burnett, R. Roth, M. Edwards, C. J. Williams and C. W. Clark,  _Bogoliubov approach to superfluidity of atoms in an optical lattice_ ,  Journal of Physics B: Atomic, Molecular and Optical Physics 36(5), 825 (2003),  10.1088/0953-4075/36/5/304.
  * [35] F. Caleffi, M. Capone, C. Menotti, I. Carusotto and A. Recati,  _Quantum fluctuations beyond the Gutzwiller approximation in the Bose-Hubbard model_ ,  Phys. Rev. Research 2, 033276 (2020),  10.1103/PhysRevResearch.2.033276.
  * [36] F. Caleffi, M. Capone, I. de Vega and A. Recati,  _Impurity dephasing in a bose–hubbard model_ ,  New Journal of Physics 23(3), 033018 (2021),  10.1088/1367-2630/abe080.
  * [37] D. B. M. Dickerscheid, D. van Oosten, P. J. H. Denteneer and H. T. C. Stoof,  _Ultracold atoms in optical lattices_ ,  Phys. Rev. A 68, 043623 (2003),  10.1103/PhysRevA.68.043623.
  * [38] D. Huerga, J. Dukelsky and G. E. Scuseria,  _Composite boson mapping for lattice boson systems_ ,  Phys. Rev. Lett. 111, 045701 (2013),  10.1103/PhysRevLett.111.045701.
  * [39] I. Frérot and T. Roscilde,  _Entanglement entropy across the superfluid-insulator transition: a signature of bosonic criticality_ ,  Phys. Rev. Lett. 116, 190401 (2016),  10.1103/PhysRevLett.116.190401.
  * [40] J. Catani, L. De Sarlo, G. Barontini, F. Minardi and M. Inguscio,  _Degenerate bose-bose mixture in a three-dimensional optical lattice_ ,  Phys. Rev. A 77, 011603 (2008),  10.1103/PhysRevA.77.011603.
  * [41] B. Gadway, D. Pertot, R. Reimann and D. Schneble,  _Superfluidity of interacting bosonic mixtures in optical lattices_ ,  Phys. Rev. Lett. 105, 045303 (2010),  10.1103/PhysRevLett.105.045303.
  * [42] T. Fukuhara, P. Schauß, M. Endres, S. Hild, M. Cheneau, I. Bloch and C. Gross,  _Microscopic observation of magnon bound states and their dynamics_ ,  Nature 502(7469), 76 (2013),  10.1038/nature12541.
  * [43] P. A. Lee, N. Nagaosa and X.-G. Wen,  _Doping a mott insulator: Physics of high-temperature superconductivity_ ,  Rev. Mod. Phys. 78, 17 (2006),  10.1103/RevModPhys.78.17.
  * [44] P. Ao and S. T. Chui,  _Binary Bose-Einstein condensate mixtures in weakly and strongly segregated phases_ ,  Phys. Rev. A 58, 4836 (1998),  10.1103/PhysRevA.58.4836.
  * [45] P. Chen and M.-F. Yang,  _Quantum phase transitions in a two-species hard-core boson Hubbard model in two dimensions_ ,  Phys. Rev. B 82, 180510 (2010),  10.1103/PhysRevB.82.180510.
  * [46] S. D. Huber, E. Altman, H. P. Büchler and G. Blatter,  _Dynamical properties of ultracold bosons in an optical lattice_ ,  Phys. Rev. B 75, 085106 (2007),  10.1103/PhysRevB.75.085106.
  * [47] D. Podolsky and S. Sachdev,  _Spectral functions of the Higgs mode near two-dimensional quantum critical points_ ,  Phys. Rev. B 86, 054508 (2012),  10.1103/PhysRevB.86.054508.
  * [48] D. Podolsky, A. Auerbach and D. P. Arovas,  _Visibility of the amplitude (Higgs) mode in condensed matter_ ,  Phys. Rev. B 84, 174522 (2011),  10.1103/PhysRevB.84.174522.
  * [49] D. Pekker and C. Varma,  _Amplitude/Higgs Modes in Condensed Matter Physics_ ,  Annual Review of Condensed Matter Physics 6(1), 269 (2015),  10.1146/annurev-conmatphys-031214-014350.
  * [50] D. Yamamoto, T. Ozaki, C. A. R. Sá de Melo and I. Danshita,  _First-order phase transition and anomalous hysteresis of Bose gases in optical lattices_ ,  Phys. Rev. A 88, 033624 (2013),  10.1103/PhysRevA.88.033624.
  * [51] B. Capogrosso-Sansone, i. m. c. G. m. c. Söyler, N. Prokof’ev and B. Svistunov,  _Monte Carlo study of the two-dimensional Bose-Hubbard model_ ,  Phys. Rev. A 77, 015602 (2008),  10.1103/PhysRevA.77.015602.
  * [52] A. J. Leggett,  _Quantum liquids: Bose condensation and Cooper pairing in condensed-matter systems_ ,  Oxford university press (2006).
  * [53] C. Cohen-Tannoudji, J. Dupont-Roc and G. Grynberg,  _Photons and atoms: introduction to quantum electrodynamics_ ,  Wiley-VCH Verlag GmbH & Co. KGaA, Weinheim, Germany (1997).
  * [54] J.-P. Blaizot and G. Ripka,  _Quantum theory of finite systems_ ,  MIT Press, Cambridge, Massachusetts (1985).
  * [55] C. J. Pethick and H. Smith,  _Bose-Einstein condensation in dilute gases_ ,  Cambridge University Press, 2 edn.,  10.1017/CBO9780511802850 (2008).
  * [56] Y. Castin,  _Coherent atomic matter waves_ ,  Lecture Notes of Les Houches Summer School, édité par CWR Kaiser et F. David (Springer-Verlag, Berlin/Heidelberg, 2001) (2001).
  * [57] Y. Castin and R. Dum,  _Low-temperature Bose-Einstein condensates in time-dependent traps: beyond the $U(1)$ symmetry-breaking approach_,  Phys. Rev. A 57, 3008 (1998),  10.1103/PhysRevA.57.3008.
  * [58] A. J. Leggett,  _On the superfluid fraction of an arbitrary many-body system at T = 0_ ,  Journal of Statistical Physics 93(3/4), 927 (1998),  10.1023/b:joss.0000033170.38619.6c.
  * [59] F. Carlini and S. Stringari,  _Spin drag and fast response in a quantum mixture of atomic gases_ ,  Phys. Rev. A 104, 023301 (2021),  10.1103/PhysRevA.104.023301.
  * [60] S. Sachdev,  _Quantum phase transitions_ ,  Handbook of Magnetism and Advanced Magnetic Materials (2007).
  * [61] M. E. Pezzoli, F. Becca, M. Fabrizio and G. Santoro,  _Local moments and magnetic order in the two-dimensional anderson-mott transition_ ,  Phys. Rev. B 79, 033111 (2009),  10.1103/PhysRevB.79.033111.
  * [62] S. Zhou, L. Liang and Z. Wang,  _Dynamical slave-boson mean-field study of the Mott transition in the Hubbard model in the large- $z$ limit_,  Phys. Rev. B 101, 035106 (2020),  10.1103/PhysRevB.101.035106.
  * [63] A. Recati and F. Piazza,  _Breaking of Goldstone modes in a two-component Bose-Einstein condensate_ ,  Phys. Rev. B 99, 064505 (2019),  10.1103/PhysRevB.99.064505.
  * [64] F. Caleffi,  _Quantum fluctuations beyond the Gutzwiller approximation in the Bose-Hubbard model_ ,  Master’s thesis, University of Trento (2018).
  * [65] A. Syrwid, E. Blomquist and E. Babaev,  _Dissipationless vector drag - Superfluid spin Hall effect_ ,  Phys. Rev. Lett. 127, 100403 (2021),  10.1103/PhysRevLett.127.100403.
  * [66] P. Fazekas and K. Penc,  _Generalized Gutzwiller ansatz for the half-filled Hubbard chain_ ,  International Journal of Modern Physics B 2(5), 1021 (1988),  10.1142/S0217979288000858,  https://doi.org/10.1142/S0217979288000858.
  * [67] P. Fazekas,  _Variational study of the pair hopping model_ ,  International Journal of Modern Physics B 3(12), 1765 (1989),  10.1142/S0217979289001135,  https://doi.org/10.1142/S0217979289001135.
  * [68] D.-S. Lühmann,  _Cluster Gutzwiller method for bosonic lattice systems_ ,  Phys. Rev. A 87, 043619 (2013),  10.1103/PhysRevA.87.043619.
  * [69] K. Suthar, H. Sable, R. Bai, S. Bandyopadhyay, S. Pal and D. Angom,  _Supersolid phase of the extended Bose-Hubbard model with an artificial gauge field_ ,  Phys. Rev. A 102, 013320 (2020),  10.1103/PhysRevA.102.013320.
  * [70] S. Mondal, A. Kshetrimayum and T. Mishra,  _Two-body repulsive bound pairs in a multibody interacting Bose-Hubbard model_ ,  Phys. Rev. A 102, 023312 (2020),  10.1103/PhysRevA.102.023312.
  * [71] M. Malakar, S. Ray, S. Sinha and D. Angom,  _Phases and collective modes of bosons in a triangular lattice at finite temperature: a cluster mean field study_ ,  Phys. Rev. B 102, 184515 (2020),  10.1103/PhysRevB.102.184515.
  * [72] S. Baier, M. J. Mark, D. Petter, K. Aikawa, L. Chomaz, Z. Cai, M. Baranov, P. Zoller and F. Ferlaino,  _Extended Bose-Hubbard models with ultracold magnetic atoms_ ,  Science 352(6282), 201 (2016),  10.1126/science.aac9812,  https://www.science.org/doi/pdf/10.1126/science.aac9812.
  * [73] V. W. Scarola and S. Das Sarma,  _Quantum phases of the extended Bose-Hubbard Hamiltonian: possibility of a supersolid state of cold atoms in optical lattices_ ,  Phys. Rev. Lett. 95, 033003 (2005),  10.1103/PhysRevLett.95.033003.
  * [74] J. Koepsell, J. Vijayan, P. Sompet, F. Grusdt, T. A. Hilker, E. Demler, G. Salomon, I. Bloch and C. Gross,  _Imaging magnetic polarons in the doped Fermi–Hubbard model_ ,  Nature 572(7769), 358 (2019),  10.1038/s41586-019-1463-1.
  * [75] J. Koepsell, D. Bourgund, P. Sompet, S. Hirthe, A. Bohrdt, Y. Wang, F. Grusdt, E. Demler, G. Salomon, C. Gross and I. Bloch,  _Microscopic evolution of doped Mott insulators from polaronic metal to Fermi liquid_ ,  Science 374(6563), 82 (2021),  10.1126/science.abe7165,  https://www.science.org/doi/pdf/10.1126/science.abe7165.
  * [76] J. Gavoret and P. Nozières,  _Structure of the perturbation expansion for the Bose liquid at zero temperature_ ,  Annals of Physics 28(3), 349 (1964),  https://doi.org/10.1016/0003-4916(64)90200-3.
  * [77] Y. A. Nepomnyashchii and A. A. Nepomnyashchii,  _Infrared divergence in field theory of a Bose system with a condensate_ ,  Sov. Phys. - JETP (Engl. Transl., United States) 48(3) (1978).
  * [78] D. S. Petrov and G. V. Shlyapnikov,  _Interatomic collisions in a tightly confined Bose gas_ ,  Phys. Rev. A 64, 012706 (2001),  10.1103/PhysRevA.64.012706.
  * [79] D. S. Petrov, M. Holzmann and G. V. Shlyapnikov,  _Bose-Einstein condensation in quasi-2D trapped gases_ ,  Phys. Rev. Lett. 84, 2551 (2000),  10.1103/PhysRevLett.84.2551.
  * [80] A. Rançon and N. Dupuis,  _Universal thermodynamics of a two-dimensional Bose gas_ ,  Phys. Rev. A 85, 063607 (2012),  10.1103/PhysRevA.85.063607.
  * [81] D. J. Scalapino, S. R. White and S. C. Zhang,  _Superfluid density and the Drude weight of the Hubbard model_ ,  Phys. Rev. Lett. 68, 2830 (1992),  10.1103/PhysRevLett.68.2830.
  * [82] D. J. Scalapino, S. R. White and S. Zhang,  _Insulator, metal, or superconductor: the criteria_ ,  Phys. Rev. B 47, 7995 (1993),  10.1103/PhysRevB.47.7995.


# Heuristic approach to trajectory correlation functions in bounded regions
with Lambert scattering walls

T. Rao  R. Golub  Physics Department, North Carolina State University,
Raleigh, North Carolina 27695, USA

###### Abstract

The behavior of spins undergoing Lamor precession in the presence of time
varying fields is of interest to many research fields. The frequency shifts
and relaxation resulting from these fields are related to their power spectrum
and can be determined from the Fourier Transform of the auto-correlation
functions of the time varying field. Using the method of images [C. M. Swank,
A. K. Petukhov, and R. Golub, Phys. Lett. A 376, 2319 (2012)] calculated the
position-position auto-correlation function for particles moving in a
rectangular cell with specular scattering walls. In this work we present a
heuristic model that extends this work to the case of Lambert scattering
walls. The results of this model are compared to simulation and show good
agreement from the ballistic to diffusive regime of gas collisions, for both
square and general rectangular cells. This model requires three parameters,
two of which describe the distribution of images in the case of a square cell,
and one of which describes the asymmetry in the mixing of the x and y
components of the velocity in the case of non-square rectangular cells.

## Introduction

Systems of spins moving under the influence of static and time varying
magnetic fields are of wide-ranging scientific interest. It was Bloembergen,
Purcell, and Pound [1] who first showed that the relaxation time of the spins
is given by the power spectrum of the fluctuating field evaluated at the Lamor
frequency. In general, the power spectra can be written in terms of the
Fourier transform of the auto-correlation function of the fluctuating field.
See the introductions of [2, 3] for a brief history of the field.

One of the applications of these techniques is in next generation neutron
electric dipole moment (nEDM) searches, which require measurements of spin
dynamics at the nanohertz level. Additionally the motional magnetic field
resulting from the interaction between the spins and the electric field in
such experiments results in “false EDM” systematic error that needs to be
accounted for [4]. Redfield [5], Slichter [6], and McGregor[7] have given
formal derivations of the relation between relaxation and the auto-correlation
functions of the fields, and the method was applied to the “false EDM”[4]
systematic error in nEDM searches [8]. General methods for calculating auto-
correlation functions for particles diffusing in inhomogenous fields can be
found in [9, 10].

Additionally, methods for calculating the trajectory auto-correlation
functions, or the frequency shift resulting from the motional magnetic field
have been developed for a variety of geometries, wall conditions, and
collision rates. Barabanov et al [11] used a calculation of the velocity-
velocity correlation function for indiviual ballistic trajectories to obtain
the motional field frequency shift in the case of a circular cell with
specular scattering walls. This result was expanded to include the effect of
gas collisions by noting the spectrum was a sum of hamonic oscillators, and
including a damping term to model the gas collisions. The work of Golub,
Steyerl et al [12, 13] similarly considered cylindrical cells in terms of
individual ballistic trajectories, but now with diffuse Lambert scattering
walls. Additionally, the case of square and more general rectangular cells for
specular and Lambert scattering walls were solved for. However, noticeable
qualitative differences between their calculation and simulation exist in the
case of non-square rectangular cells with Lambert scattering walls.
Specifically, simulation shows two peaks in the spectrum, while their
calculation only shows one. Clayton [14] also considers the case of
reactangular cells, but in the diffusive limit rather than the ballistic
limit, and uses the condtional probabality of finding the particle at a
position inside the cell to calculate the trajectory correlation functions.
This work covers the 1D, 2D, and 3D cells but only works when gas collisions
are frequent enough that the diffusion equation is valid. This work was
expanded upon by Swank et al [15, 16] giving solutions for 1D, 2D, and 3D
rectangular cells that apply from the ballistic to diffusive limit, but only
in the case of specular wall collisions.

The method in [15, 16] utilizes the method of images by applying the
conditional probability of finding a particle at position $r$, given initial
position $r_{0}$, and travel time $\tau$, for an unbounded domain [17]. As
explained in Fig. 1 going from the real cell to the image cells can be
represented by a triangle wave periodic extension
$\vec{\tilde{r}}_{0}(r_{0})$.

In this work we extend the technique of Swank et al to non specular wall
conditions, specifically the case of Lambert scattering walls. We alter the
position of the images in an attempt to alter the scattering angle at the
wall. This is accomplished by introducing a distribution of reflecting wall
positions. The properties of the image distribution required to reproduce
Lambert scattering is determined by comparing the resulting spectrum of the
trajectory correlation function to simulation. In this way we develop an
heuristic model of 2D rectangular cells with Lambert scattering walls that’s
valid from the ballistic to diffusive regime. Previously, Lambert scattering
walls have only been considered in the ballistic limit [12],

FIG. 1: The physical region (shaded) and a few of the periodic image cells for
the case of specular scattering walls described in [15]. The shown image
positions all correspond to the same source inside the physical cell.
Trajectories between an image and an observation point correspond to a
trajectory between the source point and the observation point that has
undergone the same number of wall collisions as image cell walls interesected
by the image trajectory _._ The walls of the physical cell are located at
$x=\pm L_{x}/2$ and $y=\pm L_{y}/2$. The field in the physical region, in this
case $B(x)=x$, is mirrored at the boundary between image cells so that images
see the same field as the source. This results in a periodic extension of the
position inside the cell $\tilde{x}(x)$, in this case a traingle wave, over
the infinite domain.

giving results that agree with simulation in the case of cylindrical and
square cells, but not in the case of rectangular cells.

## The Position-Position Trajectory Correlation Function for Lambert
Scattering Walls

From [15] Eq. (10) we see we can write the position-position trajectory
correlation function as

$\left\langle\overrightarrow{\widetilde{r}}_{0}(t)\cdot\overrightarrow{r}(t+\tau)\right\rangle=\frac{1}{L_{x}L_{y}}\int\int
d^{2}r_{0}d^{2}r\overrightarrow{\widetilde{r}}_{0}(r_{0})\cdot\overrightarrow{r}P(\vec{r}-\vec{r_{0}},\tau)$
(1)

with the infinite domain conditional probability $P(\vec{r}-\vec{r_{0}},\tau)$
from [17]. The method of images allows the use of the infinite domain
conditional probability. The conditional probability of a trajectory in a
finite domain is replaced by the infinite domain conditional probability
between the image position $r_{0}$, and the observation point inside the
physical cell $r$, for travel time $\tau$. The dot product is evaluated using
the periodic extension of the physical cell $\vec{\tilde{r}}_{0}(r_{0})$ for
the source position.

From Eq. (1) the 1D correlation function from [15, 16] for a 2D cell can be
written as

$\left\langle\widetilde{x}_{0}(t)x(t+\tau)\right\rangle=\frac{1}{L_{x}}\intop_{-\infty}^{\infty}\intop_{-L_{x}/2}^{L_{x}/2}dx_{0}dx\widetilde{x}_{0}(x_{0})xP(\vec{r}-\vec{r_{0}},\tau)$
(2)

To get the Spectrum we must take the Fourier transform of Eq. (1). From [17]
we know the Fourier-Laplace transform of the 2D infinite domain conditional
probability is given by

$P(\vec{q},s)=\frac{t_{c}}{\sqrt{(t_{c}vq)^{2}+(1+st_{c})^{2}}-1}$ (3)

where $t_{c}$ is the mean gas collisions time, and $v$ is the particle speed.
If we use an even extension of the conditional probability to negative $\tau$
the Fourier transform of $P(\vec{r}-\vec{r_{0}},\tau)$ is a real even function
equal to $2\mathrm{Re}[\intop_{0}^{\infty}P(r,t)e^{-i\omega t}dt]$. This is
simply twice the real part of the Laplace transform of $P(r,t)$ for
$s=i\omega$ giving us

$S_{xx}(\omega)=\frac{2}{L_{x}}\mathrm{Re}\left[\intop_{-\infty}^{\infty}\intop_{-L_{x}/2}^{L_{x}/2}dx_{0}dx\widetilde{x}_{0}(x_{0})xP(\vec{r}-\vec{r_{0}},i\omega)\right]$
(4)

where $S_{xx}(\omega)$ is the spectrum of 1D correlation function, and
$P(\vec{r}-\vec{r_{0}},s=i\omega)$ is the Laplace transform of
$P(\vec{r}-\vec{r_{0}},\tau)$. Additionally we can rewrite
$P(\vec{r}-\vec{r_{0}},i\omega)$ as $\int
d\vec{q}P(\vec{q},i\omega)\frac{e^{i\vec{q}\cdot(\vec{r}-\vec{r}_{0})}}{2\pi}$,
where $P(\vec{q},i\omega)$ is the Fourier-Laplace transform given in Eq. (3).
Since the conditional probability for an infinite domain has no angular
dependence we can choose $\vec{q}=q\hat{x}$ without loss of generality. This
in turn results in

$S_{xx}(\omega)=\frac{2}{2\pi L_{x}}$

$\mathrm{\times
Re}\left[\intop_{-\infty}^{\infty}dq\intop_{-L_{x}/2}^{L_{x}/2}dx\intop_{-\infty}^{\infty}dx_{0}P(q,i\omega)e^{iq(x-x_{0})}\widetilde{x}_{0}x\right]$
(5)

or equivalently

$S_{xx}(\omega)=\frac{2}{2\pi L_{x}}$

$\times\mathrm{Re}\left[\intop_{-\infty}^{\infty}dqP(q,i\omega)\intop_{-L_{x}/2}^{L_{x}/2}dxe^{iqx}x\intop_{-\infty}^{\infty}dx_{0}e^{iqx_{0}}\widetilde{x}_{0}\right]$
(6)

For specular walls, a trajectory’s wall reflection is accounted for by a
single image position, since the reflected angle is exactly determined by the
incoming angle. For rough walls the reflected angle is random and follows a
cosine distribution. Because of this a single image is no longer sufficient to
describe the scattering. Instead there exists a distribution of image
positions that correspond to the possible reflected trajectories (different
observation points), or equivalently, each image position can be thought of as
corresponding to multiple possible sources inside the physical cell. If the
specular image is altered by changing the size of the image cells, the
$\widetilde{x}_{0}$ corresponding to a given image is also altered. We attempt
to find a distribution of image cell sizes that generate the correct
distribution of $\widetilde{x}_{0}$ to reproduce Lambert wall scattering when
calculating the trajectory correlation function. The distribution of image
cell sizes that corresponds to Lambert scattering can then be used to
calculate an effective extension to the x axis that can be used in place of
the specular periodic extension.__

In [15, 16] the image positions for specular walls are accounted for by a
triangle wave extension of the particle’s position outside of the cell (See
Fig. 1)

$\widetilde{x}_{0}(x_{0})=\sum_{n=odd}\frac{-i^{n}2L_{x}}{\pi^{2}n^{2}}e^{i\frac{\pi
n}{L_{x}}x_{0}}$ (7)

So for Lambert scattering walls we need an appropriate ensemble average of
triangle wave extensions that will reproduce the wall scattering. This amounts
to taking an ensemble average of triangle waves with different periods $l$,
which would correspond to the length of a side of the physical cell in the
case of specular scattering (Eq. (8)). Additionally we average the x and y
dimensions separately since we want to include all combinations of $l_{y}$ and
$l_{x}$. Note that for a square cell the symmetry in the real cell suggests
the distribution of $l_{x}$ and $l_{y}$ should be the same. Consequently, in
the case of a square cell, although the individual image cells being averaged
over are rectangles, $l_{x}$ and $l_{y}$ are statistically independent, the
effective cell will still be a square.

$\widetilde{x}_{0}(x_{0})=\sum_{odd}\left\langle\frac{-i^{n}2l}{\pi^{2}n^{2}}e^{i\frac{\pi
n}{l}x_{0}}\right\rangle$ (8)

If we assume a Gaussian distribution in $l$ with average value $L_{modx}$, and
with width $\sigma$, Eq. (8) becomes

$\widetilde{x}_{0}(x_{0})=\sum_{odd}\frac{-i^{n}2}{\pi^{2}n^{2}}\intop_{0}^{\infty}le^{i\frac{n\pi}{l}x_{0}}\frac{e^{-\frac{(l-L_{modx})^{2}}{2\sigma^{2}}}}{\sqrt{2\pi}\sigma}dl$
(9)

Two examples comparing the x axis extension for Lambert scattering walls given
by Eq. (9) to the specular extension given by Eq. (7) are plotted in Fig. 2.

(a)  (b)

FIG. 2: Specular and rough wall extensions used to account for wall scattering
for a square and rectangular cell. (a) The specular and rough wall extensions
for a square cell with $L_{x}=2$, $L_{modx}=2.13$, and $\sigma=0.393$. Note
the region between -1 and 1 corresponding to the physical cell is barely
altered from the specular result. (b) The specular and rough wall extensions
for the long side of a rectangular cell with an aspect ratio of 4,
$L_{modx}=10.08$, and $\sigma=3.93$. The length of the sides are normalized to
1/2 of the length of the shorter side. Note the region between -4 and 4
corresponding to the physical cell begins to significantly deviate from the
specular case near the walls of the physical cell.

In order to calculate the spectrum of the correlation function, Eq. (6), we
need to take the Fourier transform of Eq. (9) giving

$\widetilde{x}_{0}(q)=\sum_{odd}\frac{-i^{n}2}{\pi^{2}n^{2}}\intop_{-\infty}^{\infty}\intop_{0}^{\infty}l\frac{e^{-\frac{(l-L_{modx})^{2}}{2\sigma^{2}}}}{\sqrt{2\pi}\sigma}e^{i\frac{n\pi}{l}x_{0}}e^{-iqx_{0}}dx_{0}dl$
(10)

evaluating the $x_{0}$ integral gives us

$\widetilde{x}_{0}(q)=\sum_{odd}\frac{-i^{n}2}{\pi^{2}n^{2}}\intop_{0}^{\infty}l\frac{e^{-\frac{(l-L_{modx})^{2}}{2\sigma^{2}}}}{\sqrt{2\pi}\sigma}2\pi\delta(\frac{n\pi}{l}-q)dl$
(11)

we then use the substitution $u=\frac{n\pi}{l}$ and evaluate the delta
function

$\widetilde{x}_{0}(q)=\sum_{odd}-i^{n}\frac{2\sqrt{2\pi}}{\sigma}\intop_{0}^{\infty}e^{-\frac{(\frac{n\pi}{u}-L_{modx})^{2}}{2\sigma^{2}}}\delta(u-q)\frac{1}{u^{3}}du\\\
=\sum_{odd}-i^{n}\frac{2\sqrt{2\pi}}{\sigma\left|q\right|^{3}}e^{-\frac{(\frac{n\pi}{q}-L_{modx})^{2}}{2\sigma^{2}}}$
(12)

The magnitude of $q$ used in the denominator of Eq. (12) is needed to account
for the case of negative $n$ and consequently negative $q$. We now define
$q_{x0}=n\pi/L_{modx}$ and rewrite this as

$\widetilde{x}_{0}(q)=\sum_{odd}-i^{n}\frac{4\pi}{L_{modx}q^{2}}\frac{n\pi}{\sqrt{2\pi}}\frac{e^{-\frac{n^{2}\pi^{2}(q_{{}_{x0}}-q)^{2}}{2(\sigma^{2}q^{2}q_{x0}^{2})}}}{(\left|q\right|\sigma
q_{x0})}$ (13)

From Eq. (13) we now have $\int dxe^{-iqx_{0}}\widetilde{x}_{0}(x_{0})$

Additionally we can evaluate the x integral in Eq. (6) to get

$\intop_{-L_{x}/2}^{L_{x}/2}dxe^{iqx}x=-\frac{i}{q^{2}}[qL_{x}\cos(\frac{qL_{x}}{2})-2\sin(\frac{qL_{x}}{2})]$
(14)

and we can now, using Eq. (13) and Eq. (14), write $S_{xx}$ , Eq. (6), as

$S_{xx}=\frac{4}{L_{x}L_{modx}}\mathrm{Re}\left[\sum_{odd}i^{n+1}\intop_{-\infty}^{\infty}dqP(q,i\omega)\frac{1}{q^{4}}[qL_{x}\cos(\frac{qL_{x}}{2})-2\sin(\frac{qL_{x}}{2})]\frac{n\pi}{\sqrt{2\pi}}\frac{e^{-\frac{n^{2}\pi^{2}(q_{{}_{x0}}-q)^{2}}{2(\sigma^{2}q^{2}q_{x0}^{2})}}}{(\left|q\right|\sigma
q_{x0})}\right]$ (15)

For a square $S_{xx}=S_{yy}$ and therefore the 2D spectrum $S_{rr}$ is simply
$2S_{xx}$. For a more general rectangle we could write

$S_{rr}=\frac{4}{L_{x}L_{modx}}\mathrm{Re}\left[\sum_{odd}i^{n+1}\intop_{-\infty}^{\infty}dqP(q,i\omega)\frac{1}{q^{4}}[qL_{x}\cos(\frac{qL_{x}}{2})-2\sin(\frac{qL_{x}}{2})]\frac{n\pi}{\sqrt{2\pi}}\frac{e^{-\frac{n^{2}\pi^{2}(q_{{}_{x0}}-q)^{2}}{2(\sigma_{x}^{2}q^{2}q_{x0}^{2})}}}{(\left|q\right|\sigma_{x}q_{x0})}\right]$

$+\frac{4}{L_{y}L_{mody}}\mathrm{Re}\left[\sum_{odd}i^{n+1}\intop_{-\infty}^{\infty}dqP(q,i\omega)\frac{1}{q^{4}}[qL_{y}\cos(\frac{qL_{y}}{2})-2\sin(\frac{qL_{y}}{2})]\frac{n\pi}{\sqrt{2\pi}}\frac{e^{-\frac{n^{2}\pi^{2}(q_{{}_{y0}}-q)^{2}}{2(\sigma_{y}^{2}q^{2}q_{y0}^{2})}}}{(\left|q\right|\sigma_{y}q_{y0})}\right]$
(16)

similarly to the case of specular walls [15]. However, unlike specular wall
collisions which preserve the magnitude of $v_{x}$ and $v_{y}$ individually,
Lambert wall collisions mix the x and y components only preserving the overall
magnitude of the velocity, in the case of elastic collisions considered here.
For a square cell the symmetry between x and y means on average the x and y
components aren’t effected by the velocity mixing. However, for rectangular
cells, $L_{y}>L_{x}$ and the velocity mixing is no longer symmetric. Each wall
collision results in change in the y velocity relative to the x velocity which
needs to be accounted for. Since the solution of a square cell is given by
$S_{xx}$ we can account for this effect by altering the velocity used in
calculating $S_{yy}$ relative to the nominal velocity in $S_{xx}$. We note the
velocity only appears in the conditional probability and assume a Gaussian in
effective speed for calculating $S_{yy}$ giving

$P(q,i\omega)\rightarrow\intop_{0}^{\infty}P(v,q,i\omega)\frac{e^{-\frac{(v-v_{nom})^{2}}{2\sigma_{v}^{2}}}}{\sqrt{2\pi}\sigma_{v}}dv$
(17)

Where the average speed $v_{nom}$ is determined from the ballistic time
$t_{b}$ for a square cell with sides of length $L_{x}$. The ballistic time is
taken to be the minimum time required to travel from the center of a square
cell to one of the walls. This then gives $v_{nom}=\frac{L_{x}}{2t_{b}}$. For
simplicity we will take $L_{x}=2$, and measure times in units of $t_{b}$ so
that $v_{nom}=1$.

However a Gaussian distribution cannot be exactly correct as the speed cannot
be less than zero, and a Gaussian allows for negative values if the standard
deviation is large enough. We will show that the width of the Gaussian is
$\approx 0.4$ in the ballistic limit for a 2x8 cell. Which results in $\approx
0.6$ % of the Gaussian’s area being less than zero. By changing the integral
so the distribution is cutoff for $v<0.1v_{nom}$ we have now have

$\intop_{0.1}^{\infty}P(v,q,i\omega)\frac{e^{-\frac{(v-1)^{2}}{2\sigma_{v}^{2}}}}{\sqrt{2\pi}\sigma}dv=P_{v}(q,i\omega,\sigma_{v})$
(18)

and problems at small $v$ due to the the unphysical Gaussian tail are removed
while only effecting $\approx 1.2$ % of the area of the probability
distribution for a 2x8 cell. Above $v$=0.1 the Gaussian approximation works
well.

Applying this change to the conditional probability we then have

$S_{rr}=\frac{4}{L_{x}L_{modx}}\mathrm{Re}\left[\sum_{odd}i^{n+1}\intop_{-\infty}^{\infty}dqP(q,i\omega)\frac{1}{q^{4}}[qL_{x}\cos(\frac{qL_{x}}{2})-2\sin(\frac{qL_{x}}{2})]\frac{n\pi}{\sqrt{2\pi}}\frac{e^{-\frac{n^{2}\pi^{2}(q_{{}_{x0}}-q)^{2}}{2(\sigma_{x}^{2}q^{2}q_{x0}^{2})}}}{(\left|q\right|\sigma_{x}q_{x0})}\right]$

$+\frac{4}{L_{y}L_{mody}}\mathrm{Re}\left[\sum_{odd}i^{n+1}\intop_{-\infty}^{\infty}dqP_{v}(q,i\omega,\sigma_{v})\frac{1}{q^{4}}[qL_{y}\cos(\frac{qL_{y}}{2})-2\sin(\frac{qL_{y}}{2})]\frac{n\pi}{\sqrt{2\pi}}\frac{e^{-\frac{n^{2}\pi^{2}(q_{{}_{y0}}-q)^{2}}{2(\sigma_{y}^{2}q^{2}q_{y0}^{2})}}}{(\left|q\right|\sigma_{y}q_{y0})}\right]$
(19)

Additionally we can recover the specular limit from Eq. (15) by taking
$\sigma_{x}$ goes to zero and $L_{modx}$ goes to $L_{x}$. Taking this limit
$\frac{n\pi}{\sqrt{2\pi}}\frac{e^{-\frac{n^{2}\pi^{2}(q_{{}_{0}}-q)^{2}}{2(\sigma^{2}q^{2}q_{0}^{2})}}}{(\left|q\right|\sigma
q_{0})}\rightarrow\delta(q_{x0}-q)$, giving

$S_{xx}=-\frac{8}{L_{x}^{2}}\sum_{odd}\mathrm{Re}\left[P(q_{x0},i\omega)\left(\frac{1}{q_{x0}}\right)^{4}i^{n+1}\sin(n\pi/2)\right]$

$=\frac{8}{L_{x}^{2}}\sum_{odd}\mathrm{Re}\left[P(q_{x0},i\omega)\right]\left(\frac{1}{q_{x0}}\right)^{4},$

the expected specular result. In the limit of large gas collision rate
($\lambda=t_{b}/t_{c}=1/t_{c}$) the wall collisions will be suppressed by the
gas collisions, and the Lambert scattering walls and specular walls results
should be equal. Therefore $\sigma$ in this limit should tend toward zero and
$L_{modx}$ to $L_{x}$. Additionally for rectangular cells in this limit we
should see that $\sigma_{v}$ goes to zero (as the velocity mixing from wall
collisions becomes negligible) so that $P_{v}\rightarrow P$.

## Model of Gas Collision Rate Dependence of the Parameters of the Gaussians

One way to consider the gas collision rate dependence of the correlation
function is that a trajectory that undergoes a gas collision but not a wall
collision should be accurately represented by the specular solution, as it
doesn’t interact with the walls. Therefore as the rate of gas collisions
increase we should expect to see more trajectories like this and the ensemble
average of triangle wave extensions from Eq. (8) will move toward a narrower
distribution closer to being centered on $L$. Therefore it’s reasonable to
think $L_{mod}$ and $\sigma$ should go like the fraction of wall collisions to
total collisions, which can be written as

$N_{w}/N=\frac{N_{w}}{N_{w}+N_{g}}=\frac{1}{1+N_{g}/N_{w}}$ (20)

where $N_{w}$ is the number of wall collisions, N is the total number of
collisions, and $N_{g}$ is the number of gas collisions. In order to
understand the gas collisions’ contribution let us consider the number of gas
collisions per wall collisions in a square cell with sides $L$. Describing the
gas collisions as a 2D random walk the number of gas collisions that occur
while traveling a distance $d$ is $N_{g}=(d/l_{mfp})^{2}$, where $l_{mfp}$ is
the mean free path. Given that $v=1$, and the gas collision rate
$\lambda=1/t_{c}$, then $l_{mfp}=vt_{c}=1/\lambda$. On average, a wall
collision will occur after traveling a distance $L/2$ in 1D. Therefore
adjusting for the the 4 walls (and therefore 4 wall collisions) of a square
cell, the number of gas collisions per wall collision in the 1D correlation
function for a 2D cell ($S_{xx}$ or $S_{yy}$) is

$\frac{N_{w}}{N}=\frac{1}{1+\frac{1}{4}\left(\frac{L\lambda}{2}\right)^{2}}$
(21)

where L is the length of cell in the specified direction

We can now calculate the expected $\sigma$ for different aspect ratios as a
function of $\lambda$ if we take the value of $\sigma$ in the ballistic limit
and scale it by $N_{w}/N$ giving

$\sigma(\lambda)=\sigma_{0}\frac{N_{w}}{N}$ (22)

Where $\sigma_{0}$ is the ballistic limit value of $\sigma$. Similarly, the
expected $L_{mod}(\lambda)$ would be

$L_{mod}(\lambda)=(L_{mod0}-L)\frac{N_{w}}{N}+L=\text{{\hbox{\Delta
L_{0}}}}\frac{N_{w}}{N}+L$ (23)

where $L_{mod0}$ is the ballistic limit value of $L_{mod}$, and $\Delta
L_{0}=L_{mod0}-L$.

However, when extending this model of the $\lambda$ dependence to work with
$\sigma_{v}$ we need to remember that the width of the $v$ Gaussian comes from
asymmetry in the mixing of the x and y components of the velocity, so that the
Gaussian gets narrower not simply as wall collisions are reduced but with an
additional factor of the relative number of gas collisions between the
calculation of $S_{yy}$ (a square with sides $L_{y}$) and $S_{xx}$ (a square
with sides $L_{x})$ for a given $\lambda$ giving

$\sigma_{v}(\lambda)=\sigma_{v0}\frac{1}{1+\frac{1}{4}a^{2}\left(\frac{aL_{x}\lambda}{2}\right)^{2}}$
(24)

Where $a=\frac{L_{y}}{L_{x}}$, and $\sigma_{v0}$ is the ballistic limit value
of $\sigma_{v}$.

## Model of Aspect Ratio Dependence in the Ballistic Limit

Although the previous model provides a way to find the spectrum of the
correlation function as a function of $\lambda$ given the ballistic limit
correlation function; $\sigma_{0}$, $L_{mod0}$, and $\sigma_{v0}$ still need
to be determined for each value of the aspect ratio ($a$). By understanding
how $\sigma_{0}$ , $L_{mod0}$, and $\sigma_{v0}$ scale with the aspect ratio
of the physical cell it is possible to write the spectrum of the correlation
function in a way that depends only on three parameters ($\sigma_{square}$,
$L_{modsquare}$, and $\sigma_{v,a=\infty}$).

###  $L_{mod0}$ scaling in terms of the square cell

For a square cell, scaling the cell to a larger square should preserve the
change in $L_{x}$ $\Delta L_{x}$ as a fraction of $L_{x}$ in order to keep the
same relative change in image positions from the specular case. In the case of
a rectangular cell where $L_{y}>L_{x}$, with aspect ratio integer
$a=\frac{L_{y}}{L_{x}}$, we would need $a$ $L_{x}\times L_{x}$ square cells
along the y direction, as shown in Fig. 3, to make our rectangle. Therefore we
might expect to see $\Delta L_{y}/L_{y}=a\Delta L_{x}/L_{x}$ for this more
general case. Assuming this is true we then have

$L_{mody0}=\frac{\Delta L_{y0}}{L_{y}}L_{y}+L_{y}=a\frac{\Delta
L_{x0}}{L_{x}}L_{y}+L_{y}\\\ =a^{2}\Delta L_{x0}+aL_{x}$ (25)

for a rectangular cell in the ballistic limit, where $\Delta
L_{x0}=L_{modsquare}-L$, and $L_{modsquare}$ is the average image cell wall
length for a square cell with sides $L_{x}$ in the ballistic limit.

FIG. 3: Diagram of a single 2x8 periodic cell constructed from 4 2x2 image
cells. The “real” boundary of the cell is red. The boundary of each 2x2 cell
alternates between black or yellow to show the overlap of the boundary between
neighboring cells. Note that 6 of the 4 2x2 squares’ walls used to make the
rectangular cell are internal to the rectangular cell.

From this same idea we would then also have that $\Delta L_{x}/L_{x}$ is
independent of the aspect ratio. Also note that in these arguments $a$ is
always a positive integer as we are constructing the rectangular cell in units
of the $L_{x}\times L_{x}$cell.

###  $\sigma_{y0}$ scaling in terms of the square cell

For Lambert scattering walls each wall collision mixes the x and y components
of a trajectory. Therefore wall collisions along the x (or y) axis also
contribute to $\sigma_{y}$ (or $\sigma_{x}$). So in determining $\sigma_{y0}$,
we start by considering square space $L_{y}\times L_{y}$ in size made from
$L_{x}\times L_{x}$ image cells in an attempt to include the collisions along
the x axis that contribute to $\sigma_{y0}$. Doing so we would then have

$\sigma_{y0}=a^{2}\sigma_{square}$ (26)

Where $\sigma_{square}$ is the standard deviation of the image cell wall
length for the 2x2 square cell in the ballistic limit. However, this
expression fails to take into account the fact that scaling the standard
deviation as if we had $a^{2}$ $L_{x}\times L_{x}$ square cells results in
effectively over estimating the contribution of wall collisions in the y
direction. This scaling adds $2(a-1)$ unphysical internal walls to each
physical rectangular cell in the space. See Fig. 3. If we assume every wall of
a $L_{x}\times L_{x}$ cell equally contributes $1/4\sigma_{square}$ to
$\sigma_{y0}$, then since there are $a$ rectangular physical cells over the
$L_{y}\times L_{y}$ space, we find we have overestimated $\sigma_{y}$ by
$1/2a(a-1)\sigma_{square}$. Therefore the correct $\sigma_{y}$ is

$\sigma_{y0}=a^{2}\sigma_{square}-\frac{1}{2}(a^{2}-a)\sigma_{square}\\\
=\frac{\sigma_{square}}{2}(a^{2}+a)$ (27)

###  $\sigma_{x0}$ scaling in terms of the square cell

Similarly, simply taking $\sigma_{x0}=\sigma_{square}$, will also overestimate
the contribution to $\sigma_{x0}$ made by wall collisions in the y direction.
For $\sigma_{y0}$ we saw the over estimation was
$\frac{1}{2}(a^{2}-a)\sigma_{square}$ over a $L_{y}\times L_{y}$ square space.
If we rescale this result for a $L_{x}\times L_{x}$ space, we find the
overestimate is $\frac{1}{2}\frac{(a^{2}-a)}{a^{2}}\sigma_{square}$.
$\sigma_{x0}$is then given by

$\sigma_{x0}=\sigma_{square}\left(1-\frac{\frac{1}{2}(a^{2}-a)}{a^{2}}\right)\\\
=\frac{\sigma_{square}}{2}\frac{(a^{2}+a)}{a^{2}}$ (28)

###  $\sigma_{v0}$ scaling as a measure of the asymmetry of the cell

As discussed earlier $\sigma_{v0}$ comes from the asymmetry of the cell. In
the case where $L_{y}\gg L_{x}$ the situation is effectively 1D and only the y
component of the velocity contributes to transport of particles across the
cell, but with a spread in the effective velocity in the y direction of the
particles determined by the number of collisions in the x direction.
Conversely, in the square cell x and y must be weighted equally so the
effective velocity in the y direction is always its nominal value of 1. Given
the value at $a=\infty$, $\sigma_{v0}$ could therefore reasonably be scaled by
a measure of the asymmetry of the cell giving

$\sigma_{v0}(a)=\sigma_{v,a=\infty}\frac{A_{yy}-A_{xy}}{A_{yy}}=\sigma_{v,a=\infty}\frac{a^{2}L_{x}^{2}-aL_{x}^{2}}{a^{2}L_{x}^{2}}\\\
=\sigma_{v,a=\infty}\frac{a^{2}-a}{a^{2}}$ (29)

where $A_{yy}$ is the area of a $L_{y}\times L_{y}$ square and $A_{xy}$ is the
area of a $L_{x}\times L_{y}$ rectangle.

## Comparison to simulation

Now that we have a model of the $\lambda$ and aspect ratio dependence, we
combine the Eqs. (22)-(24) with Eqs. (25), (27)-(29) to get

$\sigma_{y}=\sigma_{square}\frac{a^{2}+a}{2\left(1+\frac{1}{4}\left(\frac{aL_{x}\lambda}{2}\right)^{2}\right)}$
(30)

$\sigma_{x}=\sigma_{square}\frac{a^{2}+a}{2a^{2}\left(1+\frac{1}{4}\left(\frac{L_{x}\lambda}{2}\right)^{2}\right)}$
(31)

$L_{mody}=a^{2}\frac{L_{modsquare}-L_{x}}{1+\frac{1}{4}\left(\frac{aL_{x}\lambda}{2}\right)^{2}}+aL_{x}$
(32)

$L_{modx}=\frac{L_{modsquare}-L_{x}}{1+\frac{1}{4}\left(\frac{L_{x}\lambda}{2}\right)^{2}}+L_{x}$
(33)

$\sigma_{v}=\sigma_{v,a=\infty}\frac{a^{2}-a}{a^{2}}\frac{1}{1+\frac{1}{4}a^{2}\left(\frac{aL_{x}\lambda}{2}\right)^{2}}$
(34)

Leaving us with 3 free parameters $\sigma_{square}$, $L_{modsquare}$, and
$\sigma_{v,a=\infty}$. These parameters are determined by fitting to
simulation.

The simulation calculates trajectories for a uniform distribution of particles
in the cell. Over each timestep trajectories have constant velocity, with
checks for gas collisions at the end of each timestep. Gas collisions rotate
the trajectory to an uniformly distributed random angle, but do not change the
magnitude of the velocity. The rate of gas collisions is determined by an
exponential distribution with characteristic timescale $t_{c}$. The walls of
the 2D rectangular cell are Lambert scattering with the scattering angle at
each wall collision calculated by $\theta=2\arcsin(\sqrt{Rand()})-\pi/2$,
where $\theta$ is the scattering angle relative to the normal of the wall and
ranges from $-\pi/2$ to $\pi/2$, and $Rand()$ is a function that returns a
uniformly distributed random number between 0 and 1, thereby giving a cosine
distributed $\theta$. Wall collision times are determined by checking if
between the i and the i+1 timestep the trajectory will intersect with the
wall, and if so calculating a wall collision at the ith time. Small timesteps
relative to the ballistic time are needed to accurately calculate the
trajectories in the ballistic limit. The simulated trajectories are used to
calculate the position-position correlation function. The spectrum is then
determined by taking the fast Fourier transform. The simulation only considers
$\tau>0$, where the full spectrum includes negative tau, therefore the real
part of the simulation is equal to 1/2 of $S_{rr}$. Like in our model, times
are measured in units of $t_{b}$ and the speed of the particles is fixed to 1.
Additionally the center of the cell is located at the origin with walls at
$\pm L/2$.

Fitting our model to a 4000 particle simulation over 20000 0.0125$t_{b}$
timesteps of a 2x2 square cell in the ballistic limit ($\lambda=0.02$) with
Lambert scattering walls as shown in Fig. 4 we determine
$\sigma_{square}=0.393$ and $L_{modsquare}=2.13$. When fitting our model we
only use the $n=\pm 1$ terms, assuming the contribution from
$\left|n\right|>1$ terms is negligible. To find the value of
$\sigma_{v,a=\infty}$ that best fits the simulation, our model (using
$\sigma_{square}=0.393$ and $L_{modsquare}=2.13$) is fit to a 16000 particle
20000 timestep ballistic limit simulation for a 2x8 rectangular cell ($a=4$)
as shown in Fig. 5. From this fit we find $\sigma_{v,a=\infty}=0.54$. In all
plots comparing the model and simulation only 1/2 of the $S_{rr}$ calculated
from the model is plotted, and only $n=\pm 1$ terms of the calculation are
used.

FIG. 4: Simulation of the spectrum of the position-position auto-correlation
function for a 2x2 square with Lambert scattering walls in the ballistic limit
($\lambda=0.02$) along with a a fit to the data using the model described by
Eq. (19) and Eqs. (30)-(34). From this we find the best fit values for
$\sigma_{square}$ and $L_{modsquare}$ are 0.393 and 2.13 respectively. FIG. 5:
Simulation of the spectrum of the position-position auto-correlation function
for a 2x8 rectangle with Lambert scattering walls in the ballistic limit
($\lambda=0.02$) along with a a fit to the data using the model described by
Eq. (19) and Eqs. (30)-(34). From this we find the best fit value for
$\sigma_{v,a=\infty}=0.54$, given $\sigma_{square}=0.393$ and
$L_{modsquare}=2.13$.

Using this model and the parameter values found from out fitting the $a=1$ and
$a=4$ cells, we find good agreement in the ballistic limit for the entire
range of aspect ratios from $a=1$ to $a=6$ as shown in Fig. 6.

FIG. 6: Simulations of the spectrum of the position-position auto-correlation
function between $a=1$ and $a=6$ in the ballistic limit ($\lambda=0.02$) along
with the model described in by Eq. (19) and Eqs. (30)-(34) for the same aspect
ratios. Model curves are plotted using $\sigma_{v,a=\infty}=0.54$,
$\sigma_{square}=0.393$, and $L_{modsquare}=2.13$.

Similarly we see good agreement going from the ballistic to diffusive limit
for all aspect ratios, especially for the case of the square cell, with the
largest disagreement between model and simulation occurring around the
$S_{xx}$ peak at higher aspect ratios. Fig. 7 shows comparisons of our model
to simulation from the ballistic to diffusive limit for $a=1$ through 4.
Including the $\left|n\right|>1$ terms in the initial fitting may result in
slightly different values of the fitted parameters. For rectangular cells
where the spectrum is dominated by $S_{yy}$ this change should be most
noticeable in $\sigma_{v,a=\infty}$. Additionally since we only consider
$\left|n\right|=1$ the possibility of the fitted parameters having $n$
dependence is not explored.

(a)

(b)

(c)

(d)

FIG. 7: Simulations of the spectrum of the position-position auto-correlation
function from the ballistic to diffusive regime along with the model described
in by Eq. (19) and Eqs. (30)-(34) . Model curves are plotted using
$\sigma_{v,a=\infty}=0.54$, $\sigma_{square}=0.393$, and $L_{modsquare}=2.13$.
Curves for aspect ratios of 1, 2, 3, and 4 are plotted in subfigures (a), (b),
(c), and (d) respectively, with $\lambda\leq 1$ on the left and $\lambda>1$ on
the right. For the $\lambda\leq 1$ curves of subfigures (b),(c), and (d) we
additionally plot a zoomed in view of of the $\lambda=0.02$ curve to show the
peak resulting from $S_{xx}$ contribution to the correlation function. This
peak is suppressed at larger $\lambda$.

## Conclusion

In this work we have expanded on the results of [15, 16] and shown how the
method of images can be used to model Lambert scattering walls for rectangular
cells with reasonable accuracy from the ballistic to diffusive regime for
aspect ratios between 1 and 6. This is done by replacing the periodic
extension in the method of images used in [15, 16] by a Gaussian distribution
of extensions characterized by a standard deviation given by $\sigma_{square}$
and an average $L_{mod}$. For rectangular cells a Gaussian spread in the speed
is also needed in the $S_{yy}$ spectrum characterized by standard deviation
$\sigma_{v}$ in order to account for the effects of an asymmetric cell on the
x and y components of the velocity after a wall collision. It also should be
noted that the Gaussian distributions used in this model are only
approximations. The true distributions shouldn’t allow negative values of $l$
or $v$. The resulting correlation function is given by Eq. (19) and Eqs.
(30)-(34). Comparing to simulation we find $\sigma_{square}=0.393$,
$L_{modx}=2.13$, and $\sigma_{v,a=\infty}=0.54$. With only three free
parameters, two of which describe the wall scattering in a square cell, and a
third that describes asymmetry in the mixing of x and y components of the
velocity due to wall collisions, this models allows us to describe a
remarkably wide range of aspect ratios and gas collision rates, and reproduces
the feature of two peaks in the ballistic limit spectrum seen in simulations
of rectangular cells. The larger peak coming from the $S_{yy}$ contribution to
the spectrum, and the much smaller peak between $\omega=1.2$ and 1.4 coming
from $S_{xx}$. Although further work is needed to understand the differences
between model and simulation in the $S_{xx}$ peak shape. Given that this model
works by finding the distribution of triangle periodic extensions that
reproduces Lambert scattering, it may be possible to reproduce other wall
scattering conditions by finding the appropriate values of $\sigma_{square}$,
$L_{modx}$, and $\sigma_{v,a=\infty}$. As such this work has applications for
transport in mesoscopic systems with with a variety of boundary conditions
ranging from purely specular to the purely Lambert scattering case discussed
in this work. Additionally the trajectory correlation functions in this work
have applications to the problem of spin relaxation and frequency shifts for
polarized gasses in inhomogeneous magnetic fields. As well as calculations of
the motional field “false nEDM” frequency shift seen in nEDM experiments.

## References

  * [1] N. Bloembergen, E. Purcell, and R. Pound, Relaxation effects in nuclear magnetic resonance absorption, Phys. Rev. 73, 679 (1948)
  * [2] R. Golub, R. M. Rohm, and C. M. Swank, Reexamination of relaxation of spins due to a magnetic field gradient: Identity of the Redfield and Torrey theories, Phys. Rev. A 83, 023402 (2011)
  * [3] G. Pignol, M. Guigue, A. Petukhov, and R. Golub, Frequency shifts and relaxation rates for spin-1/2 particles moving in electromagnetic fields, Phys. Rev. A 92, 053407 (2015)
  * [4] J. M. Pendlebury, W. Heil, Y. Sobolev, P. G. Harris, J. D. Richardson, R. J. Baskin, D. D. Doyle, P. Geltenbort, K. Green, M. G. D. van der Grinten, P. S. Iaydjiev, S. N. Ivanov, D. J. R. May, and K. F. Smith, Geometric-phase-induced false electric dipole moment signals for particles in traps, Phys. Rev. A 70, 032102 (2004)
  * [5] A. Redfield, On the theory of relaxation processes, IBM J. Res. Dev. 1, 19 (1957)
  * [6] C. P. Slichter, Principles of Magnetic Resonance, 3rd ed. (Springer, Berlin-Heidelberg, 1996)
  * [7] D. D. McGregor, Transverse relaxation of spin-polarized 3 He gas due to a magnetic field gradient, Phys. Rev. A 41, 2631 (1990)
  * [8] S. K. Lamoreaux and R. Golub, Detailed discussion of a linear electric field frequency shift induced in confined gases by a magnetic field gradient: Implications for neutron electric-dipole-moment experiments, Phys. Rev. A 71, 032104 (2005)
  * [9] A. K. Petukhov, G. Pignol, D. Jullien, and K. H. Andersen, Polarized 3He as a Probe for Short-Range SpinDependent Interactions, Phys. Rev. Lett. 105, 170401 (2010)
  * [10] S. M. Clayton, Spin relaxation and linear-in-electric-field frequency shift in an arbitrary, time-independent magnetic field, J. Magn. Reson. 211, 89 (2011)
  * [11] A. L. Barabanov, R. Golub, and S. K. Lamoreaux, Electric dipole moment searches: Effect of linear electric field frequency shifts induced in confined gases, Phys. Rev. A 74, 052115 (2006)
  * [12] R. Golub, A. Steyerl, C. Kaufman, and G. Müller, Geometric phases in electric dipole searches with trapped spin-1/2 particles in general fields and measurement cells of arbitrary shape with smooth or rough walls, Phys. Rev. A 92, 06123 (2015)
  * [13] A. Steyerl, C. Kaufman, G. Müller, S. S. Malik, A. M. Desai, and R. Golub Phys. Rev. A 89, 052129 (2014)
  * [14] S. M. Clayton, Spin relaxation and linear-in-electric-field frequency shift in an arbitrary, time-independent magnetic field, J. Magn. Reson. 211, 89 (2011)
  * [15] C. M. Swank, A. K. Petukhov, and R. Golub, Correlation functions for restricted Brownian motion from the ballistic through to the diffusive regimes, Phys. Lett. A 376, 2319 (2012)
  * [16] C. M. Swank, A. K. Petukhov, and R. Golub, Random walks with thermalizing collisions in bounded regions: Physical applications valid from the ballistic to diffusive regimes Phys. Rev. A 93, 062703 (2016)
  * [17] J. Masoliver, J. M. Porra, and G. H. Weiss, Some two- and threedimensional persistent random walks, Phys. A (Amsterdam, Neth.) 193, 469 (1993)