Density Shift of Optical Lattice Clock via Multi-Bands Sampling Exact Diagonalization Method

Time, as one of the fundamental quantities, plays an important role in physics. The precise measurement of time is key to the frontier researches including measurement of fundamental constants Blatt et al. (2008), detecting gravitational waves Kolkowitz et al. (2016), testing general relativity Takamoto et al. (2020); Bothwell et al. (2022); Zheng et al. (2022) and searching for the dark matter Derevianko and Pospelov (2014); Arvanitaki et al. (2015). As one of the most accurate time-frequency measurement devices, the optical lattice clock (OLC) platform has made great progress due to the experimental efforts in manipulating alkaline-earth or alkaline-earth-like atoms. Recently, OLC of Ye’s group in JILA has reached an ultra-stable level with a fractional frequency measurement uncertainty $7.6\times 10^{-21}$, so that the gravitational redshift at millimeter-scale can be resolved Bothwell et al. (2022).

To suppress the quantum projection noise, a large number of atoms need to be prepared in the ground state at each optical lattice site T. Akatsuka and Katori (2011). However, the collisions between atoms can degrade the precision of clocks Ludlow et al. (2008), so the corresponding quantitative calculation becomes vital for quantum metrology. On the other hand, understanding the collisional behavior is also important for both quantum simulations Gorshkov et al. (2010); Foss-Feig et al. (2010a, b); Cazalilla et al. (2009); Zhang et al. (2014); Chen and Wu and quantum information Gorshkov et al. (2009); Daley et al. (2008); Hayes et al. (2007) based on OLC. In the early work Blatt et al. (2009), only the s-wave scattering was considered due to the inhomogeneous excitation. Then, people found that the p-wave scattering in fact dominates the density shift Rey et al. (2014), and proposed the inter-atomic collisional model which describes the deep OLC under low longitudinal (z) temperature (around $1\mu K$). In the calculation, only the lowest longitudinal band is taken into account because the longitudinal energy gap is much larger than the temperature effect. Meanwhile, under the collective approximation, the frequency shift is found to be linear with the atom density, and such prediction is checked by the experiments of Ye’s group Martin et al. (2013).

However, sufficient cooling can not always be guaranteed. Considering the rapid progress of the transportableKoller et al. (2017) and space OLCOriglia et al. (2018) which might be difficult to reach the low temperature compared with the lab, it is important to extend the collisional model to multi-bands case, so that the density shift can be estimated with high precision. In this paper, we first derive the multi-bands collisional model of OLC. Then, we developed an advanced algorithm which combines Monte Carlo sampling and exact diagonalization together to calculate the density shift of OLC. Thus, we name it as Multi-Band Sampling Exact Diagonalization (MBSED). In comparison with the collective approximation, the MBSED method presents higher accuracy, especially at higher temperatures. Therefore, our research can provide guidance and a benchmark for OLC working at relatively higher temperatures.

The paper is organized as follows. In Sec. II, we discuss the effective spin model which further takes into account the longitudinal multi-band effect of quantum many-body systems in OLCs. In Sec. III, the algorithm of multi-band sampling exact diagonalization (MBSED) is discussed. In Sec. IV, we compare the MBSED with the collective approximation for both Ramsey and Rabi spectroscopy. At last, Sec. V includes the conclusions and outlook.

We considered an OLC that uses the ¹S${}_{0}(|g\rangle)\leftrightarrow^{3}$P${}_{0}(|e\rangle)$ clock transition in nuclear spin-polarized ⁸⁷Sr atoms. The atoms are loaded in a deep optical lattice constructed by laser with magic wavelength Ye et al. (2008), so that they feel the same lattice potential with inter-site tunneling strongly suppressed. Then, each site can be taken as a 3D slightly anharmonic oscillator with $N$ atoms staying at the different motional bands. Considering the pairwise interaction between atoms through the s-wave and p-wave channels Gorshkov et al. (2009, 2010); Rey et al. (2009); Yu and Pethick (2010); Hazzard et al. (2011); Bishof et al. (2011), The Hamiltonian can be written as:

where $\hat{\psi}_{\alpha}(\bm{r})$ is a fermionic field operator at position $\bm{r}$ for atoms with mass $m$ in electronic state $\alpha$=$g$ or $e$ while $\hat{\rho}_{\alpha}(\bm{r})=\hat{\psi}_{\alpha}^{\dagger}(\bm{r})\hat{\psi}_{\alpha}(\bm{r})$ is the corresponding density operator, and $a_{eg}^{-}$ is the scattering length describing collisions between two atoms in the anti-symmetric electronic state, $\frac{1}{\sqrt{2}}(|ge\rangle-|eg\rangle)$. The p-wave scattering volumes $b_{gg}^{3}$, $b_{ee}^{3}$, $b_{eg}^{3}$ represent three possible electronic symmetric states ($|gg\rangle$, $|ee\rangle$, and $\frac{1}{\sqrt{2}}(|eg\rangle+|ge\rangle)$) respectively. $\omega_{L}$ is the frequency of probing laser and $\bm{k}$ is its wave vector. $\omega_{0}$ is the atomic transition frequency and $\Omega_{0}$ is the bare Rabi frequency.

The previous works Martin et al. (2013); Rey et al. (2014) focus on the system at the low temperature, so that only the lowest longitudinal band needs to be considered. In contrast, turning to the higher temperature, the influence of higher longitudinal bands can not be ignored. Therefore, different from previous studies Rey et al. (2014), the expansion of the field operators in a non-interacting harmonic eigenbasis is not limited to the lowest longitudinal band :

where $\hat{c}_{\alpha\vec{n}}$ annihilates a fermion in motional mode $\vec{n}=\{n_{x},n_{y},n_{z}\}$ and electronic state $\alpha$. Because all trap frequencies are much greater than the characteristic interaction energy, the motional degrees of freedom are effectively frozen Rey et al. (2014). Considering the energy conservation law and anharmonicity of the trap, the external states of atoms could only keep the same or exchange after the collision. Then the atom distribution in the external states is fixed and obtained by sampling as the Boltzmann distribution, so that the effective model can be greatly simplified Rey et al. (2014). After mapping into a spin-1/2 model with $|g\rangle\rightarrow|\downarrow\rangle$ and $|e\rangle\rightarrow|\uparrow\rangle$, the many-body interacting Hamiltonian can be written as:

where $\hat{S}^{\gamma}_{i}$ ($\gamma=x,y,z$) denotes spin operator for $i$th external state, $\delta=\omega_{L}-\omega_{0}$ is the laser detuning from atomic resonance, and $\Omega_{i}$ is the Rabi frequency. The coefficients of the effective spin interactions are explicitly expressed as

where $G^{S}_{i,j}$, $G^{P}_{i,j}$ describes the interaction strength arising from s-wave and p-wave respectively. The derivation of the model in detail can be found in Appendix. A. Although the Hamiltonian $\hat{H}_{s}$ only keeps the spin degree of freedom, the total Hilbert space is still extremely large. One solution is the collective approximation Martin et al. (2013), which takes the average of all the coefficients to replace mode dependent ones. In this approximation, the total spin is confined in the Dicke space with total spin $S=N/2$, see Appendix B for more information about collective approximation. This method is applicable when the fluctuation of all coefficients are small. However, it is no longer valid at higher temperature, because it is necessary to take the multi-bands due to the stronger thermal fluctuations. The strength of the thermal fluctuation can also be reflected by the variation of the Rabi frequency $\Delta\Omega/\overline{\Omega}$. Because the external states of the atoms are sampled according to the Boltzmann distribution, the mean value Rabi frequency can be directly calculated as $\overline{\Omega}=\sum_{i}^{N}\Omega_{i}/N$ and so is the standard deviation $\Delta\Omega$. As demonstrated in Fig.1, $\Delta\Omega/\overline{\Omega}$ rapidly increases while the longitudinal temperature becomes higher, and it is larger than $\approx 0.3$ at $T_{r}\gtrsim 3\mu$K. To include the influence of the higher band at higher temperature, we develop the MBSED method.

The algorithm of the MBSED method is pretty straightforward, and the flow chart diagram is shown in Fig.2. The main processes can be concisely described as follows: (1) numbers of samples are obtained according to the Boltzmann distribution; (2) because the number of particles in each site is not too much, so that time evolution of the quantum state can be calculated via the exact diagonalization. (3) the probability of the excited state can be estimated by taking the mean value of all the samples.

We first use the Monte Carlo method to generate samples that follow the Boltzmann distribution. The corresponding partition function can be represented as:

in which $\omega_{z}$ ($\omega_{r}$) is the longitudinal (transverse) trap frequency, $T_{z}$ ($T_{r}$) is the longitudinal (transverse) temperature, and $n_{r}=n_{x}+n_{y}$ labels the transverse motional state with degeneracy $n_{r}+1$. The number of eigenstates $N_{z}$ ($N_{r}$) in longitudinal (transverse) direction can be estimated with $U_{r}/\hbar\omega_{z}$(($\omega_{r}$) where $U_{r}$ is the depth of lattice potential energy. Typically, $N_{z}$ is around $\approx 5$ and $N_{r}$ is around $\approx 1000$. However, the motional index $(n_{z},n_{r})$ of each atom should be guaranteed that the total energy is less than $U_{r}$. And then, the motional states of all the atoms are sampled according to the distribution Eq.5. Notice that, because of the degeneracy in the transverse direction, the index $(n_{x},n_{y})$ should be randomly selected after fixing $(n_{z},n_{r})$. Meanwhile, the samples with atoms located in the same motional state $(n_{x},n_{y},n_{z})$ have to be kicked out due to Pauli exclusion principle. Such processes are extremely rare at high temperature, and less than one per thousand at 1$\mu$K, so its influence can be omitted.

✓After the motional index $\vec{n}$ is sampled, the corresponding coefficients of Hamiltonian (Eq.3) on the basis of harmonic eigenstates can be directly calculated (see Appendix A for details). Then, all the eigenstates $|\psi\rangle_{m}$ with eigenenergy $E_{m}$ can be obtained by using the exact diagonalization. Assuming the system is prepared at the initial state $|\psi\rangle_{0}$ (i.e. the ground state $|\psi\rangle_{G}=|\downarrow\downarrow...\downarrow\rangle$), the time-dependent wave function turns out to be

However, the full diagonalization requires the atom number $N$ should be $\lesssim 20$, so it is better to find a way to truncate the Hilbert space. The initial state $|\psi\rangle_{G}$ stays at the subspace with total spin $S=\frac{N}{2}$, but evolves into other subspaces due to the interaction and inhomogeneous excitations. The good thing is that the variation of $\langle S\rangle$ is very slow, so that we can project the Hamiltonian to the subspaces with total spin equal to $\{\frac{N}{2},\frac{N}{2}-1,...,\frac{N}{2}-m\}$ when the evolution time is not so long. The selection of the truncation number $m$ depends on parameters of OLC such as temperatures and misalignment angle.

The physical observables can be numerically calculated by taking the mean value of all the samples, such as the probability of the excited state. If the deviation of the mean value is converged to a certain precision, we output the final result, otherwise, the sampling will continue.

The MBSED method is very useful to numerically simulate the spectroscopy with high precision in an OLC beyond the collective regime. To demonstrate its efficiency, the densities of both Ramsey spectroscopy and Rabi spectroscopy are estimated. First, the values of scattering length $b_{ee}$ and $b_{eg}$ should be fitted out by using the MBSED method. Based on the experimental data from Ye’s groupMartin et al. (2013), the fitting scattering lengths are $b_{ee}=150.19a_{B}$ and $b_{eg}=192.34a_{B}$ ($a_{B}$ is the Bohr radius). More details related to the fitting process are shown in Appendix C. Notice that, the scatting lengths we obtained are strongly different from Ref. Martin et al. (2013) (different sign), and we think it may result from the different derivation of the effective model. Therefore, we include the derivation in detail in the supplementary material sup , although it has no effect on our conclusion.

In this paper, we extended the effective spin Hamiltonian which describes the atomic collision of OLC from single band to the multi-band case, so that the OLC system at higher temperatures can be described. To deal with the extended Hamiltonian, we developed a numerical algorithm-MBSED, which combines the Monte-Carlo method to sample the distribution of atoms and exact diagonalization to simulate the time evolving of the OLC system. After implementing the MBSED method on both Ramsey and Rabi spectroscopies, we found the collective approximation method is only valid for the system at low temperatures. For the Ramsey spectroscopy, we found that the linear relation between density shift and atom number still holds, and the special excitation fraction where density shift equals zero decreases with increasing temperature. For the Rabi case, we found that the relation between density shift and atom number is no longer linear. Meanwhile, there is no special excitation fraction where the density shift equals zero, although the density shift becomes smaller when the excitation fraction is higher. Thus we suggest using the highest excitation of the Rabi spectrum for precision measurement.

The quantum many-spin model has been extensively studied in the field of condensed matter physics. In order to solve such quantum many-body problems, numerous numerical algorithms are developed, including exact diagonalization, quantum Monte Carlo, matrix product states, and so on. However, in the field of OLC, there is no trial that borrows these algorithms into solving the quantum many-body Hamiltonian. Thus, our work provides a bridge between OLC and quantum many-body computation. Then, a lot of questions could be asked in future work. One direction is checking how average entanglement entropy grows with the atom number to decide whether the model could deal with the matrix product state or not. Another direction would be further considered in-elastic collision between atoms which has been ignored in this paper.

We wish to thank A. M. Rey for the discussions. This work is supported by the National Science Foundation of China under Grants No. 12274045 and China Postdoctoral Science Foundation Funded Project No. 2020M673118. X.-F. Z. acknowledges funding from the National Science Foundation of China under Grants No. 12274046, No. 11874094 and No.12147102, Fundamental Research Funds for the Central Universities Grant No. 2021CDJZYJH-003. T. Wang acknowledges funding supported by the Program of State Key Laboratory of Quantum Optics and Quantum Optics Devices(No:KF202211).