Dissipative cooling of spin chains by a bath of dipolar particles

The geometry of our system is as follows. The two atomic species are trapped by the same lattice laser beams, but see a different potential depth due to their different electronic structures. We will only consider situations where the fermions are effectively arranged into disconnected 1D chains, mostly out of computational time limits, while the bosons retain 3D coherence (see section 4.3). The formalism is nevertheless generic, and dissipation imposed on 2D or 3D lattice structures would deserve further scrutiny. The quantization axis z is defined along the magnetic field. The spin chain axis z' can differ from z (see figure 2).

Zoom In Zoom Out Reset image size

Several restrictions are made for the sake of simplicity, that are not fundamental requirements. (i) We consider that the external magnetic field is strong enough to maintain the bath polarization: ${g}_{B}{\mu }_{{\rm{B}}}B\gt {k}_{B}T,J$. (ii) We restrict ourselves to pseudo-spin 1/2 chains, where only two neighboring Zeeman states of the fermions are ever populated, and assume that the Zeeman energy difference between these two states Δ can be made negligible, ${\rm{\Delta }}\ll J$, or irrelevant (see in section 4.2 considerations on how to impose this). As $J\ll U$ this implies ${\rm{\Delta }}\ll U$: all single-particle energies are smaller than U, such that with the proper chemical potential the occurrence of holes and doublons remains negligible.

We now describe the interaction between the dipolar bosons and the spin chain. In the fermionic Mott regime for large enough U/t, with one atom per site, we make the simplified assumption that ${k}_{B}T\ll U$ ensures negligible residual number fluctuations [46, 47, 49, 50] (perfect Fock state on each site). Under this approximation the density and spin operators for the fermions factorize, and we can label fermions by their site. The density distribution of the jth fermion ${n}_{F}^{j}(\overrightarrow{r})$ is constant, set by the ground band Wannier function on site j. For the bosons, as the BEC is energetically constrained to the polarized state, spin and density also factorize. This enables us to write the coupling between spin chain and condensate as:

$$\begin{eqnarray}&&{H}_{{\rm{int}}}=\displaystyle \sum _{j}\int {{\rm{d}}}^{3}{{r}{\rm{d}}}^{3}{r}^{{\prime} }{n}_{F}^{j}({\overrightarrow{r}}^{{\prime} }){\hat{n}}_{B}(\overrightarrow{r}){V}_{{dd}}^{j}(\overrightarrow{r},{\overrightarrow{r}}^{{\prime} }).\end{eqnarray} \tag{ 4 }$$

Here, ${\hat{n}}_{B}$ is the bath density operator in second quantization. ${V}_{{dd}}^{j}$ is the interaction V_dd of equation (2), acting on the fermionic spin ${\overrightarrow{{\rm{\Sigma }}}}_{j}$ of the fermion on the jth site under the influence of a boson of spin $\overrightarrow{S}$ at position $\overrightarrow{r}$. It includes the added constraint of fixed BEC magnetization: the S⁺ and S⁻terms are neglected, and the spin operator S^z is replaced by the spin length S of this species. We use Fourier transformation and write:

$$\begin{eqnarray}&&{V}_{{dd}}^{j}(\overrightarrow{r},{\overrightarrow{r}}^{{\prime} })=\int \displaystyle \frac{{{\rm{d}}}^{3}k}{{(2\pi )}^{3}}{V}_{{dd}}^{j}(\overrightarrow{k})\exp ({\rm{i}}\overrightarrow{k}\cdot (\overrightarrow{r}-{\overrightarrow{r}}^{{\prime} }))\end{eqnarray} \tag{ 5 }$$

with

$$\begin{eqnarray}&&{V}_{{dd}}^{j}(\overrightarrow{k})=\displaystyle \frac{{\mu }_{0}}{3}{g}_{B}{g}_{F}{\mu }_{{\rm{B}}}^{2}\left(3\displaystyle \frac{{k}_{z}}{{k}^{2}}({\overrightarrow{{\rm{\Sigma }}}}_{j}\cdot \overrightarrow{k})-{{\rm{\Sigma }}}_{j}^{z}\right)S.\end{eqnarray} \tag{ 6 }$$

As in dipolar relaxation, total intrinsic spin is not conserved. The magnetization-changing collisions are tied to the spin–orbit coupling inherent to dipolar interactions [51].

We now describe the evolution of the spin chain under the assumption that the BEC behaves as a bath. Proper justification of this will be given once the formalism is entirely laid out, in section 3.6. The processes that transfer energy from the spin chain to the BEC connect an initial spin state $| {i}_{{\rm{spin}}}\rangle$ to a final spin state $| {f}_{{\rm{spin}}}\rangle$ with energy difference ${E}_{i}-{E}_{f}={E}_{{if}}$. Starting from an initial bath state $| {i}_{{\rm{bath}}}\rangle$, the corresponding rates are:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{i\to f}=\displaystyle \frac{2\pi }{{\hslash }}\displaystyle \sum _{| {f}_{{\rm{bath}}}\rangle }| \langle {f}_{{\rm{spin}}};{f}_{{\rm{bath}}}| {H}_{{\rm{int}}}| {i}_{{\rm{spin}}};{i}_{{\rm{bath}}}\rangle {| }^{2}\,\delta ({E}_{{if}}+{E}_{{if}}^{{\rm{bath}}}).\end{eqnarray} \tag{ 7 }$$

This equation involves a summation over all possible final bath states $| {f}_{{\rm{bath}}}\rangle$, with bath energy difference ${E}_{{if}}^{{\rm{bath}}}$. We express all equations for a finite box of volume V with periodic boundary conditions, so that this summation is discrete. Regarding the bath, we start our derivation in the zero-temperature limit, only to generalize to finite T in the later section 3.5. Quantum depletion is neglected (see section 4.3). Thus, the initial bath state is for now described by a coherent state $| {\rm{BEC}}\rangle$ with average atom number N₀.

We now derive an explicit expression of the coupling element of equation (7). Using equations (4) and (5), we have:

$$\begin{eqnarray}\langle {f}_{{\rm{spin}}};{f}_{{\rm{bath}}}| {H}_{{\rm{int}}}| {i}_{{\rm{spin}}};{i}_{{\rm{bath}}}\rangle & = & \displaystyle \sum _{j}\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{3}{k}_{1}}{{(2\pi )}^{3}}\langle {f}_{{\rm{spin}}}| {V}_{{dd}}^{j}({\overrightarrow{k}}_{1})| {i}_{{\rm{spin}}}\rangle \\ & & \cdot \,{\displaystyle \int }_{V}{{\rm{d}}}^{3}{r}^{{\prime} }{{\rm{e}}}^{-{\rm{i}}{\overrightarrow{k}}_{1}\cdot {\overrightarrow{r}}^{{\prime} }}{n}_{F}^{j}({\overrightarrow{r}}^{{\prime} })\\ & & \cdot \,{\displaystyle \int }_{V}{{\rm{d}}}^{3}r\langle {f}_{{\rm{bath}}}| {{\rm{e}}}^{{\rm{i}}{\overrightarrow{k}}_{1}\cdot \overrightarrow{r}}{\hat{n}}_{B}(\overrightarrow{r})| {\rm{BEC}}\rangle .\end{eqnarray} \tag{ 8 }$$

There, the corresponding physical process is a transfer of momentum ${\overrightarrow{k}}_{1}$ from the fermions to the bosons. The lattice potential constrains the fermions to their unchanged Wannier states, whose finite momentum widths allow for this momentum transfer : ${n}_{{F}}^{j}({\overrightarrow{r}}^{{\prime} })=| {{w}}_{{F}}({\overrightarrow{r}}^{{\prime} }-{\overrightarrow{R}}_{j}){| }^{2}$, where ${{w}}_{{F}}({\overrightarrow{r}}^{{\prime} }-{\overrightarrow{R}}_{j})$ is the ground band Wannier function centered on the site coordinate ${\overrightarrow{R}}_{j}$. We calculate the bath matrix element using the discrete basis of plane waves in cubic volume V with periodic boundary conditions: ${\hat{n}}_{B}(\overrightarrow{r})={\sum }_{{\overrightarrow{k}}_{2}}{\sum }_{{\overrightarrow{k}}_{3}}\tfrac{{{\rm{e}}}^{{\rm{i}}({\overrightarrow{k}}_{3}-{\overrightarrow{k}}_{2})\cdot \overrightarrow{r}}}{V}{\psi }^{\dagger }({\overrightarrow{k}}_{2})\psi ({\overrightarrow{k}}_{3})$. Writing ${F}[f](\overrightarrow{k})$ the Fourier transform of a function $f(\overrightarrow{r})$, we obtain:

$$\begin{eqnarray}&&\langle {f}_{{\rm{spin}}};{f}_{{\rm{bath}}}| {H}_{{\rm{int}}}| {i}_{{\rm{spin}}};{i}_{{\rm{bath}}}\rangle =\displaystyle \int \displaystyle \frac{{{\rm{d}}}^{3}{k}_{1}}{{(2\pi )}^{3}}\displaystyle \sum _{j}{{\rm{e}}}^{-{\rm{i}}{\overrightarrow{k}}_{1}\cdot {\overrightarrow{R}}_{j}}\langle {f}_{{\rm{spin}}}| {V}_{{dd}}^{j}({\overrightarrow{k}}_{1})| {i}_{{\rm{spin}}}\rangle \\ &&\,\,\cdot \,{F}[{| {{w}}_{{F}}| }^{2}]({\overrightarrow{k}}_{1})\displaystyle \sum _{{\overrightarrow{k}}_{2}}\langle {f}_{{\rm{bath}}}| {\psi }^{\dagger }({\overrightarrow{k}}_{2})\psi ({\overrightarrow{k}}_{2}-{\overrightarrow{k}}_{1})| {\rm{BEC}}\rangle .\end{eqnarray} \tag{ 9 }$$

This expression explicitly relates the change in collective spin state to the interference of phonon radiation patterns from the different fermions, as an operator acting collectively on all the spins of the chain appears:

$$\begin{eqnarray}&&{V}_{{\rm{col}}}({\overrightarrow{k}}_{1})=\displaystyle \sum _{j}{{\rm{e}}}^{-{\rm{i}}{\overrightarrow{k}}_{1}\cdot {\overrightarrow{R}}_{j}}{V}_{{dd}}^{j}({\overrightarrow{k}}_{1})\end{eqnarray} \tag{ 10 }$$

(see figure 2). The operator V_col defines selection rules regarding which collective spin states are coupled. For these rules, we diagonalize the effective Hamiltonian H_mag together with the total pseudo-spin length operator ${\overrightarrow{{\rm{\Sigma }}}}_{{\rm{eff}},{\rm{tot}}}^{2}={({\sum }_{j}{\overrightarrow{{\rm{\Sigma }}}}_{{\rm{eff}},j})}^{2}$ and with both the total spin and pseudo-spin projection operators, e.g. ${{\rm{\Sigma }}}_{{\rm{tot}}}^{z}={\sum }_{j}{{\rm{\Sigma }}}_{j}^{z}$. In this eigenstate basis, if $| {i}_{{\rm{spin}}}\rangle$ and $| {f}_{{\rm{spin}}}\rangle$ differ by up to one in pseudo-spin length ${{\rm{\Sigma }}}_{{\rm{eff}},{\rm{tot}}}$ and 0 in projection ${{\rm{\Sigma }}}_{{\rm{tot}}}^{z}$

$$\begin{eqnarray}&&\langle {f}_{{\rm{spin}}}| {V}_{{\rm{col}}}(\overrightarrow{k}){| {i}_{{\rm{spin}}}\rangle }_{{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{tot}}}^{z}=0}\propto \left(3\displaystyle \frac{{k}_{z}^{2}}{{k}^{2}}-1\right)\cdot \displaystyle \sum _{j}{{\rm{e}}}^{-{\rm{i}}\overrightarrow{k}\cdot {\overrightarrow{R}}_{j}}\langle {f}_{{\rm{spin}}}| {{\rm{\Sigma }}}_{j}^{z}| {i}_{{\rm{spin}}}\rangle \end{eqnarray} \tag{ 11 }$$

and the corresponding two-body operator is Σ^z S^z. If $| {i}_{{\rm{spin}}}\rangle$ and $| {f}_{{\rm{spin}}}\rangle$ differ by up to one in pseudo-spin length Σ_{eff, tot} and ±1 in projection ${{\rm{\Sigma }}}_{{\rm{tot}}}^{z}$

$$\begin{eqnarray}&&\langle {f}_{{\rm{spin}}}| {V}_{{\rm{col}}}(\overrightarrow{k}){| {i}_{{\rm{spin}}}\rangle }_{{\rm{\Delta }}{{\rm{\Sigma }}}_{{\rm{tot}}}^{z}=\pm 1}\propto 3\displaystyle \frac{{k}_{z}}{{k}^{2}}\displaystyle \frac{{k}_{x}\mp {{\rm{i}}{k}}_{y}}{2}\cdot \displaystyle \sum _{{j}}{{\rm{e}}}^{-{\rm{i}}\overrightarrow{k}\cdot {\overrightarrow{R}}_{{j}}}\langle {f}_{{\rm{spin}}}| ({{\rm{\Sigma }}}_{{j}}^{x}\pm {\rm{i}}{{\rm{\Sigma }}}_{{j}}^{y})| {i}_{{\rm{spin}}}\rangle \end{eqnarray} \tag{ 12 }$$

and the corresponding two-body operator is Σ^± S^z. Otherwise, the coupling vanishes.

The coupling term (equation (9)) is sensitive to the potentials seen by the bosons. Here, we assume that the bosons in the lattice retain large coherence over the 3D sample, which we check independently (see section 4.3). We develop accordingly the BEC matrix element in equation (9) by transforming annihilation and creation operators from the plane wave basis $\psi (\overrightarrow{q})$ to the Bloch wave basis $\phi (\overrightarrow{q})$, and then to the Bogoliubov excitation basis $b(\overrightarrow{q})$ and BEC state ${\hat{c}}_{{\rm{BEC}}}$.

The first transformation, between plane waves and Bloch waves, is based on the decomposition of the Bloch band structure of our lattice, of axes (x', y', z'). We assume a cubic lattice with same period along all three axes, and reciprocal lattice wavevector length K. Then

$$\begin{eqnarray}&&\phi (\overrightarrow{q})=\displaystyle \sum _{{n}_{x}\in {\mathbb{Z}}}\displaystyle \sum _{{n}_{y}\in {\mathbb{Z}}}\displaystyle \sum _{{n}_{z}\in {\mathbb{Z}}}{a}_{{n}_{x}}^{* }[{q}_{{x}^{{\prime} }}]{a}_{{n}_{y}}^{* }[{q}_{{y}^{{\prime} }}]{a}_{{n}_{z}}^{* }[{q}_{{z}^{{\prime} }}]\,\psi (\overrightarrow{q}+\overrightarrow{n}K),\end{eqnarray} \tag{ 13 }$$

where $\overrightarrow{q}={q}_{{x}^{{\prime} }}{\overrightarrow{e}}_{{x}^{{\prime} }}+{q}_{{y}^{{\prime} }}{\overrightarrow{e}}_{{y}^{{\prime} }}+{q}_{{z}^{{\prime} }}{\overrightarrow{e}}_{{z}^{{\prime} }}$ is in the first Brillouin zone, and we define the notation $\overrightarrow{n}={n}_{x}{\overrightarrow{e}}_{{x}^{{\prime} }}+{n}_{y}{\overrightarrow{e}}_{{y}^{{\prime} }}+{n}_{z}{\overrightarrow{e}}_{{z}^{{\prime} }}$. The lattice can have different depths along the three axes, such that the ${a}_{{n}_{x}},{a}_{{n}_{y}},{a}_{{n}_{z}}$ implicitly differ. We omit in our notations the band indices: in this work only the ground band is ever populated given the relevant energies. We generalize the notation: ${a}_{\overrightarrow{n}}(\overrightarrow{q})={a}_{{n}_{x}}[{q}_{{x}^{{\prime} }}]{a}_{{n}_{y}}[{q}_{{y}^{{\prime} }}]{a}_{{n}_{z}}[{q}_{{z}^{{\prime} }}]$.

The second step is the Bogoliubov transformation [52]:

$$\begin{eqnarray}&&\phi (\overrightarrow{q})=u(\overrightarrow{q})b(\overrightarrow{q})+{v}^{* }(-\overrightarrow{q}){b}^{\dagger }(-\overrightarrow{q})+{\delta }^{K}(\overrightarrow{q},\overrightarrow{0}){\hat{c}}_{{\rm{BEC}}},\end{eqnarray} \tag{ 14 }$$

where δ^K is a Kronecker symbol. The prefactors $(u(\overrightarrow{q}),{v}^{* }(-\overrightarrow{q}))$ depend on kinetic and interaction terms (${E}_{\overrightarrow{q}},{V}_{\overrightarrow{q}}$, respectively) that can be strongly affected by the lattice potential:

$$\begin{eqnarray}&&u(\overrightarrow{q}),v(\overrightarrow{q})=\pm \sqrt{\displaystyle \frac{1}{2}\left(\pm 1+\displaystyle \frac{{E}_{\overrightarrow{q}}+{V}_{\overrightarrow{q}}}{\sqrt{{E}_{\overrightarrow{q}}({E}_{\overrightarrow{q}}+2{V}_{\overrightarrow{q}})}}\right)},\end{eqnarray} \tag{ 15 }$$

${E}_{\overrightarrow{q}}$ is the non-interacting energy for wavevector $\overrightarrow{q}$ of the Bloch wave. The interaction term ${V}_{\overrightarrow{q}}$ depends on wavevector due to the dipolar interaction within the polarized dipolar condensate. In the absence of a lattice potentials for bosons, ${V}_{\overrightarrow{q}}=g{\rho }_{0}+\tfrac{{\mu }_{0}}{3}{({g}_{B}{\mu }_{{\rm{B}}}S)}^{2}f(\overrightarrow{q}){\rho }_{0}$, with ρ₀ = N₀/V the mean bosonic density, g the 3D pseudopotential interaction parameter, and $f(\overrightarrow{q})=(3{q}_{z}^{2}/{q}^{2}-1)$. We generalize this in the presence of the lattice potential seen by the bosons, using a tight binding approximation: we introduce boson Wannier functions ${{{w}}_{{B}}}^{j}(\overrightarrow{r})={{w}}_{{B}}(\overrightarrow{r}-{\overrightarrow{R}}_{j})$ on each site j, and assume that ${{{w}}_{{B}}}^{j,* }(\overrightarrow{r}){{{w}}_{{B}}}^{k}(\overrightarrow{r})$ can be approximated by ${\delta }^{K}(j,k)| {{{w}}_{{B}}}^{j}(\overrightarrow{r}){| }^{2}$. We obtain:

$$\begin{eqnarray}&&{V}_{\overrightarrow{q}}=\tilde{g}{\rho }_{0}+\displaystyle \frac{{\mu }_{0}}{3}{({g}_{B}{\mu }_{{\rm{B}}}S)}^{2}{\rho }_{0}\displaystyle \sum _{\overrightarrow{n}\in {{\mathbb{Z}}}^{3}}{| {F}[{| {{w}}_{{B}}| }^{2}](\overrightarrow{q}+\overrightarrow{n}K)| }^{2}f(\overrightarrow{q}+\overrightarrow{n}K),\end{eqnarray} \tag{ 16 }$$

where $\tilde{g}$ is the bath contact interaction parameter renormalized by the lattice. For a 3D lattice confinement such that the bosons experience on-site interaction U_B, this interaction parameter reads $\tilde{g}={U}_{B}{(\lambda /2)}^{3}$, with λ/2 the lattice spacing. The Bogoliubov quasi-particles have an anisotropic dispersion relation $\epsilon (\overrightarrow{q})=\sqrt{{E}_{\overrightarrow{q}}({E}_{\overrightarrow{q}}+2{V}_{\overrightarrow{q}})}$.

We finally evaluate the coupling element equation (9), using the basis transformations of bosonic operators, equations (13) and (14). The coupling Hamiltonian H_int is of second order in boson field operators, which we just decomposed between Bogoliubov modes and BEC state mode. Consequently, the initial and final bath states can only be coupled if they differ by zero, one or two Bogoliubov annihilation and creation operators. The dominant terms for dissipation are those differing by one excitation, which also contain one BEC operator with macroscopic expectation value $\sqrt{{N}_{0}}$. Thus, for T = 0 and neglecting quantum depletion, the state $| {f}_{{\rm{bath}}}\rangle$ only differs from $| {\rm{BEC}}\rangle$ by the creation of a Bogoliubov excitation, parametrized by its wavevector $\overrightarrow{q}$:

$$\begin{eqnarray}&&| {f}_{{\rm{bath}}}(\overrightarrow{q})\rangle ={b}^{\dagger }(\overrightarrow{q})| {\rm{BEC}}\rangle .\end{eqnarray} \tag{ 17 }$$

We assume the BEC to occupy the zero-momentum Bloch wave state. Then, we get without explicit tight binding approximation

$$\begin{eqnarray}&&\langle {f}_{{\rm{spin}}};{f}_{{\rm{bath}}}(\overrightarrow{q})| {H}_{{\rm{int}}}| {i}_{{\rm{spin}}};{i}_{{\rm{bath}}}\rangle =\displaystyle \frac{\sqrt{{N}_{0}}}{V}\displaystyle \sum _{\overrightarrow{n}\in {{\mathbb{Z}}}^{3}}\displaystyle \sum _{\overrightarrow{m}\in {{\mathbb{Z}}}^{3}}\langle {f}_{{\rm{spin}}}| {V}_{{\rm{col}}}(\overrightarrow{q}+\overrightarrow{m}K)| {i}_{{\rm{spin}}}\rangle \\ &&\,\,\cdot \,({a}_{\overrightarrow{n}}^{* }[\overrightarrow{q}]{a}_{\overrightarrow{n}-\overrightarrow{m}}[0]{u}^{* }(\overrightarrow{q})+{a}_{\overrightarrow{n}}^{* }[0]{a}_{\overrightarrow{n}-\overrightarrow{m}}[-\overrightarrow{q}]{v}^{* }(\overrightarrow{q}))\,{F}[{| {{w}}_{{F}}| }^{2}](\overrightarrow{q}+\overrightarrow{m}K).\end{eqnarray} \tag{ 18 }$$

The tight binding approximation may also be applied along some dimensions to simplify equation (18). Applying it to all three dimensions leads to

$$\begin{eqnarray}&&\langle {f}_{{\rm{spin}}};{f}_{{\rm{bath}}}(\overrightarrow{q})| {H}_{{\rm{int}}}| {i}_{{\rm{spin}}};{i}_{{\rm{bath}}}\rangle =\displaystyle \frac{\sqrt{{N}_{0}}}{V}({u}^{* }(\overrightarrow{q})+{v}^{* }(\overrightarrow{q}))\\ &&\cdot \displaystyle \sum _{\overrightarrow{n}\in {{\mathbb{Z}}}^{3}}\langle {f}_{{\rm{spin}}}| {V}_{{\rm{col}}}(q+\overrightarrow{n}K)| {i}_{{\rm{spin}}}\rangle {F}[{| {{w}}_{{F}}| }^{2}](q+\overrightarrow{n}K)\cdot {F}[{| {{w}}_{B}| }^{2}](-\overrightarrow{q}-\overrightarrow{n}K).\end{eqnarray} \tag{ 19 }$$

From equation (19) we can see that the concentration of bosons onto Wannier functions by the lattice has an impact on the inter-species dipole coupling³ , despite the long-range nature of the dipolar interaction. Due to the anisotropic nature of V_col, the lattice can be beneficial or harmful to the coupling strength, as the reciprocal space summation ${\sum }_{\overrightarrow{n}}$ can involve terms of opposite signs, depending on the weights imposed by the Wannier functions.

At this stage we conclude our calculation for T = 0. To reveal explicitly the role of box volume, we replace the discrete summation on $| {f}_{{\rm{bath}}}(\overrightarrow{q})\rangle$ by a continuous one on $\overrightarrow{q}$.

$$\begin{eqnarray}{{\rm{\Gamma }}}_{i\to f}^{T=0} & = & \displaystyle \frac{2\pi }{{\hslash }}{\displaystyle \int }_{\mathrm{Brillouin}\mathrm{Zone}}{{\rm{d}}}^{3}q\displaystyle \frac{V}{{(2\pi )}^{3}}\delta (\epsilon (\overrightarrow{q})-{E}_{{if}})| \langle {f}_{{\rm{spin}}};{f}_{{\rm{bath}}}(\overrightarrow{q})| {H}_{{\rm{int}}}| {i}_{{\rm{spin}}};{i}_{{\rm{bath}}}\rangle {| }^{2}\\ & = & \,\displaystyle \frac{2\pi }{{\hslash }}{\displaystyle \int }_{0}^{\pi }\sin (\theta ){\rm{d}}\theta {\displaystyle \int }_{0}^{2\pi }{\rm{d}}\phi \,{q}^{2}(\theta ,\phi ,{E}_{{if}})\displaystyle \frac{V}{{(2\pi )}^{3}}\displaystyle \frac{1}{| \partial \epsilon (\overrightarrow{q})/\partial q{}_{| \theta ,\phi }| }\\ \,\, & & \cdot \,| \langle {f}_{{\rm{spin}}};{f}_{{\rm{bath}}}(\overrightarrow{q})| {H}_{{\rm{int}}}| {i}_{{\rm{spin}}};{i}_{{\rm{bath}}}\rangle {| }^{2},\end{eqnarray} \tag{ 20 }$$

where the 3D Dirac function results in an integration on the 2D anisotropic wavevector surface $\overrightarrow{q}(\theta ,\phi ,{E}_{{if}})$ defined by conservation of energy: $\epsilon (\overrightarrow{q})={E}_{{if}}$. In this paper, these energies will always lie in the ground band for the bosons.

At finite temperature, the initial bath state is not a pure state anymore. The theory above can then be generalized. We use a bath eigenstate basis $| \{{n}_{q}\}\rangle ={\prod }_{q}{({b}_{q}^{\dagger })}^{{n}_{q}}| {\rm{BEC}}\rangle$ consisting of Fock states of the Bogoliubov modes labeled q and of energy _q. The statistical probability of any such eigenstate $| {\{{n}_{q}\}}^{\mathrm{init}}\rangle$ is

$$\begin{eqnarray}&&{p}({\{{n}_{q}\}}^{\mathrm{init}})=\displaystyle \frac{{{\rm{e}}}^{-\displaystyle \sum _{q}{\epsilon }_{q}{n}_{q}^{\mathrm{init}}/{k}_{B}T}}{\displaystyle \sum _{\{{n}_{q}\}}{{\rm{e}}}^{-\displaystyle \sum _{q}{\epsilon }_{q}{n}_{q}/{k}_{B}T}}.\end{eqnarray} \tag{ 21 }$$

The rate between two spin chain states $| {i}_{{\rm{spin}}}\rangle ,| {f}_{{\rm{spin}}}\rangle$ with respective energies differing by E_i − E_f = E_if now writes:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{i\to f}=\displaystyle \sum _{{\{{n}_{q}\}}^{{\rm{init}}}}{p}({\{{n}_{q}\}}^{{\rm{init}}})\displaystyle \sum _{{\{{n}_{q}\}}^{{\rm{final}}}}\delta ({E}_{{if}}+{E}_{{if}}^{{\rm{bath}}}){| \langle {\{{n}_{q}\}}^{{\rm{final}}};{f}_{{\rm{spin}}}| {H}_{{\rm{int}}}| {\{{n}_{q}\}}^{{\rm{init}}};{i}_{{\rm{spin}}}\rangle | }^{2},\end{eqnarray} \tag{ 22 }$$

where ${E}_{{if}}^{{\rm{bath}}}={\sum }_{q}({n}_{q}^{\mathrm{init}}-{n}_{q}^{\mathrm{final}}){\epsilon }_{q}$. As explained in section 3.3, the dominant terms in the decomposition of H_int are those involving one Bogoliubov operator and the BEC mode.

If the spin energy diminishes (E_f < E_i), the most contributing $| {\{{n}_{q}\}}^{\mathrm{final}}\rangle$ differ from $| {\{{n}_{q}\}}^{\mathrm{init}}\rangle$ by one added excitation (phonon emission). Bose enhancement factors appear as excitation states are not necessarily initially empty: ${n}_{q}^{\mathrm{init}}\geqslant 0$. After statistical averaging over $| {\{{n}_{q}\}}^{\mathrm{init}}\rangle$, as explicit in equation (22), the rate ${{\rm{\Gamma }}}_{i\to f}$ is modified by finite temperature T through an additional factor:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{i\to f}^{T\gt 0}={{\rm{\Gamma }}}_{i\to f}^{T=0}\cdot \langle {n}_{q}+1\rangle ,\end{eqnarray} \tag{ 23 }$$

where $\langle {n}_{q}\rangle =\tfrac{1}{{{\rm{e}}}^{{E}_{{if}}/{k}_{B}T}-1}$ depends on the spin energy released in the phonon.

If the spin energy increases (E_f > E_i), the most contributing $| {\{{n}_{q}\}}^{\mathrm{final}}\rangle$ differ from $| {\{{n}_{q}\}}^{\mathrm{init}}\rangle$ by one less excitation. The rate ${{\rm{\Gamma }}}_{i\to f}$, associated to phonon absorption, is now non-zero. It scales with the inverse (phonon emission) process as:

$$\begin{eqnarray}&&{{\rm{\Gamma }}}_{i\to f}^{T\gt 0}={{\rm{\Gamma }}}_{f\to i}^{T=0}\cdot \langle {n}_{q}\rangle .\end{eqnarray} \tag{ 24 }$$

This asymetry between the processes pumping energy out (equation (23)) or in (equation (24)) the spin chain leads to thermalization between the spin chain collective degrees of freedom and the bath motional temperature. This is of particular interest for a sufficiently cold BEC, ${k}_{B}T\lesssim J$, such that $\langle {n}_{q}\rangle \lesssim 1$.

We conclude this section by defining properly the initial assumption of treating the BEC as a bath, (i) with which energy exchange processes are irreversible (justifying the Fermi golden rule), and (ii) such that its state is negligibly affected by interaction with the chain.

• (i)  
  For a given process $i\to f$, dephasing of the populated phonon modes is required to prevent coherent oscillations of the energy between chain and bath. As for spontaneous emission of light, this is satisfied even for an isolated BEC provided the number N_modes of excitation modes with energies close enough from E_if is much larger than one. We estimate N_modes by integrating the shell of wavevectors corresponding to excitations resonant with E_if, within an energy band width given by the coupling strength:

  $$\begin{eqnarray}&&{N}_{{\rm{modes}}}\simeq \displaystyle \frac{2\pi }{{\hslash }}{\int }_{\mathrm{Brillouin}\mathrm{Zone}}{{\rm{d}}}^{3}q\displaystyle \frac{V}{{(2\pi )}^{3}}\delta (\epsilon (\overrightarrow{q})-{E}_{{if}})| \langle {f}_{{\rm{spin}}};{f}_{{\rm{bath}}}(\overrightarrow{q})| {H}_{{\rm{int}}}| {i}_{{\rm{spin}}};{i}_{{\rm{bath}}}\rangle | .\end{eqnarray} \tag{ 25 }$$

  This amounts to a product of density of state and coupling element, taking into account the anisotropy of the coupling. For fixed BEC mean density N₀/V, matrix elements of H_int scale as $1/\sqrt{V}$, such that N_modes scales as $\sqrt{V}$ and incoherent dynamics is always justified at the thermodynamic limit. Once in the incoherent regime, the rate is independent on volume. For small clouds that would not satisfy this density of state criterion, irreversibility requires dissipation or damping of the phonons by other means, either inherent to finite temperature BECs [53–56], or external such as evaporation. The latter is experimentally relevant, as BEC temperatures lower than J are typically obtained in traps with low depths.
• (ii)  
  Furthermore, throughout the numerical results, the bath excitations will be assumed to be thermally distributed at a fixed temperature. Precisely, we assume that the emissions and absorptions of phonons affect negligibly the populations $\langle {n}_{q}\rangle$. For a strongly damped bath, this may be ensured by a sufficient heat capacity. For weakly damped baths as most BECs, this will be ensured either by evaporation, or by N_modes being larger than the number of phonons emitted within the same energy band. The latter criterion improves with the ratio of boson atom number to fermion atom number.

Using one of the most magnetic dipolar species (Cr, Er, Dy) as the bath, dipolar interactions will drive dissipation even on spin chains from the weakly dipolar alkali species. More exotic situations with actual dipole–dipole interactions within the spin chain could be studied, e.g. using an isotope of the strongly dipolar species for the chain as well.

In the prospect of studies of the t–J model, constituting spin chains from fermionic species in the Mott regime is very attractive. Lithium or Potassium are the most prominent in recent experiments, together with alkaline-earth species that however have much too weak Landé factors in the ground state for our present purpose. We will here discuss ⁴⁰K in its ground hyperfine state, with spin F = 9/2, and Landé factor g_F = 2/9. Fermionic Lithium would be more problematic due to a combination of complications⁴ .

A key point is the optical lattice potential seen by the bath species. It generally differs from that seen by the spin chain species. As mentioned earlier (section 3.3), a lattice of opposite sign may often lead to poor inter-species dipolar coupling, concentrating the bosons half-a-space away from the fermions. For Dy or Er in association with K spins, wavelengths > 767 nm would generate lattices of the same sign for both species. However larger J and large inter-species coupling can be obtained at shorter visible wavelengths, using the red-detuned vicinity of 'narrow' lines that would dominate over the effect of the wide lines at 421 and 401 nm, respectively, for Dy and Er. We suggest for example one narrow line of Dy at 626.1 nm (140 kHz) and one narrow line of Er at 582.8 nm (190 kHz) (see figure 3). The intervals 622.5–624.3 nm for Dy and 578.2–580.2 nm for Er are attractive, as in these intervals the boson lattice depths can be tuned from 0 to ∼0.4 times the fermion lattice depth while light scattering remains dominated by the fermions. Here we will use ¹⁶⁴Dy, which has the strongest magnetic dipole moment μ = g_B S μ_B = 10μ_B, but the shorter lattice wavelengths suggested for Er translates in higher super-exchange and fairly similar spin cooling rates despite the lower moment 7μ_B.

Zoom In Zoom Out Reset image size

In the following, we will consider ⁴⁰K in the presence of ¹⁶⁴Dy [58], in an optical lattice with wavelength λ_L ≈ 623 nm. We effectively create a 2D array of 1D spin chains by using a lattice, as seen by K, of depth 25 E_r^K along two axes x', y' (where ${E}_{r}^{K}={h}^{2}/(2{m}_{K}{\lambda }_{L}^{2})$ is the potassium recoil energy), and $3.5\,{E}_{r}^{K}$ along one axis z'. The spontaneous light emission rate is ∼0.2/s. The on-site interaction energy $U=7.5\,{t}_{{z}^{{\prime} }}$, where ${t}_{{z}^{{\prime} }}$ is the tunneling energy along the weak axis z'. Fermions can then be prepared in the 1D Mott regime with unit filling [45, 59, 60]. The super-exchange energy is, along the weak axis, J = h × 632 Hz, while it is along the transverse axes h ×  0.07 Hz, thus effectively creating decoupled 1D spin chains. The lattice depth for the bosons, on the other hand, sensitively depends on the lattice wavelength within a few nanometers, and will be used as an independent optimization variable in section 4.3.

The actual implementation requires to respect a certain hierarchy of energy scales, involving the magnetic field B, the BEC temperature T, and the super-exchange energy of fermions J (see figure 1). First, the bath temperature $T\lt J/{k}_{B}\simeq 30$ nK⁵ is such that low energy spin chain states are favored by coupling to the bath. Second, (J, k_BT) < g_B μ_B B ensures that the bath remains polarized. This is largely satisfied with fields in the mG range. Although the scheme might work also for a BEC with free magnetization, this study is beyond the scope of our calculations. Third, to study antiferromagnetic models the fermion Zeeman shift Δ should not favor trivially polarized spin chains. There are several ways to ensure this, with different consequences.

• (i)  
  The fermionic Zeeman states of interest may be kept quasi-degenerate by ${g}_{F}{\mu }_{{\rm{B}}}B\ll J$. Combined with our prior assumption of a polarized BEC, this implies ${g}_{F}\ll {g}_{B}$, which is not well satisfied for most of the atomic species combinations. To circumvent this, we propose to add micro-wave dressing to the spin states [64]: the linear Zeeman effect can be combined with an engineered quadratic shift $\alpha {m}_{F}^{2}$, such that the detuning between the two lowest Zeeman states is ${\rm{\Delta }}=| {E}_{{m}_{F}=-F}-{E}_{{m}_{F}=-F+1}| =| {g}_{F}{\mu }_{{\rm{B}}}B-\alpha ({F}^{2}-{(F-1)}^{2})| \ll J$. To emulate pseudo-spin 1/2 dynamics in a larger spin species, one further requires $| {E}_{{m}_{F}=-F+1}-{E}_{{m}_{F}=-F+2}| \gt J,{k}_{B}T$ to prevent population of the other Zeeman states. This implies $B\gt \tfrac{J}{{g}_{F}{\mu }_{{\rm{B}}}}\tfrac{2F-1}{2}\simeq 8$ mG. This protocol, and the effect of a residual Δ, are discussed further in appendix A. The calculations below assume such a setting. Varying Δ would enable studies of the Heisenberg antiferromagnetic model with free magnetization for varying bias field, complementing the more usual fixed magnetization setting in experiments.
• (ii)  
  Another possibility is to ensure that the energy released by magnetization-changing collisions, Δ, lies in the energy gaps within the boson band structure. Dipolar relaxation is then prevented by energy conservation, as seen in [51, 57]. The thermalization of the spin chain would then occur at fixed magnetization (i.e. fixed spin projection along $\overrightarrow{B}$). This approach requires a much lower degree of field control, and in particular releases the constraint that ${\rm{\Delta }}\ll J$. Cooling of the chain is slightly less efficient, as the magnetization-changing processes are suppressed. Nevertheless, processes given by equation (11) are still at play and may thermalize the collective spin at fixed magnetization. The magnetization is set initially by creating the appropriate spin imbalance using standard manipulation techniques.

The lattice potential experienced by the bosons has a crucial effect on phonon emission, strongly tied to the anisotropic character of the dipolar interaction. We first illustrate the importance of this anisotropy without the lattice potential. Figure 4 presents phonon emission diagrams for different orientations of a two-fermions chain in the triplet states. We make two observations on the emission directionality. (a) It is strongly spin dependent, as can be inferred by equations (11) and (12). (b) It is suppressed orthogonally to the spin chain, as the phase terms ${{\rm{e}}}^{-{\rm{i}}\overrightarrow{q}\cdot {\overrightarrow{R}}_{j}}$ in V_col are close to 1. Indeed, V_col then produces an identical rotation of all spins in the chain, which does not modify the collective spin and therefore is unable to drive any dissipative dynamics⁶ .

Zoom In Zoom Out Reset image size

The bath lattice potential introduces the following effects. (i) It can strongly influence the anisotropy of the phonon dispersion relations, and thus of the resonant phonon wavevector surface (see figure 5). This affects both the density of states and the phase terms ${{\rm{e}}}^{-{\rm{i}}\overrightarrow{q}\cdot {\overrightarrow{R}}_{j}}$ in V_col discussed above. (ii) It can generate either constructive or destructive interference through ${\sum }_{\overrightarrow{n}}[..]{V}_{{\rm{col}}}(\overrightarrow{q}+\overrightarrow{n}K)$ terms of equations (18) and (19), as the dipolar interaction changes sign depending on the orientation of $\overrightarrow{q}+\overrightarrow{n}K$.

Zoom In Zoom Out Reset image size

The lattice potential depth for bosons is thus an important control tool, with two main functions: (i) tuning the boson dispersion relation, to explore high-density of states, and also to stabilize against the bath dynamical instabilities of strongly magnetic species [67, 68] (for Dy, ${\epsilon }_{{dd}}={\mu }_{0}{\mu }^{2}m/12\pi {{\hslash }}^{2}{a}_{{sc}}\simeq 1.4$); (ii) maximizing the constructive interference and thus enhancing the coupling. Optionally, as mentioned in section 4.2, band gaps may also be used to enforce fixed magnetization dynamics with non-zero fermion Zeeman shifts.

In systematic studies of the cooling rates as a function of system parameters, for ${\theta }_{z-{z}^{{\prime} }}\ne 0$ the rates typically were more sensitive than for ${\theta }_{z-{z}^{{\prime} }}=0$ to the other parameters, including to the released energy E_if. Looking for conclusions that are robust with chain length, we will focus on ${\theta }_{z-{z}^{{\prime} }}=0$. We find that in such a case where the spin chain is oriented along the magnetic field axis, the lattice transverse to $\overrightarrow{B}$ mostly contributes to enhancing magnetization-conserving couplings, while the longitudinal lattice mostly contributes in stabilizing the dipolar gas, but quickly becomes detrimental to the couplings. Deeper transverse lattices are preferable until one has to consider the deconfinement transition from 1D tubes or 1D Mott states to 3D BEC [69–72], and the quantum depletion that arises before it [73].

Fully independent optimization of longitudinal (z') and transverse (x', y') boson lattice depths is possible using two laser sources, given the strong wavelength dependence of the boson light shift⁷ . We set our parameters to k_BT = 0.3 J, a mean density of ${\rho }_{0}=3\times {10}^{19}$/m³, i.e. a filling of 0.9 bosonic atoms per 3D site, in a transverse lattice depth $12{E}_{r}^{{\rm{Dy}}}$, and a longitudinal lattice depth $3.5{E}_{r}^{{\rm{Dy}}}$ (stabilization against dynamical instabilities operates from $\approx 2.5{E}_{r}^{{\rm{Dy}}}$). The first gap of the band structure is h × 5.2 kHz wide.

Based on 1D physics considerations neglecting the longitudinal lattice, the 1D Lieb–Liniger parameter $\gamma ={{mg}}_{1d}/{{\hslash }}^{2}{n}_{1d}\approx 1.3$, while the transverse tunneling is ${t}_{{\rm{perp}}}\approx 6\times {10}^{-2}{\mu }_{1d}$, which according to [71] and consistently with the experiment [72] would ensure being deep in the 3D coherent regime. From the 3D lattice perspective, we have $U/(2{t}_{{x}^{{\prime} }}+2{t}_{{y}^{{\prime} }}+2{t}_{{z}^{{\prime} }})\simeq 0.9$, far from the critical value ≃5.8 where the Mott transition could be reached. Following [73], we estimate quantum depletion based on Bogoliubov theory in the homogeneous system limit to be ${N}_{{qd}}/{N}_{0}=1/{N}_{0}\int {{\rm{d}}}^{3}{qV}/{(2\pi )}^{3}v{(\overrightarrow{q})}^{2}\simeq 5 \%$. The transverse lattice depth is thus kept conservatively low to ensure the BEC 3D coherence.

In this setting, the surface of resonant phonon wavevectors is reduced to two sheets, orthogonal to z' (figure 5). Due to the diverging density of states at the edges of the Brillouin zone, the phonon radiation pattern is heavily affected, as visible in figure 6. As the lattice is more confining in the (x', y') directions, the terms of ${\sum }_{\overrightarrow{n}}[...]{F}[{| {{w}}_{{B}}| }^{2}](-\overrightarrow{q}-\overrightarrow{n}K){V}_{{\rm{col}}}(\overrightarrow{q}+\overrightarrow{n}K)$ in equation (19) bear most weight for $\overrightarrow{n}$ in the (x', y') plane, in which V_col has constant sign. Thus, the summation is constructive.

Zoom In Zoom Out Reset image size

Finally, we estimate with these settings that for the representative process $| {T}_{0}\rangle \to | S\rangle$ in a length-2 chain, the number of contributing phonon modes is, for V = (20 μm)³, N_modes ≃ 1.7, and for V = (100 μm)³, N_modes ≃ 19. The incoherent description is thus justified without the need to consider evaporation for large but realistic BECs.

We now present the dissipative evolution of spin chains with the above parameters. To get a hint at how our simple picture scales towards the large system limit, we will repeat our calculation for varying 1D chain lengths L. Although our theory is valid for arbitrary chain length, we are limited by the computation time of the radiation pattern and anisotropic density of states, which have to be evaluated for all transition energies between spin chain eigenstates. The evolutions shown below rely on a three step approach. (i) Find all spin chain eigenstates, in a basis that also diagonalizes total pseudo-spin length and total spin projection along the polarizing field. (ii) Compute the rates between all of these states—this being the computationally limiting part as L is increased. We use for efficiency a tight binding approximation in the H_int matrix element along the two strong axes (x', y') (see appendix B), and along all axes for the ${V}_{\overrightarrow{q}}$ and thus the Bogoliubov description of the bath. (iii) Given all rates, simulate the dissipative evolution from a statistical distribution of the eigenstate probabilities.

First, we show in figure 7 the evolution with time of the probabilities of the four eigenstates from a length-2 spin chain, for various angles ${\theta }_{z-{z}^{{\prime} }}$ of the chain axis with respect to $\overrightarrow{B}$. The cooling dynamics shows very strong sensitivity with the alignment, as a consequence of the anisotropy of the dipolar interaction. For parallel alignment (${\theta }_{z-{z}^{{\prime} }}=0$), the fixed magnetization transition $| {T}_{0}\rangle \to | S\rangle$ has the fastest rate, with timescale ≈0.9 s. As will be demonstrated later, preparing the initial chain at total magnetization close to 0 (i.e., starting from a balanced spin mixture in the chain) may then provide the fastest dissipative route to low energy antiferromagnetic chains.

Zoom In Zoom Out Reset image size

We demonstrate in figure 8 that the spin chain reaches thermal equilibrium with the bath. The plot shows the occupation probabilities of the various spin chain eigenstates after several evolution times, for a chain of length L = 7. We start the evolution from a state with minimal magnetization, i.e. all states with total magnetization ± 1/2 have same probability, while all the other states have zero initial probability. This corresponds to a state of infinite spin temperature with the constraint of ± 1/2 magnetization. We observe that the evolution converges towards a distribution of probabilities exponential with respect to energy, i.e., a thermal distribution in the many-body states at free magnetization. We furthermore apply an exponential fit to the equilibrated distribution and extract a spin temperature of k_BT_S = 0.3 J which corresponds to the bath temperature. Two seconds of evolution already lead close to a thermal state at spin temperature k_BT_S ≃0.45 J.

Zoom In Zoom Out Reset image size

To extend these observations to longer chains, we define global, 'thermodynamic' quantities describing the spin chain: the magnetic energy E, rescaled by its initial value, and the global von Neumann spin entropy ${ \mathcal S }=-{\sum }_{i}{p}_{i}\mathrm{log}({p}_{i})$, where p_i are the probabilities of the collective eigenstates. In figure 9, we plot these quantities as a function of time, for chain lengths L ranging from 2 to 7. We actually compare the evolution from samples prepared at minimal initial magnetization, 0 or ± 1/2 depending on the chain length parity (left column), and prepared with free initial magnetization (right column). The dissipative evolution itself frees the magnetization, which is responsible for the initial entropy increase for samples prepared at fixed magnetization. Figures 9(a), (b) and (d), (e) show that the evolutions of energy and entropy for increasingly large chains tend to converge with L, for both preparation protocols. The equilibration rates ${{\rm{\Gamma }}}_{{\rm{equil}}}={\rm{d}}/{\rm{d}}{t}\mathrm{ln}(E-E(t\to \infty ))$ for both preparation protocols (figures 9(c), (f)), evaluated at early times, seem to quickly converge to similar and stable (size-independent) values. These observations suggest that long, mesoscopic chains may evolve with similar equilibration rates.

Zoom In Zoom Out Reset image size

We observe in these plots, already after 2 s, entropies per particle and spin temperatures that compare very favorably with the most recent experiments, e.g. [16, 60]. Ultimately, provided the equilibration rate is sufficient to overcome imperfections such as light scattering, a key point is that the outcome is tied to the attainable BEC temperature; and indeed, as mentioned in section 4.2, extremely cold BECs can be produced experimentally, in bulk and in optical lattices [61–63]. It thus decouples from the usual challenge to produce low entropy spinful fermions in lattices, tackled e.g. by [74, 75].

To conclude, we show in figure 10 spin correlations for the largest spin chain (L = 7) before and after the evolution, and compare them to correlations at zero-temperature. We use the following correlator definition, normalized to one for full correlation:

$$\begin{eqnarray}&&C(i,j)=4\langle {{\rm{\Sigma }}}_{{\rm{eff}},i}^{z}{{\rm{\Sigma }}}_{{\rm{eff}},j}^{z}\rangle .\end{eqnarray} \tag{ 26 }$$

The flat negative correlation of the initial state is a finite size bias from the constrained initial magnetization ±1/2. For example, if the first site has magnetization −1/2, the L −1 other sites have a positive average magnetization per site $\propto 1/(L-1)$ that reflects upon the correlator. At the final temperature of 0.3 J we observe alternating spin correlations already similar to the ground state, spanning over the whole sample. This correlator is accessible on experiments, especially using quantum gas microscope techniques discriminating two spin states, as in [76, 77].

Zoom In Zoom Out Reset image size

The Bloch states and their energies are computed by diagonalizing the non-interacting motion of the bosons and of the fermions in a truncated momentum space, ranging from −30.5K to +30.5K, along each of the three directions of space. The Wannier functions constructed from this diagonalization are by default truncated in Fourier space to a support of [ −4K, 4K] along each dimension.

We use twice in this paper a tight binding approximation, introducing the bosons Wannier functions and assuming ${{{w}}_{{B}}}^{j,* }(\overrightarrow{r}){{{w}}_{{B}}}^{k}(\overrightarrow{r})\equiv {\delta }^{K}(j,k)| {{{w}}_{{B}}}^{j}(\overrightarrow{r}){| }^{2}$. When lowering the bath lattice depth along some direction of space, this becomes improper, such that the bath interaction with the spin chain, equation (18), may be poorly approximated by equation (19). In the conditions of section 4, we find that for a length-2 chain, the $| {T}_{0}\rangle \to | S\rangle$ transition rate differs between the two expressions by 12%. The agreement is improved to 1% if we apply the tight binding approximation only along the transverse axes (x', y'):

$$\begin{eqnarray}&&\langle {f}_{{\rm{spin}}};{f}_{{\rm{bath}}}(\overrightarrow{q})| {H}_{{\rm{int}}}| {i}_{{\rm{spin}}};{i}_{{\rm{bath}}}\rangle =\displaystyle \frac{\sqrt{{N}_{0}}}{V}\displaystyle \sum _{{n}_{z}\in {\mathbb{Z}}}\displaystyle \sum _{\overrightarrow{m}\in {{\mathbb{Z}}}^{3}}\langle {f}_{{\rm{spin}}}| {V}_{{\rm{col}}}(\overrightarrow{q}+\overrightarrow{m}K)| {i}_{{\rm{spin}}}\rangle \\ &&\,\,({a}_{{n}_{z}}^{* }[{q}_{{z}^{{\prime} }}]{a}_{{n}_{z}-{m}_{z}}[0]{u}^{* }(\overrightarrow{q})+{a}_{{n}_{z}}^{* }[0]{a}_{{n}_{z}-{m}_{z}}[-{q}_{{z}^{{\prime} }}]{v}^{* }(\overrightarrow{q}))\\ &&\,\,{F}[{| {{w}}_{{F}}| }^{2}](\overrightarrow{q}+\overrightarrow{m}K)\cdot {F}[{| {{w}}_{{B}}^{\perp }| }^{2}](-({q}_{{x}^{{\prime} }}+{m}_{x}K){\overrightarrow{e}}_{{x}^{{\prime} }}-({q}_{{y}^{{\prime} }}+{m}_{y}K){\overrightarrow{e}}_{{y}^{{\prime} }}),\end{eqnarray} \tag{ 27 }$$

where we use a 3D Wannier wavefunction for the fermions ${{w}}_{{F}}({x}^{{\prime} },{y}^{{\prime} },{z}^{{\prime} })$ but a 2D Wannier wavefunction for the bosons along the two tight axes, ${{{w}}_{{B}}}^{\perp }({x}^{{\prime} },{y}^{{\prime} })$.

Furthermore, to reduce computation time and study long chains, we calibrate on the length-2 chain calculation the acceptable truncation to Wannier functions in Fourier space, that reduces the reciprocal space summation ${\sum }_{{n}_{z}}{\sum }_{\overrightarrow{m}}$ while not affecting the rates within a 1% tolerance.

Each rate ${{\rm{\Gamma }}}_{i\to f}$ requires to numerically sample the spherical coordinates (θ, ϕ) of the emitted excitation wavevectors. We design an approximately homogeneous, isotropic sampling around any point on the sphere as follows.

• (i)  
  The polar angle θ is varied by a fixed step : ${\theta }_{n}=\tfrac{\pi }{{N}_{\theta }}(n+\tfrac{1}{2})$, with n ∈ [0, N_θ − 1].
• (ii)  
  All perimeters on the unit sphere, defined by a polar angle θ_n, are sampled with same angular spacing. This means that the azimuthal angle ϕ sampling depends on θ_n. We furthermore sample each perimeter from a different, random starting point. This sums up as, for each θ_n, an azimuthal sampling ${\phi }_{k}=\mathrm{random}[0,2\pi ]+\tfrac{2\pi }{{N}_{\phi }({\theta }_{n})}k$, with ${N}_{\phi }({\theta }_{n})=\mathrm{Integer}({N}_{\phi }^{0}\sin ({\theta }_{n}))$ and k ∈ [0, N_ϕ(θ_n) −1].

We observe that approximately isotropic sampling around any point on the sphere is obtained for N_ϕ⁰ ≃ 2N_θ. Then, we find that the cooling rates have converged to a $6\times {10}^{-2}$ stability using N_θ = 50. Figures 4–6 are constructed with N_θ varying from 50 to 200, figures 7–10 with N_θ = 50.