Rydberg Atoms/Quantum phase transitions

The Rydberg many-body Hamiltonian

We are now interested in the situation where we have a large number of atoms that can be either in the ground state $|g\rangle$ or in a single Rydberg state $|r\rangle$. In contrast to our previous scenarios, we will be interested in the thermodynamic limit, i.e., both the number of atoms and the volume of the system going to infinity. Therefore, the system will be much larger than the blockade radius, and our previous techniques to solve the problem no longer apply.

We will base our analysis on several assumptions:

1. We will not consider any movement of the atoms, we will treat them as being stationary. This requires to laser cool the atoms to very low temperatures (at most a few $\mu {\text{K}}$), this regime is called a frozen Rydberg gas.
2. We will neglect spontaneous emission events, as we assume that the dynamics is much faster than the decay rate of the Rydberg state or the effective decay rate of any intermediate state.
3. We assume that the interaction is described by a van der Waals interaction on all scales. This might be questionable on very short scales, but for such short distances, the interaction strength will be so large that the precise behavior is essentially irrelevant.
4. To simplify the analysis, we assume the atoms to be loaded into an optical lattice, with one atom per lattice site.

Under these assumptions, we can write the many-body Hamiltonian as

$H=-\Delta \sum \limits _{i}|r\rangle \langle r|_{i}+{\frac {\Omega }{2}}\sum \limits _{i}\left(|g\rangle \langle r|_{i}+{\text{H.c.}}\right)-\sum \limits _{i<j}{\frac {C_{6}}{|\mathbf {r} _{i}-\mathbf {r} _{j}|^{6}}}|r\rangle \langle r|_{i}|r\rangle \langle r|_{j}.$

As the many-body system is composed of a large number of two level systems, it is often convenient to rewrite the Hamiltonian in terms of a spin 1/2 description using Pauli matrices. In particular, we have

${\begin{aligned}|r\rangle \langle r|_{i}&={\frac {1+\sigma _{z}^{(i)}}{2}}\\\left(|g\rangle \langle r|_{i}+{\text{H.c.}}\right)&=\sigma _{x}^{(i)},\end{aligned}}$

which allows us to write the Hamiltonian as

$H=-{\frac {\Delta }{2}}\sum \limits _{i}\sigma _{z}^{(i)}+{\frac {\Omega }{2}}\sum \limits _{i}\sigma _{x}^{(i)}-\sum \limits _{i<j}{\frac {C_{6}}{|\mathbf {r} _{i}-\mathbf {r} _{j}|^{6}}}{\frac {1}{4}}\left[1+\sigma _{z}^{(i)}\right]\left[1+\sigma _{z}^{(j)}\right],$

where we have dropped an irrelevant constant energy shift. We can absorb the terms in the interaction containing only a single $\sigma _{z}$ into the definition of the detuning, obtaining

$H=-\sum \limits _{i}{\bigg (}{\frac {\Delta }{2}}+\underbrace {\sum \limits _{j<i}{\frac {C_{6}}{2|\mathbf {r} _{i}-\mathbf {r} _{j}|^{6}}}} _{\Lambda }{\bigg )}\sigma _{z}^{(i)}+{\frac {\Omega }{2}}\sum \limits _{i}\sigma _{x}^{(i)}-\sum \limits _{i<j}{\frac {C_{6}}{4|\mathbf {r} _{i}-\mathbf {r} _{j}|^{6}}}\sigma _{z}^{(i)}\sigma _{z}^{(j)}.$

For primitive lattices, the lattice sum $\Lambda$ is identical for all atoms. In the case of a one-dimensional lattice, we can express its value using the nearest-neighbor interaction energy $V$, finding

$\Lambda _{1D}={\frac {V}{2}}\sum \limits _{j}{\frac {1}{j^{6}}}={\frac {V}{2}}\zeta (6)={\frac {V}{2}}{\frac {\pi ^{6}}{945}},$

where we have used the Riemann $\zeta$ function. In higher dimensions, the value of $\Lambda$ depends on the type of the lattice. For a two-dimensional square lattice, we can express the result as as

$\Lambda _{\mathbf {Z} ^{2}}=2V\beta (3)\zeta (3),$

where $\beta (x)$ is the Dirichlet $\beta$ function .

In the following, we will first focus on the ground state properties of the Rydberg many-body Hamiltonian. This focus will allow us to classify the behavior of the system in terms of various thermodynamic quantum phases. It is possible to prepare the ground state of the Rydberg many-body Hamiltonian for arbitrary values of $\Delta$, $\Omega$ and $V$ by performing an adiabatic preparation. For this, we first start with $\Delta$ being very large, such that the ground state of the many-body Hamiltonian corresponds to all atoms being polarized into their electronic ground state $|g\rangle$. Then, the laser parameters $\Delta$ and $\Omega$ are varied slowly until they reach their desired values, such that the system adiabatically follows the ground state of the many-body Hamiltonian.

The strategy to theoretically study the Rydberg many-body problem depends on whether nearest-neighbor interactions are important. If the blockade radius is smaller than the next-nearest neighbor distance, then the fast decay of the van der Waals interaction allows us to neglect the terms involving next-nearest neighbor and beyond. In the presence of a detuning from the Rydberg state, we can modify the condition for the Rydberg blockade as

${\frac {|C_{6}|}{r_{b}^{6}}}-2\Delta =\Omega .$

If the blockade radius is small enough, then we can write the Hamiltonian as

$H=g\sum \limits _{i}\sigma _{x}^{(i)}+h\sum \limits _{i}\sigma _{z}^{(j)}+V\sum \limits _{\langle ij\rangle }\sigma _{z}^{(i)}\sigma _{z}^{(j)},$

where we have introduced the quantities $g=\Lambda -\Delta /2$ and $g=\Omega /2$. This Hamiltonian describes an Ising model in a transverse field $g$ and a longitudinal field $h$. Especially the purely transverse case with $h=0$ has been studied in great detail in the past.

The transverse Ising model

Let us consider the two possible limits of the transverse Ising model. For $g\to \infty$, we can ignore the interaction, and the problem reduces to a single spin in a magnetic field. Consequently, the ground state is given by

$|\phi _{g}\rangle =|\leftarrow \rangle ^{N}=\prod \limits _{i=1}^{N}{\frac {1}{\sqrt {2}}}(|\uparrow \rangle _{i}-|\downarrow \rangle _{i}).$

In the limit of zero external field, $g\to 0$, we just have to minimize the interaction energy. As the interaction is antiferromagnetic ($V>0$, we have to arrange the spins in an alternating pattern, assuming that the lattice does not lead to geometric frustration. To make the analysis easier, we map the antiferromagnetic Ising model onto the ferromagnetic model by applying a unitary transformation that flips every second spin. After this unitary transformation, the Hamiltonian reads

$H=g\sum \limits _{i}\sigma _{x}^{(i)}-V\sum \limits _{\langle ij\rangle }\sigma _{z}^{(i)}\sigma _{z}^{(j)}.$

In the limit of vanishing $g$, we have two distinct fully polarized ground states, $|\phi _{g}^{(1)}\rangle =|\downarrow \downarrow \downarrow \cdots \rangle$ and $|\phi _{g}^{(2)}\rangle =|\uparrow \uparrow \uparrow \cdots \rangle$. Note that the Hamiltonian does not distinguish between up and down spins (a so-called $\mathbb {Z} _{2}$ symmetry), but the two possible ground states do. This symmetry breaking is a manifestation of the system being in two different quantum phases for $g\gg V$ (ferromagnet) and $g\ll V$ (paramagnet).

To understand more about this quantum phase transition, we have to resort to approximative solutions to the problem, as an exact diagonalization of the Hamiltonian is intractable for anything but a few spins. A simple but powerful technique is mean field theory. For this, we assume that the $\sigma _{z}$ operators occurring in the interaction term are well described by their expectation values and a small fluctuation $\varepsilon$, according to

$\sigma _{z}^{(i)}=\langle \sigma _{z}^{(i)}\rangle +\varepsilon _{i}.$

Inserting this expression into the interaction term, we find

$\sigma _{z}^{(i)}\sigma _{z}^{(j)}=\left[\langle \sigma _{z}^{(i)}\rangle +\varepsilon _{i}\right]\left[\langle \sigma _{z}^{(j)}\rangle +\varepsilon _{j}\right]=\langle \sigma _{z}^{(i)}\rangle \langle \sigma _{z}^{(j)}\rangle +\langle \sigma _{z}^{(j)}\rangle \varepsilon _{i}+\langle \sigma _{z}^{(i)}\rangle \varepsilon _{j}+\varepsilon _{i}\varepsilon _{j}=\langle \sigma _{z}^{(i)}\rangle \langle \sigma _{z}^{(j)}\rangle +\langle \sigma _{z}^{(j)}\rangle \varepsilon _{i}+\langle \sigma _{z}^{(i)}\rangle \varepsilon _{j}+\varepsilon _{i}\varepsilon _{j}+\langle \sigma _{z}^{(i)}\rangle \langle \sigma _{z}^{(j)}\rangle -\langle \sigma _{z}^{(i)}\rangle \langle \sigma _{z}^{(j)}\rangle =\langle \sigma _{z}^{(i)}\rangle \sigma _{z}^{(j)}+\langle \sigma _{z}^{(j)}\rangle \sigma _{z}^{(i)}-\langle \sigma _{z}^{(i)}\rangle \langle \sigma _{z}^{(j)}\rangle +\varepsilon _{i}\varepsilon _{j}.$

The mean field approximation now consists of neglecting the term containing the product of two $\varepsilon$ operators, which describe the quadratic fluctuation around the mean field. The Hamiltonian then decomposes into a sum of single site Hamiltonians $H=\sum _{i}H_{i}-zN\langle \sigma _{z}\rangle ^{2}$ given by

$H_{i}=-z\langle \sigma _{z}\rangle V\sigma _{z}+g\sigma _{x},$

where we have introduced the coordination number $z$ and we have assumed that the magnetization corresponding to the expectation value $\langle \sigma _{z}\rangle$ is identical for all sites. The ground state energy for $H_{i}$ can then be calculated as

$E_{0}=-{\sqrt {g^{2}+(z\langle \sigma _{z}\rangle V)^{2}}}.$

For our mean-field approach to be valid, we also have to demand that the solution is self-consistent, i.e., that the magnetization $\langle \sigma _{z}\rangle$ assumed for all the other sites is identical to the value for site $i$. For the ground state, we find $\langle \sigma _{z}\rangle$ to be given by

$\langle \sigma _{z}\rangle ={\frac {z\langle \sigma _{z}\rangle V}{E_{0}}}=-{\frac {z\langle \sigma _{z}\rangle V}{\sqrt {g^{2}+(z\langle \sigma _{z}\rangle V)^{2}}}}.$

Solving for $\langle \sigma _{z}\rangle$ yields three solutions,

$\langle \sigma _{z}\rangle =0$

$\langle \sigma _{z}\rangle =\pm {\sqrt {1-\left({\frac {g}{zV}}\right)^{2}}}.$

The first solution describes the paramagnetic phase and is the only solution for $g/zV\equiv {\tilde {g}}>1$. For ${\tilde {g}}<1$, the ground state energy is always lower for the other two solutions describing the ferromagnetic phase. As the ground state energy does not depend on the sign of $\langle \sigma _{z}\rangle$, both solutions are degenerate, as expected. The critical point of the quantum phase transition is therefore at ${\tilde {g}}=1$. We can also look into the behavior close to the phase transition. In particular, we find that the magnetization behaves as

$\langle \sigma _{z}\rangle \sim (1-{\tilde {g}})^{1/2}+O([1-{\tilde {g}}]^{3/2}).$

From this expression, we can identify the critical exponent $\beta =1/2$. However, as we have obtained this result within mean field theory, this value is only correct when the coordination number $z$ is sufficiently large. For a $d$-dimensional square lattice we have $z=2d$, and for $d\geq 3$, the mean field exponents are correct. In contrast, in the one-dimensional case, where the transverse Ising model can be solved exactly , we have $\beta =1/8$, while in 2D, Monte-Carlo simulations show a critical exponent of $\beta =0.32644$ . These difference underline the importance of the fluctuations around the mean field in low-dimensional systems.