Efficient variational approach to dynamics of a spatially extended bosonic Kondo model

TL;DR

This paper introduces an efficient variational method for analyzing the dynamics of a quantum spin impurity coupled to a bosonic bath, enabling analytical treatment of complex out-of-equilibrium phenomena in impurity models.

Contribution

The authors develop a variational approach using bosonic Gaussian states and symmetry considerations to simplify and analytically describe the dynamics of a bosonic Kondo impurity model.

Findings

Analytical equations for variational time evolution derived

Method applicable to a wide class of quantum impurity problems

Comparison with Monte Carlo shows new features from orbital motion

Abstract

We develop an efficient variational approach to studying dynamics of a localized quantum spin coupled to a bath of mobile spinful bosons. We use parity symmetry to decouple the impurity spin from the environment via a canonical transformation and reduce the problem to a model of the interacting bosonic bath. We describe coherent time evolution of the latter using bosonic Gaussian states as a variational ansatz. We provide full analytical expressions for equations describing variational time evolution that can be applied to study in- and out-of-equilibrium phenomena in a wide class of quantum impurity problems. In the accompanying paper [Y. Ashida {\it et al.}, Phys. Rev. Lett. 123, 183001 (2019)], we present a concrete application of this general formalism to the analysis of the Rydberg Central Spin Model, in which the spin-1/2 Rydberg impurity undergoes spin-changing collisions in a dense cloud of two-component ultracold bosons. To illustrate new features arising from orbital motion of the bath atoms, we compare our results to the Monte Carlo study of the model with spatially localized bosons in the bath, in which random positions of the atoms give rise to random couplings of the standard central spin model.

Equations

$$\hat{H}=\sum _{n\alpha }{ϵ}_{n\alpha }{\hat{b}}_{n\alpha }^{†}{\hat{b}}_{n\alpha }+{\hat{\mathbf{S}}}_e\cdot {\hat{\mathbf{S}}}_{\mathrm{env}}+h_z{\hat{S}}_e^z,$$

$${\hat{S}}_{\mathrm{env}}^a=\frac{1}{2}\sum _{mn\alpha \beta }{g}_{mn}^a{\hat{b}}_{m\alpha }^{†}{\sigma }_{\alpha \beta }^a{\hat{b}}_{n\beta }$$

$$\hat{U}=\frac{1}{\sqrt{2}}(1+i{\hat{\sigma }}_e^y{\hat{\mathrm{P}}}_{\mathrm{env}}).$$

$$[\hat{\tilde{H}},{\hat{\sigma }}_e^x]=0,$$

$$\hat{\tilde{H}}$$

$${\hat{\tilde{H}}}_0$$

$${\hat{\tilde{H}}}_1$$

$$\mathbf{\varphi }$$

$$\Gamma$$

$$\hat{\mathbf{x}}$$

$$\hat{\mathbf{p}}$$

$$|{\Psi }_{\mathrm{GS}}⟩=e^{i\theta }{e}^{\frac{i}{2}{\hat{\mathbf{\psi }}}^{\mathrm{T}}\sigma \mathbf{\varphi }}{e}^{-\frac{i}{4}{\hat{\mathbf{\psi }}}^{\mathrm{T}}\xi \hat{\mathbf{\psi }}}|0⟩\equiv {\hat{U}}_{\mathrm{GS}}|0⟩,$$

$$\Gamma$$

$$\gamma$$

$$i\hslash {\partial }_t|\Psi (t)⟩=\hat{\tilde{H}}|\Psi (t)⟩,$$

$$i\hslash {\partial }_t|{\Psi }_{\mathrm{GS}}(t)⟩={\hat{P}}_{\partial }\hat{\tilde{H}}|{\Psi }_{\mathrm{GS}}(t)⟩,$$

$${\partial }_t|{\Psi }_{\mathrm{GS}}(t)⟩={\hat{U}}_{\mathrm{GS}}$$

$$\hat{\tilde{H}}|{\Psi }_{\mathrm{GS}}(t)⟩={\hat{U}}_{\mathrm{GS}}$$

$${\mathcal{H}}_{\varphi }$$

$${\mathcal{H}}_{\Gamma }$$

$$\hslash d_t\mathbf{\varphi }$$

$$\hslash d_t\Gamma$$

$$⟨{\hat{\tilde{H}}}_0{⟩}_{\mathrm{GS}}$$

$$⟨{\hat{\tilde{H}}}_1{⟩}_{\mathrm{GS}}$$

$${\tilde{h}}_0=\mathrm{diag}({ϵ}_{n\alpha })+({\sigma }_e^x/4){\sigma }^x\otimes g_{mn}^x.$$

$${\Sigma }_g={\sigma }_e^x{\sigma }^z\otimes g_{mn}^z-i{\sigma }^y\otimes g_{mn}^y$$

$$\Omega$$

$$\hat{\mathbf{b}}$$

$$⟨{\hat{\mathrm{P}}}_{\mathrm{env}}{⟩}_{\mathrm{GS}}=\frac{1}{\sqrt{\det ({\Gamma }_B/2)}}{e}^{-\frac{1}{2}{\mathbf{\varphi }}^{\mathrm{T}}{\Gamma }_B^{-1}(1+\Lambda )\mathbf{\varphi }}.$$

$${\mathcal{H}}_{\varphi }$$

$${\mathcal{H}}_{\Gamma }$$

$$\Phi$$

$$\mathcal{G}$$

$$A(\omega )=\mathrm{Re}\left[{\int }_0^{\infty }dt{e}^{i\omega t}S(t)\right],$$

$$S(t)=⟨\Psi (0)|e^{-i\hat{\tilde{H}}t/\hslash }|\Psi (0)⟩.$$

Full text

I Introduction

Understanding dynamics of a quantum system coupled to a reservoir is a ubiquitous problem in modern physics. Models of this type can be used as a starting point for describing many important phenomena in condensed matter physics, atomic and molecular physics, and optics. Furthermore, coherent control and manipulation of quantum systems coupled to many-body reservoirs are the primary prerequisite for a successful implementation of quantum information processing. Current constraints on the applications of quantum technologies arise from practical limitations on coherence of quantum systems due to their interaction with external environments.

In the context of many-body physics, a problem of particular importance is the model of quantum impurity interacting with a bosonic bath. This class of models can be used to describe a large number of physical systems, including electrons coupled to phonons Landau and Pekar (1948); Feynman (1955); Devreese and Peters (1984), spins coupled to a dissipative environment Leggett et al. (1987), the Bose Kondo problem Falco et al. (2004); Florens et al. (2006); Foss-Feig and Rey (2011); Flottat et al. (2015), electrons in semiconductors Furdyna and Kossut (1988); Lo et al. (1994); Vukmirović et al. (2012), and 3He-4He mixtures Bardeen et al. (1967). Renewed interest in these problems comes from recent experiments with ultracold atoms Schirotzek et al. (2009); Fröhlich et al. (2011); Koschorreck et al. (2012); Meinert et al. (2017); Cetina et al. (2015, 2016); Scazza et al. (2017); Hu et al. (2016); Jørgensen et al. (2016); Yan et al. (2019) that allow one to create many-body systems with tunable parameters and to obtain detailed characterizations of the dynamics. Theoretical studies lead to significant advances in understanding of polaron physics Cucchietti and Timmermans (2006); Schmidt et al. (2012); Efimkin and MacDonald (2017); Tempere et al. (2009); Novikov and Ovchinnikov (2010); Casteels et al. (2011a, b); Rath and Schmidt (2013); Vlietinck et al. (2013); Li and Das Sarma (2014); Shchadilova et al. (2016); Levinsen et al. (2015); Ardila and Giorgini (2015, 2016); Christensen et al. (2015); Ashida et al. (2018a); Peña Ardila et al. (2018) beyond the conventional paradigm established in the previous studies of solid-state systems Alexandrov (2007); Devreese and Alexandrov (2009). An important common feature of these systems is that bath modes are delocalized. Typical examples include phonons in crystal lattices or Bose-Einstein condensates (BEC). In the Bose Kondo model Zhu et al. (2004), the bosonic spinful bath consists of paramagnons of the nearby antiferromagnetic phase and is thus delocalized.

Another important class of quantum baths features localized bath modes, such as nuclear spins that interact with the electron spin on a quantum dot Khaetskii et al. (2003); Schliemann et al. (2003); Coish and Loss (2004). A paradigmatic model for describing such systems is the central spin model Gaudin (1976); Prokof’ev and Stamp (2000); Taylor et al. (2003); Dukelsky et al. (2004); Yuzbashyan et al. (2005); Bortz and Stolze (2007); Zhang et al. (2007); Chen et al. (2007); Lee et al. (2008); Bortz et al. (2010); Witzel et al. (2012); Kessler et al. (2012); Faribault and Schuricht (2013); Rowlands and Lamacraft (2018), which couples a central two-level system to localized (i.e., immobile) environmental spins. This model is also commonly used to describe dynamics of spin qubits realized with phosphorus impurities in silicon Witzel et al. (2010); Tyryshkin et al. (2012) and nitrogen-vacancy centers in diamond Hanson et al. (2008). A particularly intriguing feature of the central spin model is its quantum integrability Gaudin (1976), which implies the existence of an extensive number of integrals of motion. Several powerful techniques have been developed to analyze dynamics of this model, including Bethe ansatz Gaudin (1976); Amico et al. (2001); Bortz and Stolze (2007), the Chebyshev-expansion method Dobrovitski and De Raedt (2003), perturbative approaches Coish and Loss (2004); Fischer and Breuer (2007); Ferraro et al. (2008), effective-Hamiltonian methods Yao et al. (2006); Cywinski et al. (2009), cluster correlation expansion Yang and Liu (2008), the nonperturbative master-equation approach Barnes et al. (2012) and a self-consistent Holstein-Primakoff approximation Kessler et al. (2012). Relation between the central spin model and the BCS theory of superconductivity has also been actively explored Barankov et al. (2004); Yuzbashyan et al. (2005).

Until now the above two classes of many-body systems, namely, quantum impurity problems with mobile and localized bath modes have been analyzed separately. In particular, it is not obvious how to generalize the previous theoretical approaches developed in the central spin problem so as to handle the case of delocalized (i.e., mobile) environmental spins. In this paper, we develop a variational approach to study a generalized type of a quantum impurity problem in which environmental bosons are mobile and spinful while interaction between the localized spin and the surrounding mobile spins is spatially extended. A particularly intriguing aspect of this problem is the interplay between spin dynamics and orbital motion of host particles. A flip of the impurity spin dramatically changes the interaction between the impurity and environmental bosons, and vice versa, orbital motion of bosons affects spin-exchange processes because the interaction depends on the positions of atoms. This unique interplay thus compounds the difficulties of including both orbital motion of mobile bosons and nonlocal spin interactions characteristic of the central spin model. To solve such a class of challenging problems, we describe a new variational approach that first decouples the impurity by employing the canonical transformation Ashida et al. (2018b); *YA18B (at the expense of introducing interactions among the bath particles) and then approximates coherent dynamics of the bath with the bosonic Gaussian states. We provide explicit analytical expressions for the nonlinear time-evolution equations of the vector of coherent expectation values and the covariance matrix of individual modes.

The last two decades have witnessed remarkable developments in giant Rydberg impurities in ultracold atoms, offering a new experimental platform for quantum information processing Saffman et al. (2010) and for realizing new types of strongly correlated many-body states Greene et al. (2000); Hamilton et al. (2002); Bendkowsky et al. (2009); Gaj et al. (2014); Schmidt et al. (2016); Schlagmüller et al. (2016); Camargo et al. (2018); *RS18; Kleinbach et al. (2018). In the accompanying paper Ashida et al. (2019), we point out that Rydberg excitations in ultracold atoms naturally realize the new class of the quantum impurity problem, which we term as the Rydberg Central Spin Model (RCSM). We here discuss in detail how our general theoretical approach can be applied to study nonequilibrium properties of the RCSM. While polaronic effects in Rydberg spectroscopy have recently attracted considerable attention, the focus of previous studies has been either on settings in which only a triplet scattering channel needs to be included and thus spin dynamics is frozen Greene et al. (2000); Hamilton et al. (2002); Bendkowsky et al. (2009); Gaj et al. (2014); Schmidt et al. (2016); Schlagmüller et al. (2016); Camargo et al. (2018); *RS18; Kleinbach et al. (2018), or on low-density regimes, where only few-body molecular physics is relevant Anderson et al. (2014); Saßmannshausen et al. (2015); Böttcher et al. (2016); Deiß et al. (2019); Engel et al. (2019).

This paper is organized as follows. In Sec. II, we present details of the variational approach to studying dynamics of a quantum spin coupled to a generic bosonic environment with the only requirement of having parity symmetry (e.g., symmetry of the rotation by $\pi$ around a certain axis). We introduce a canonical transformation that completely decouples the impurity from the environment. We present the formalism to approximate the time evolution of the environmental wavefunction in the transformed frame using a bosonic Gaussian state as a variational ansatz. We derive analytical expressions of the equations of motion for the variational parameters. In Sec. III, we apply our general approach to analyze nonequilibrium dynamics of the RCSM. We compute absorption spectrum that can be experimentally probed with Rydberg spectroscopy and obtain real-time dynamics of the Rydberg-electron spin. To elucidate the crucial role of the orbital motion of bath atoms, we also analyze the central spin model with random couplings by employing the Monte Carlo sampling. Finally, we summarize the results and present an outlook in Sec. IV.

II General formulation

II.1 Impurity decoupling

In this section, we formulate the variational approach in a general way so that it can be applied to a wide class of quantum impurity problems. Specifically, we consider the following class of Hamiltonians

[TABLE]

where $\hat{b}_{n\alpha}$ ($\hat{b}^{\dagger}_{n\alpha}$) are bosonic annhiliation (creation) operators of environmental modes $n=1,2,\ldots,N_{b}$ of energy $\epsilon_{n\alpha}$ with two internal degrees of freedom $\alpha=\Uparrow,\Downarrow$. We note that environmental modes are not specified here; they can be, for example, momentum eigenstates or single-particle bound states. The localized spin-1/2 operator is denoted by $\hat{{\bf S}}_{e}=\hat{\bm{\sigma}}_{e}/2$. The second term represents a Kondo-type interaction between the localized spin and the environmental spins, in which we introduce the spin-density operator $\hat{\bf S}_{\rm env}$ including couplings as

[TABLE]

with $a=x,y,z$. Here, the $g_{mn}^{a}$ are arbitrary $N_{b}\times N_{b}$ Hermitian matrices labeled by $a$ and represent a generic form of an impurity-environment coupling that can in general be long-ranged and anisotropic. We note that this class of interaction reduces to the standard Kondo coupling Kondo (1964); Hewson (1997); Kanász-Nagy et al. (2018) when it is local in space. The last term in Eq. (1) is a magnetic-field term acting on the localized spin only. This term should be understood as a difference of Zeeman energies between the impurity and the bath spins. We note that the Hamiltonian commutes with the total spin $\hat{S}_{\rm tot}^{z}=\hat{S}_{e}^{z}+\hat{S}_{\rm env}^{z}$, and thus the term proportional to it merely shifts the spectrum by a global constant. In this paper, we set $h_{z}=0$ (except the results plotted in Fig. 4). The magnetic term can be neglected if either the magnetic field is switched off or one uses atoms having the same Zeeman energies.

The first key step in our variational approach is to employ the canonical transformation Ashida et al. (2018b); *YA18B to decouple the localized spin operator. Since the Hamiltonian (1) satisfies the parity symmetry with respect to $\pi$ rotation of the entire system around $z$ axis, it commutes with the parity operator $\hat{{\rm P}}=\hat{\sigma}_{e}^{z}\hat{{\rm P}}_{{\rm env}}$ with $\hat{{\rm P}}_{{\rm env}}=e^{i\pi\hat{N}_{\Uparrow}}$ and $\hat{N}_{\Uparrow}=\sum_{n}\hat{b}_{n\Uparrow}^{\dagger}\hat{b}_{n\Uparrow}$. This conserved parity operator can be mapped onto the localized spin-1/2 operator via $\hat{U}^{\dagger}\hat{{\rm P}}\hat{U}=\hat{\sigma}_{e}^{x}$ with the unitary transformation Ashida et al. (2018b); *YA18B:

[TABLE]

Thus, transforming to the ‘corotating’ frame of the impurity via $\hat{U}$, the localized spin operator turns out to be a conserved quantity and its dynamics freezes, i.e., the system satisfies

[TABLE]

where $\hat{\tilde{H}}\equiv\hat{U}^{\dagger}\hat{H}\hat{U}$ is the transformed Hamiltonian. Specifying the parity-symmetry sector ${\rm P}=\pm 1$, we can treat the spin operator as a classical number $\sigma_{e}^{x}=\pm 1$ in the transformed frame. The resulting Hamiltonian thus only contains the environmental degrees of freedom:

[TABLE]

Here, $\hat{\tilde{H}}_{0}$ represents the quadratic part of the transformed Hamiltonian while $\hat{\tilde{H}}_{1}$ represents interactions among environmental bosons that are in general multibody due to the nonlocality of the operator $\hat{\rm P}_{\rm env}$. The latter appears at the cost of the elimination of the impurity degree of freedom and can be interpreted as effective boson-boson interactions mediated via the impurity spin.

II.2 Variational principle

The next step is to approximate the time-evolved bath state governed by the transformed Hamiltonian $\hat{\tilde{H}}$ by a tractable many-body wavefunction. To this end, we choose a bosonic Gaussian state $|\Psi_{\rm GS}\rangle$ Weedbrook et al. (2012); Shi et al. (2018) as an efficient variational state to describe the bath wavefunction in the transformed frame. The Gaussian states naturally include, as certain subclasses, variational states appropriate for the standard mean-field theory based on, e.g., the Gross-Pitaevskii equation and the Bogoliubov-de Gennes equations.

We here consider a generic Gaussian state, which is parameterized by a $4N_{b}$-dimensional real vector ${\bm{\phi}}$ and a $4N_{b}\times 4N_{b}$ real symmetric matrix $\Gamma$:

[TABLE]

where $\langle\cdots\rangle_{\rm GS}$ denotes the expectation value with respect to the Gaussian state (GS) and $\hat{\bm{\psi}}=(\hat{\bm{x}},\hat{\bm{p}})^{{\rm T}}$ denotes a vector of the quadrature operators $\hat{x}_{n\alpha}=\hat{b}_{n\alpha}^{\dagger}+\hat{b}_{n\alpha}$ and $\hat{p}_{n\alpha}=i(\hat{b}_{n\alpha}^{\dagger}-\hat{b}_{n\alpha})$ as

[TABLE]

and $\delta\hat{\bm{\psi}}=\hat{\bm{\psi}}-\bm{\phi}$ is the fluctuation from the mean-field value. Physically, nonvanishing values of the vector $\bm{\phi}$ indicate macroscopic occupations of certain environmental modes as appropriate for a description of BEC. In addition, the covariance matrix $\Gamma$ in the Gaussian state describes the squeezing of environmental bosons, enabling one to take into account quantum depletion arising from Bogoliubov excitations on top of a simple coherent state. Note that the number of parameters grows at most quadratically with the number of environmental modes.

An explicit form of the Gaussian state is given by

[TABLE]

where $\sigma=i\sigma^{y}\otimes{\rm I}_{2N_{b}}$ with ${\rm I}_{d}$ being the $d\times d$ unit matrix ($d$ is a positive integer) and $\xi$ is a $4N_{b}\times 4N_{b}$ real-symmetric matrix. A matrix dimension is $4N_{b}$ ($2N_{b}$) if a matrix acts on spinful (spinless) single-particle modes. These matrices can be related to the covariance matrix via

[TABLE]

Note that in Eq. (12) we explicitly include a phase factor $\theta$, which is necessary to obtain the absorption spectrum of the system as detailed later.

The time-evolution equation can be obtained from the time-dependent variational principle Jackiw and Kerman (1979); Kramer (2008); Shi et al. (2018). Specifically, we project the exact real-time evolution of the environmental state (in the transformed frame),

[TABLE]

on the manifold spanned by the present variational states. This procedure is equivalent to minimizing the deviation $\epsilon=\|(i\hbar\partial_{t}-\hat{\tilde{H}})|\Psi_{{\rm GS}}(t)\rangle\|^{2}$ from the exact evolution at each moment, leading to

[TABLE]

where $\hat{P}_{\partial}$ is a projection operator onto the tangent space of the variational manifold. The time derivative of the Gaussian state can be rewritten as

[TABLE]

with $c=id_{t}\theta+(i/4){\rm Tr}[\gamma^{{\rm T}}\sigma\left(\partial_{t}\gamma\right)\Gamma]$ being a scalar while the state acted by the transformed Hamiltonian can be represented as

[TABLE]

Here, $:\;:$ indicates the normal order of the bosonic operators $\hat{b}_{n\alpha}$ and $\hat{b}^{\dagger}_{n\alpha}$ and we introduce the functional derivatives of the variational energy as

[TABLE]

The operator $\hat{P}_{\partial}$ projects out the higher-order contributions of $\hat{\bm{\psi}}$, which are denoted as $\delta\hat{O}$ in Eq. (18). Using Eqs. (17) and (18) and comparing linear and quadratic terms in $\hat{\bm{\psi}}$ on both sides of Eq. (16), we obtain the set of variational equations:

[TABLE]

We emphasize that these variational equations should in general offer better results than the standard mean-field theories since the present approach is based on a general Gaussian state and thus can autonomously take into account all the possible two-point correlations on top of an arbitrary BEC state in an unbiased manner.

II.3 Equations of motion

To integrate the real-time evolutions (21) and (22) of the time-dependent variational parameters $\bm{\phi}$ and $\Gamma$, we need analytical expressions of the functional derivatives ${\cal H}_{\phi}$ and ${\cal H}_{\Gamma}$. To this end, we first express the expectation values of the quadratic and interaction Hamiltonians in the transformed frame in terms of the variational parameters as follows:

[TABLE]

where we introduce a $4N_{b}\times 4N_{b}$ real-symmetric matrix ${\cal H}_{0}={\cal S}[({\rm I}_{2}-\sigma^{y})\otimes\tilde{h}_{0}]$ with ${\cal S}[A]\equiv(A+A^{{\rm T}})/2$ being the matrix symmetrization and $\tilde{h}_{0}$ being the single-particle Hamiltonian in the transformed frame (cf. Eq. (6)):

[TABLE]

In the interaction energy (24), we introduce a $2N_{b}\times 2N_{b}$ matrix

[TABLE]

and define the $2N_{b}\times 2N_{b}$ matrix $\Omega$ including the environmental parity operator by Shi et al. (2018)

[TABLE]

Here, we represent the vector of the bosonic creation operators of different environmental modes as

[TABLE]

and introduce matrices $\Sigma_{z}=\sigma^{z}\otimes{\rm I}_{N_{b}}$, $\Lambda={\rm I}_{2}\otimes\Sigma_{z}$ and $\Gamma_{B}=({\rm I}_{4N_{b}}+\Lambda)\Gamma+{\rm I}_{4N_{b}}-\Lambda$. The expectation value of the environmental parity operator is given by

[TABLE]

Calculating the derivatives of the variational energy $\langle\hat{\tilde{H}}\rangle_{\rm GS}$ with respect to the parameters $\bm{\phi}$ and $\Gamma$ in Eqs. (19) and (20), we now obtain the analytical expressions of ${\cal H}_{\phi}$ and ${\cal H}_{\Gamma}$ as follows:

[TABLE]

where we introduce the $4N_{b}\times 4N_{b}$ matrices

[TABLE]

We note that the variational equations (21) and (22) together with the analytical formulae (32) and (33) are general and can be readily applied to studying out-of-equilibrium dynamics of a quantum system coupled to various types of environments. While we focus on the real-time evolution in this paper, the ground-state properties can also be analyzed by using Eqs. (32) and (33) to integrate the variational imaginary-time evolution whose explicit form is given in App. A.

II.4 Absorption spectrum

Out-of-equilibrium properties can be studied by analyzing the absorption spectrum Schmidt et al. (2018b)

[TABLE]

where $S(t)$ is the overlap between the initial state and a time-evolved state

[TABLE]

For the sake of simplicity, a contribution from the free time evolution (e.g., the one without Rydberg interactions in the model discussed later) is not included here as it just shifts $A(\omega)$ by a trivial constant. To calculate $S(t)$ in the present variational approach, we have to obtain an equation of motion of the phase factor $\theta$ of the Gaussian environmental state in addition to that of the mean-field vector $\bm{\phi}$ and the covariance matrix $\Gamma$. To this end, the most convenient way is to parameterize the Gaussian state in the basis of $\hat{b}_{n\alpha}$ and $\hat{b}^{\dagger}_{n\alpha}$ operators as follows:

[TABLE]

where $\Xi$, $\Xi^{\prime}$ and $\Xi^{\prime\prime}$ are $2N_{b}\times 2N_{b}$ matrices. The overlap can then be obtained from

[TABLE]

where $\delta{\bm{\beta}}$ is the difference mean-field vector in the basis of $(\hat{\bm{b}},\hat{\bm{b}}^{\dagger})^{\rm T}$ defined as follows:

[TABLE]

with $B$ being the $4N_{b}\times 4N_{b}$ matrix

[TABLE]

The time evolution of the variational parameters $\delta{\bm{\beta}}$ is obtained by integrating the equations of motion (21) to yield ${\bm{\phi}}(t)$ and by applying Eq. (43) to transform it into an expression in the basis of $\hat{b}_{n\alpha}$ and $\hat{b}^{\dagger}_{n\alpha}$ operators. Following the same procedures relying on the time-dependent variational principle as discussed in Sec. II.2, we obtain the set of time-evolution equations for the parameters $\theta$ and $\Xi$,

[TABLE]

Here we denote $\delta\bm{\phi}(t)={\bm{\phi}}(t)-{\bm{\phi}}(0)$ and define the $2N_{b}\times 2N_{b}$ matrices $h_{\Gamma}$ and $h^{\prime}_{\Gamma}$ by

[TABLE]

Integrating Eqs. (47) and (48) together with Eqs. (32) and (33), we calculate the overlap $S(t)$ from Eq. (41) whose Fourier transform provides the absorption spectrum (38).

III Application to the Rydberg Central Spin Model

III.1 Variational results with orbital motion

We now discuss the application of our general variational approach to a new class of the quantum impurity problem, in which bath bosons are mobile and spinful while interaction between the localized spin and the surrounding particles is spatially extended. This class of many-body problems can be naturally realized with Rydberg gases Ashida et al. (2019). Rydberg excitations have been powerful probes for studying new regimes of out-of-equilibrium dynamics of molecular states and quantum many-body systems Greene et al. (2000); Hamilton et al. (2002); Bendkowsky et al. (2009); Gaj et al. (2014); Schmidt et al. (2016); Schlagmüller et al. (2016); Camargo et al. (2018); *RS18; Kleinbach et al. (2018); Anderson et al. (2014); Saßmannshausen et al. (2015); Böttcher et al. (2016); Deiß et al. (2019); Engel et al. (2019). When an atom in a BEC is excited into a Rydberg state with a high principal quantum number $n$, the surrounding environmental bosons are subject to Fermi’s pseudopotential Fermi (1934) created by frequent elastic scattering of the weakly bound Rydberg electron. If the electron-atom scattering length is negative, the interaction between the Rydberg electron and the environmental atoms is attractive and the atoms can become bound to the Rydberg excitations. This can lead to giant molecular states known as Rydberg molecules whose size can be several thousand of the Bohr radius Greene et al. (2000); Bendkowsky et al. (2009).

We consider a situation in which internal dynamics of environmental bosons is restricted to two internal states $\alpha=\Uparrow,\Downarrow$. Thus, the bosons act as mobile pseudospin-1/2 particles and their scattering with the Rydberg electron separates into two channels. We denote these scattering channels as $1$ and $2$, and the corresponding pseudopotentials as $V_{1,2}$. The resulting Hamiltonian is (cf. Ref. Ashida et al. (2019))

[TABLE]

where $\hat{\Psi}_{{\bf r}\alpha}$ ($\hat{\Psi}_{{\bf r}\alpha}^{\dagger}$) are annihilation (creation) operators of bosons at position $\bf r$ with internal state $\alpha=\Uparrow,\Downarrow$, $\hat{\bf S}_{e}=\hat{\bm{\sigma}}_{e}/2$ is the spin operator of the Rydberg electron, $h_{z}$ is a magnetic field acting on the impurity spin, and $\hat{\bf S}_{\bf r}$ is the environmental spin density

[TABLE]

We introduce ${h}_{0}$ as the matrix elements of a single-particle Hamiltonian

[TABLE]

where $m$ denotes the atomic mass. The potential $V_{0}$ and the long-range Kondo coupling $g_{\bf r}^{a}$ are in general given by linear combinations of the two Rydberg potentials $V_{1,2}$, in which the coefficients depend on the internal structure of environmental atoms. As appropriate for the simplest setup of the scattering between two spin-1/2 particles, hereafter we choose $V_{0}=(3V_{1}+V_{2})/4$ and $g_{{\bf r}}^{z}=g_{\bf r}^{\parallel}\equiv V_{1}-V_{2}$ with $V_{1}$ and $V_{2}$ being the triplet- and singlet-pseudopotentials characterizing the interaction between the Rydberg electron and the surrounding atoms. As a representative case, we use the data for the potential profiles of the scattering between the electron of the 87Rb$(87s)$ state and the ground-state 87Rb atoms; all the numerical results presented below are obtained for the triplet and singlet Rydberg potentials of Rb atoms with $n=87$. Yet, in practice we envision using two-electron environmental atoms (e.g., Sr atoms), in which the hyperfine coupling can be neglected as discussed in the accompanying paper Ashida et al. (2019). The potential profiles in this two-electron atomic setup are analogous to those in the Rb setup used here Killian and Bartschat . Thus, essential features of the RCSM discussed below, which are qualitatively insensitive to these specific choices of potentials, will remain in the two-electron atomic setup. Meanwhile, a possible difference in using the two-electron atom is the anisotropy in the spin interaction Ashida et al. (2019). To address this effect, we set $g_{{\bf r}}^{x,y}=g_{\bf r}^{\perp}\equiv\sqrt{2}g_{{\bf r}}^{z}$. We neglect a temperature effect, boson-boson interactions and the ion recoil because the Rydberg potentials typically exceed those energy scales. We also omit contributions from $p$-wave scattering.

Previous works on Rydberg molecule and polaron formation in a high-density regime Greene et al. (2000); Hamilton et al. (2002); Bendkowsky et al. (2009); Gaj et al. (2014); Schmidt et al. (2016); Schlagmüller et al. (2016); Camargo et al. (2018); *RS18; Kleinbach et al. (2018) have exclusively focused on the following polarized initial state

[TABLE]

where only the triplet channel is relevant to the dynamics. Thus, the interaction part of Eq. (III.1) plays no role and the resulting Hamiltonian is quadratic. In this case, the spin dynamics is completely frozen and the orbital motion of environmental atoms solely characterizes the nonequilibrium properties of the system.

In the present work, we are interested in a different initial state with the Rydberg spin and the surrounding bosons being polarized into the opposite directions:

[TABLE]

We note that such an initial condition has already been realized in the experiments at low-density regimes Anderson et al. (2014); Saßmannshausen et al. (2015); Böttcher et al. (2016); Deiß et al. (2019); Engel et al. (2019). In this setting, the central spin-type interaction (cf. the second term on the right-hand side of Eq. (III.1)) now plays an important role and the nonequilibrium properties of the system are characterized by the interplay between the spin and orbital degrees of freedom of the environmental particles.

The main difficulty here is taking into account mesoscopic collective response of the orbital motion as well as the entanglement of the impurity and environmental spins. The main idea of the variational approach presented in this paper is to employ the disentangling transformation (3) to decouple the impurity spin and reduce the problem to that of the bath degrees of freedom with additional interactions. We then approximate the evolution of the environmental wavefunction in the transformed frame by the Gaussian state, which guarantees that all two-particle correlations are taken into account.

To apply the general formulation presented in the previous section, it is useful to introduce single-particle energy eigenstates as the computational basis:

[TABLE]

where $n=1,2,\ldots,N_{b}$ is the nodal quantum number of the single-particle state $\psi_{n{\bf r}}$. While we focus on the situation in which only states with zero angular momentum $l=0$ are relevant, the formulation can be easily generalized to include effects of nonzero angular momentum. Further details about the calculations of the eigenenergies $\epsilon_{n}$ and the wavefunctions $\psi_{n{\bf r}}$ are given in App. B. The central spin couplings in this basis are given by

[TABLE]

with $a=x,y,z$. These identifications lead to the Hamiltonian (1), and Eqs. (5), (6) and (7) after performing the decoupling transformation (3). The initial state (56) corresponds to the following initial values of the variational parameters:

[TABLE]

where $N$ is the total number of environmental bosons, which can be related to the density $\rho$, and $\psi_{{\rm ini},{\bf r}}$ is the initial lowest-energy wavefunction without Rydberg potentials (see App. B for further details). We note that this value represents the number of bosons in the box used to calculate eigenstates rather than the number of bosons in the Rydberg potential. Because the initial state (56) resides in the parity sector $\rm P=+1$, we can treat the Rydberg spin as a classical number $\sigma_{e}^{x}=+1$ in the transformed frame.

To obtain values of the time-dependent variational parameters, we integrate the variational time-evolution equations (21), (22), (47) and (48) with the use of the analytical expressions of the functional derivatives (32) and (33). We then calculate the overlap $S(t)$ from Eq. (41) and obtain its Fourier spectrum $A(\omega)$ from Eq. (38). The magnetization dynamics of the central spin

[TABLE]

can be obtained from Eq. (31). To improve the numerical accuracy, we find it useful to implement a penalty term to ensure the spin conservation (see App. C for details). A typical number $N_{b}$ of single-particle states necessary to converge the variational results is an order of several tens of states; for instance, $N_{b}\sim 80$ is sufficient for the parameters considered here.

Figure 1 shows the results obtained for the overlap $S(t)$ and the corresponding absorption spectrum $A(\omega)$ at different densities $\rho$. As shown in the top panels in Fig. 1, we find a fast decay in $S(t)$ that signals the rapid dephasing dynamics due to the creation of high-energy excitations of the environmental bosons. The fast decay of $S(t)$ in the time domain corresponds to the emergence of a Gaussian profile of the Fourier spectra $A(\omega)$ in the frequency domain in the high-density regime as shown in the bottom panels in Fig. 1. As indicated by the dashed lines in the same panels, the centers of the Gaussian profiles agree with the following mean-field shift which neglects the flip-flop coupling:

[TABLE]

From the initial condition (56), we can set the central spin to be a classical number $S_{e}^{z}=+1/2$. This approximation of neglecting the flip-flop interaction with setting $S_{e}^{z}=+1/2$ can be a useful starting point since the central spin is largely polarized during the dynamics in the case of a dense polarized spin bath Bortz and Stolze (2007). The resulting longitudinal Hamiltonian $\hat{H}_{\parallel}$ is quadratic in terms of the bosonic operators, where free bosons feel the mean-field potential $V_{\rm mean}=(V_{T}+V_{S})/2$. Thus, the emergent Gaussian feature can be understood as a consequence of a large number of atoms stochastically occupying single-particle eigenstates of the mean-field Hamiltonian $\hat{H}_{\parallel}$, which defines one of the key features of the Rydberg polaron Camargo et al. (2018); *RS18; Kleinbach et al. (2018). We note that this feature is determined entirely by the longitudinal couplings $g_{\bf r}^{z}$ and is thus not influenced by the anisotropy in the spin-interaction term. We also remark that the peaks in $A(\omega)$ in the limit of $\rho\to 0$ eventually match with bound-state energies for each of the Rydberg potentials, $V_{T}$ and $V_{S}$ Bendkowsky et al. (2009); Gaj et al. (2014); Schmidt et al. (2016).

The overlap $S(t)$ also accompanies oscillations, which can be best understood from the corresponding Fourier spectrum $A(\omega)$. For a noninteracting setup described by a quadratic Hamiltonian (such as the triplet setting (55) in Refs. Bendkowsky et al. (2009); Gaj et al. (2014); Schmidt et al. (2016); Schlagmüller et al. (2016); Camargo et al. (2018); *RS18; Kleinbach et al. (2018)), the frequency spacing of peaks matches the single-particle energy of the dominant bound state that has the largest overlap with the initial single-particle wavefunction. Correspondingly, the sharp peaks in the spectra indicate the formation of Rydberg molecular states, where one, two or more environmental atoms are bound to and localized in the outermost well of the molecular potential. In contrast, in the present case of an interacting system, the central spin interaction causes the formation of a correlated many-body bound state. While the existence of sharp peaks in $A(\omega)$ still has its root in molecular physics, the simple explanation based on single-particle energies is no longer applicable to understand the underlying rich structures such as positions and spacings of the peaks and their sensitivity to the environmental density Ashida et al. (2019). More specifically, the many-body feature manifests itself as the spin-interaction-induced renormalization of $A(\omega)$ from that of noninteracting, bare molecular states (see Fig. 2). These results are qualitatively insensitive to specific details of the Rydberg potentials and thus should remain in a different choice of a Rydberg state as long as it has no angular momentum. In particular, choosing a small principal quantum number $n$ for a Rydberg state, one can obtain a better resolution of the peaks in $A(\omega)$ whose energy scale can be an order of, e.g., MHz in contrast to a rather small energy scale in Fig. 2.

In the noninteracting triplet setting, the anisotropy does not affect values of single-particle energies of the mean-field Hamiltonian (66) and thus the spectra remain the same as in the isotropic case Ashida et al. (2019). In the present interacting problem, quantitative values of many-body bound-state energies can depend on the anisotropy since the interaction term includes the transverse couplings $g_{\bf r}^{\perp}$, however, the qualitative features remain the same.

Figure 3(a) shows the corresponding spin dynamics $m_{z}(t)$ of the Rydberg electron at different densities $\rho$. Due to a large energy cost to flip the central spin coupled to a polarized environment, only a small fraction of a many-body state with the opposite central spin $|\!\downarrow\rangle_{e}$ can be admixed. The resulting large polarization of the Rydberg spin again indicates that the flip-flop coupling does not play a significant role at the level of the mean-field feature of the absorption spectrum (cf. Eq. (65)).

Remarkably, the spin dynamics is accompanied by fast, long-lasting oscillation that is reminiscent of the nondecaying oscillations found in the central spin problem Bortz and Stolze (2007). To further study this feature, we plot the Fourier spectra $\tilde{m}_{z}(\omega)$ of the spin dynamics in Fig. 3(b). The oscillation frequencies become higher as the density $\rho$ is increased and exhibit a square root scaling $\omega=a\sqrt{\rho}$ as indicated by the dashed lines in Fig. 3(b). We find that the proportionality constant $a$ is larger than the corresponding value in the isotropic case Ashida et al. (2019) by the anisotropic factor $g_{\bf r}^{\perp}/g_{\bf r}^{\parallel}$ (which is chosen to be $\sqrt{2}$ in the present numerical calculation). This fact indicates that the frequency of the central spin oscillation is mainly characterized by the transverse spin-exchange couplings $\hat{H}_{\perp}=(\hat{S}_{e}^{+}\int d{\bf r}g_{\bf r}^{\perp}\hat{S}_{\bf r}^{-}+{\rm H.c.})/2$. The square root scaling $\sqrt{\rho}$ can then be interpreted as a BEC enhancement factor that arises when the bosonic annihilation operator $\hat{b}_{\Downarrow}$ in $\hat{H}_{\perp}$ acts on a macroscopically occupied single-particle state of environmental bosons. Another interesting feature in Fig. 3(a) is that the oscillations exhibit the collapse and revival in accordance with the motional wavepacket dynamics of bath atoms, whose time scale is roughly on the order of the binding energy. We note that a large scattering-length difference is crucial to observe these many-body features as it directly relates to the spin-exchange couplings $g_{\bf r}^{\perp}$.

Atoms coupled to the central spin via strong flip-flop interactions are also subject to the strong orbital potential and thus typically ocuupy the dominant bound state localized in the outermost well of the Rydberg potential. We note that the observed square-root scaling with density is distinguished from the conventional linear scaling found in the standard central spin problem Khaetskii et al. (2003); Schliemann et al. (2003); Coish and Loss (2004), where the environmental spins are immobile and thus a macroscopically occupied state is absent. This weaker dependence of the oscillation frequency on a density is advantageous when one attempts to experimentally realize the predicted oscillation in atom clouds of inhomogeneous density. In practice, the magnetization dynamics might be measured by the Stern-Gerlach-like experiment. For instance, one may send back the Rydberg excitation to, e.g., the 6p state whose optical transition conserves the spin number, and then the transitioned atoms can be separated by using optical dipole traps. Measuring the population of the separated atoms with, for instance, using the ion microscopes, one could obtain information about the magnetization dynamics.

The spin oscillation is strongly influenced by controlling an external magnetic field. Figures 4(a) and (b) show the spectra $\tilde{m}_{z}(\omega)$ and the real-time evolution $m_{z}(t)$ of the dynamics of the Rydberg electron spin at different magnetic fields $h_{z}$, respectively. The oscillation frequencies of the spin oscillation can be controlled by changing $h_{z}$. We find that the shift $\delta\omega$ of the oscillation frequency from the zero-field value scales linearly with magnetic field, $\delta\omega=-h_{z}$, as indicated by the dashed lines in Fig. 4(a). The amplitude of the oscillation is enhanced for a positive field $h_{z}>0$ for which the resonance is approached, while it is suppressed for the opposite sign of the field $h_{z}<0$. These findings are consistent with the previously known magnetic-field dependences in the standard central spin problem Bortz and Stolze (2007). This fact shows that it is possible to manipulate the electron spin of dense Rydberg gases in the similar way as in solid-state qubits Khaetskii et al. (2003); Schliemann et al. (2003); Coish and Loss (2004); Hanson et al. (2008); Witzel et al. (2010); Tyryshkin et al. (2012). Furthermore, our findings also suggest that the Rydberg electron may be used to prepare and manipulate a mesoscopic, delocalized spin environment in a way analogous to localized nuclear-spin environments Koppens et al. (2005).

III.2 Monte Carlo results without orbital motion

To elucidate the importance of taking into account the orbital motion of the environmental particles, we analyze the case of infinite-mass environmental particles ($m\to\infty$). In this limit, the atomic orbital motion is completely frozen, and the position operators of the environmental bosons commute with the Hamiltonian. For each set of positions of environmental atoms $\{{\bf r}_{i}\}$, we can thus reduce the problem to the ordinary central spin problem:

[TABLE]

where $V=\sum_{i=1}^{N}V_{0}({\bf r}_{i})$ and we denote the central spin couplings as $g_{i}^{\parallel}=g_{{\bf r}_{i}}^{z}$ and $g_{i}^{\perp}=g_{{\bf r}_{i}}^{x,y}$. Because of the conservation of the total spin-$z$ components, the time-evolved wavefunction is restricted to the sector of $\hat{\sigma}_{e}^{z}+2\sum_{i}\hat{S}_{i}^{z}=1-N$ and can be parametrized as

[TABLE]

The time-evolution equations of the variables $\xi$ are given by

[TABLE]

where $G=\sum_{i=1}^{N}g_{i}^{\parallel}$. This equation can be solved analytically by performing the Laplace transformation

[TABLE]

Using the initial condition $\xi_{0}(0)=1$ and $\;\xi_{i}(0)=0$ for $i=1,2,\ldots,N$, we obtain the analytical expression for $\xi_{0}$,

[TABLE]

where $\omega_{0}=(V+G/4)/\hbar$ and the contour ${\cal C}$ is chosen so that all the poles in the integral lie to its above; one can choose the contour that extends from $-\infty$ to $\infty$ along the real axis and is closed via a half circle in the upper complex plane. After performing the integration, we find

[TABLE]

Here, $w_{l}=1/[1+\sum_{i=1}^{N}g_{i}^{\perp 2}/(2\omega_{l}+g_{i}^{\parallel})^{2}]$ and $\left\{\omega_{l}\right\}_{l=1}^{N+1}$ are the poles enclosed by the integration contour, which can be obtained from the algebraic equation

[TABLE]

We note that Eq. (75) reproduces the previous result Bortz and Stolze (2007) for obtaining the Bethe roots in the case of isotropic couplings $g_{i}^{\perp}=g_{i}^{\parallel}$. The overlap is given by

[TABLE]

leading to the absorption spectrum

[TABLE]

Similarly, the central spin dynamics is given by

[TABLE]

To calculate the absorption spectrum appropriate for the initial state (56), we randomly generate sets of atomic positions {${\bf r}_{i}$} according to the initial BEC wavefunction $\prod_{i=1}^{N}\psi_{{\rm ini},{\bf r}_{i}}$ and obtain the corresponding $A^{\{{\bf r}_{i}\}}(\omega)$ from Eq. (77). This procedure is repeated to calculate the ensemble average over different atomic configurations to yield

[TABLE]

Similarly, the central spin dynamics is obtained from Eq. (78) and reads

[TABLE]

Figures 5(a) and (b) show the calculated absorption spectrum $A(\omega)$ and the corresponding central spin dynamics $m_{z}(t)$ for the infinite-mass environmental particles at typical parameters. The results are obtained by taking the ensemble average over $10^{5}$ different stochastic realizations of atomic configurations. As shown in Fig. 5(a), the absorption spectrum exhibits the Gaussian feature which is the consequence of the average over various occupations of energy eigenstates of the stochastic central spin Hamiltonian (67) (cf. Eqs. (77) and (79)). We note that, to obtain the spectrum, one still needs to solve the (integrable) many-body problem (67), which should be distinguished from the single-particle calculations performed in, e.g., the triplet setting (55). To illustrate such a quantum aspect of the present problem, we plot typical single realizations of the central spin dynamics with the fixed spin couplings (see the dashed thin curves in Fig. 5(b)). For each single realization, the spin exhibits nondecaying oscillations as previously discussed in the literature on the central spin coupled to polarized spin environments Bortz and Stolze (2007). However, after taking the ensemble average over different atomic configurations (cf. Eq. (80)), these oscillations are averaged out, and the actual spin dynamics quickly relaxes to a stationary finite value as shown by the thick curve in Fig. 5(b).

We note that the key features found in the previous subsection, in which both the spin dynamics and the orbital motion are taken into account, are absent in the present infinite-mass treatment. For instance, the peaks in the absorption spectrum do not appear as the Monte Carlo analysis cannot describe the formation of Rydberg molecular state. Also, the characteristic spin-precession dynamics is found to be absent in the infinite-mass treatment. These facts again highlight the importance of orbital motion of environmental particles to understand the dynamics of the RCSM. We remark that it would be useful to work out in a different single-particle basis such as position eigenstates for doing a direct comparison between the Monte Carlo analysis and the variational method in the large-mass regime.

IV Conclusions and outlook

We developed an efficient variational approach to solving quantum spin systems coupled to bosonic environments and applied this technique to analyze nonequilibrium properties of the Rydberg Central Spin Model (RCSM) proposed in the accompanying paper Ashida et al. (2019). The key element of this approach is the decoupling of the impurity spin using the canonical transformation (3). This transformation maps the conserved parity operator into one of the components of the impurity-spin operator and thus makes this spin component an integral of motion. In the transformed frame, we then only need to consider the wavefunction of the bosons, which we approximate by the bosonic Gaussian state. This wavefunction includes all two particle correlations and requires the number of parameters growing only quadratically with the system size. We utilized time-dependent variational principle to derive analytical expressions for the variational time-evolution equations (see Eqs. (21) and (22) together with (32) and (33)). These equations are general and can be readily applied to analyze the in- and out-of-equilibrium properties of a wide class of systems, in which quantum impurities are coupled to bosonic environments. As concrete examples, we applied our theory to analyze several dynamical aspects of the RCSM. In particular, we predicted sharp peaks in the absorption spectrum corresponding to dressed many-body molecular states, and long-lasting central spin oscillations. To elucidate the crucial role of the orbital motion of atoms in these two features, we also analyzed the situation of immobile (i.e., infinitely heavy) bath atoms located at random positions. In this limit, the problem can be solved exactly for each stochastic realization; we then used the Monte Carlo sampling to take an ensemble average over possible atomic configurations. While the results are obtained by using the Rydberg potential of Rb atoms as an illustration, the qualitative physics found here will remain the same as long as the mean-field potential $V_{\rm mean}$ is attractive and supports a bound state, and also non-$s$-wave contributions do not play significant roles. Moreover, owing to the versatility of the developed variational approach, one can readily obtain quantitative predictions appropriate for concrete setups by using specific experimental parameters.

Our study suggests several promising directions for future research. Firstly, in the context of solid-state systems, it is intriguing to apply the present variational approach to study the Bose Kondo problems Falco et al. (2004); Florens et al. (2006); Foss-Feig and Rey (2011); Flottat et al. (2015) in both in- and out-of-equilibrium regimes. In the field of ultracold atoms, our analysis can be extended to make predictions of site-resolved many-body dynamics as appropriate for experiments using quantum gas microscopy Endres et al. (2011); Cheneau et al. (2012); Schauß et al. (2012); Fukuhara et al. (2015); Kaufman et al. (2016); Miranda et al. (2015); Ashida and Ueda (2015); *YA16OL; *YA17multi; Yamamoto et al. (2016); Hollerith et al. (2018). The covariance matrix obtained in our approach can be expressed in the real-space basis, from which one can readily extract spatiotemporal dynamics of correlation functions. It is intriguing to extend the present approach to analyze dissipative dynamics of quantum systems coupled to generic bosonic baths at finite temperatures. This can be done by generalizing the variational state to Gaussian density matrices. Another interesting direction is to extend our approach to Markovian open quantum systems subject to dissipation Daley (2014); Kessler et al. (2012); Barontini et al. (2013); Ashida et al. (2018d); *ZG18 and measurements Carmichael (1993); Lee and Chan (2014); Labouvie et al. (2016); Patil et al. (2015); Ashida et al. (2017); *YA18fcs by employing the variational principle appropriate for master-equation dynamics Cui et al. (2015); Weimer (2015). Generalizing the decoupling canonical transformation, our approach can also be applied to multiple impurities Sala et al. (2018) coupled to bosonic environments. In these ways, our approach can serve as a powerful theoretical tool to study a multitude of interesting quantum many-body problems.

Secondly, it is intriguing to study the RCSM in a fermionic setup, where the interplay between the conventional central spin problem Gaudin (1976); Prokof’ev and Stamp (2000); Taylor et al. (2003); Dukelsky et al. (2004); Yuzbashyan et al. (2005); Bortz and Stolze (2007); Zhang et al. (2007); Chen et al. (2007); Lee et al. (2008); Bortz et al. (2010); Witzel et al. (2012); Kessler et al. (2012); Faribault and Schuricht (2013); Rowlands and Lamacraft (2018) and fermionic Kondo physics Hewson (1997); Riegger et al. (2018); Ono et al. (2019) should offer another distinct class of many-body problems. Thirdly, our study suggests a promising direction of using Rydberg electrons to prepare and manipulate mesoscopic spin environments. While this technique has been implemented in solid-state devices Koppens et al. (2005) with a localized nuclear-spin environment, it remains an open question to what extent this method can be generalized to setups with delocalized (i.e., mobile) bath particles. Finally, it merits further study to explore the impact of other experimental details arising from, e.g., nonzero angular momentum of the Rydberg electron and bath bosons or the hyperfine couplings. These effects are expected to be important in several types of atoms Anderson et al. (2014); Saßmannshausen et al. (2015); Böttcher et al. (2016); Engel et al. (2019). We hope that our work stimulates further studies in these directions.

Acknowledgements.

We are grateful to Shunsuke Furukawa, Tom Killian, Jesper Levinsen, Meera Parish, Masahito Ueda, and Shuhei Yoshida for fruitful discussions. Y.A. acknowledges support from the Japan Society for the Promotion of Science through Program for Leading Graduate Schools (ALPS) and Grant Nos. JP16J03613 and JP19K23424, and Harvard University for hospitality. T.S. acknowledges the Thousand-Youth-Talent Program of China. R.S. is supported by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany’s Excellence Strategy – EXC-2111 – 390814868. J.I.C. is supported by the ERC QENOCOBA under the EU Horizon 2020 program (grant agreement 742102). E.D. acknowledges support from Harvard-MIT CUA, NSF Grant No. DMR-1308435, AFOSR Quantum Simulation MURI, AFOSR grant number FA9550-16-1-0323, the Humboldt Foundation, and the Max Planck Institute for Quantum Optics.

Appendix A Imaginary-time evolution

We here provide the variational equations for the imaginary-time evolution, which can be useful to analyze ground-state properties. The exact form of the imaginary-time evolution in the transformed frame is given by

[TABLE]

Minimizing the error of the variational evolution from the exact one, we arrive at the following differential equation (cf. discussions above Eq. (16)):

[TABLE]

Using the expressions (17) and (18) and comparing the linear and quadratic terms of $\hat{\bm{\psi}}$ on both sides of Eq. (82), we obtain

[TABLE]

Integrating these variational equations together with the analytical expressions (32) and (33) of ${\cal H}_{\phi}$ and ${\cal H}_{\Gamma}$, one can analyze the ground-state properties in the limit of $\tau\to\infty$.

Appendix B Rydberg wavefunction

We here briefly outline the calculation of single-particle energy eigenstates of atoms subject to a Rydberg potential. We first note that, since the Hamiltonian is rotationally symmetric and the initial single-particle wavefunction $\psi_{{\rm ini},{\bf r}}$ has zero angular momentum, the system resides in the sector with angular momentum $l=0$ during the course of the time evolution. The resulting radial single-particle Schrödinger equation is thus given by

[TABLE]

where $r$ denotes the radial coordinate, $V(r)$ is the Rydberg potential, and we introduce the radial wavefunction by $u_{n}(r)=\sqrt{4\pi r^{2}}\psi_{n{\bf r}}$. We numerically solve Eq. (85) in a spherical box of radius $R$ by imposing the boundary conditions $u_{n}(r_{0})=u_{n}(R)=0$ with $r_{0}=2200$ and $R=10^{5}$ in the atomic units. We then obtain $N_{b}=80$ energy eigenstates from the lowest-energy solution. Figure 6 shows three typical bound states for the mean-field pseudopotential $V_{\rm mean}=(V_{T}+V_{S})/2$ for the scattering between the electron of the 87Rb$(87s)$ state and the surrounding ground-state 87Rb atoms. The state indicated by the label $\nu=1$ is localized around the outermost antinode of the potential and corresponds to the dominant bound state having the largest overlap with the initial single-particle wavefunction. We note that the contributions from the states with the principal numbers $n=1,2$, which are localized in the inner region and cutoff-dependent, are tiny due to the vanishingly small overlaps with the initial state and thus are neglected in the analysis here.

The initial single-particle wavefunction is the lowest-energy state in the absence of an external potential and, by imposing the boundary conditions above, we obtain its solution as $u_{{\rm ini}}(r)=\sqrt{2/R}\sin(\pi r/R)$ where $u_{{\rm ini}}(r)=\sqrt{4\pi r^{2}}\psi_{{\rm ini},{\bf r}}$. We relate the total number $N$ of particles in the sphere with radial $R$ to the density $\rho$ by using the value of the initial wavefunction at the origin, i.e., $\rho=N|\psi_{{\rm ini},{\bf r}={\bf 0}}|^{2}.$ Finally, we note that the initial coefficients in Eq. (62) are obtained from

[TABLE]

Appendix C Penalty term ensuring the spin conservation

In general, if the Hamiltonian $\hat{H}$ in the original frame satisfies a certain symmetry and has a conserved quantity $\hat{O}$, then the corresponding quantity $\hat{\tilde{O}}=\hat{U}^{\dagger}\hat{O}\hat{U}$ in the transformed frame is also guaranteed to be exactly conserved through the variational time-evolution equations. This follows from the fact that the time-dependent variational principle ensures the symplecticity Haegeman et al. (2011). Yet, in practice, numerical errors accumulated during the integration of highly nonlinear variational equations can lead to an apparent violation of the conservation law. To remedy this, it is useful to add a penalty term as a perturbation on the Hamiltonian so as to make sure that the conservation law is satisfied during the variational time evolution. In particular, in our problem, we ensure the spin conservation $\hat{\sigma}_{e}^{z}+\hat{\sigma}_{{\rm env}}^{z}+\hat{N}=\hat{\sigma}_{e}^{z}+2\hat{N}_{\Uparrow}=1$ by adding the penalty term:

[TABLE]

In the transformed frame, it is given by

[TABLE]

Its expectation value with respect to the Gaussian state can be expressed as

[TABLE]

where $P_{\uparrow}={\rm I}_{2}\otimes((\sigma^{z}+1)/2\otimes{\rm I}_{N_{b}})$ and ${\cal P}=2{\rm Tr}\left[P_{\uparrow}\Omega\right]-\langle\hat{{\rm P}}_{{\rm env}}\rangle_{{\rm GS}}$. Its functional derivative with respect to the mean-field vector $\bm{\phi}$ is given by

[TABLE]

where we introduce the matrix $\Upsilon_{\uparrow}$ as

[TABLE]

Similarly, the functional derivative with respect to the covariance matrix $\Gamma$ can be calculated as

[TABLE]

The perturbations (90) and (95) are added to the functional derivatives (32) and (33), respectively, and ensure the spin conservation during the variational calculations.