\section{INTRODUCTION} \label{Sec.III}
The surface acoustic waves (SAWs) are electro-mechanical vibration waves, which are usually produced and detected by interdigital transducers (IDTs) ~\cite{Datta1986,Wakana2007,Morgan2007}, and propagate along the surface of piezoelectric material. It has $10^5$ times reduction in the propagation speed compared to electromagnetic signals (typically $3000$ m/s for SAW in solid in contrast to $3.0\times10^8$ m/s for electromagnetic wave), or say, it has $10^5$ times shorter wavelengths for the same frequencies of electromagnetic signals. Thus, SAW devices play a very important role in modern electronics including signal processing for radar and mobile telephone systems~\cite{Morgan2007,Kockum2014,Guo2017}. A wide variety of acoustic devices, e.g., delay lines, bandpass filters, micro-resonators, matched filters, and sensors~\cite{Kondalkar2018,Duquesne2019,Le2019}, have been developed and are extensively used in electronic systems. In particular, SAW resonators, which are micro-electromechanical devices, are widespread in the telecommunication and sensor industries.

People have been explored new materials to fabricate SAW devices~\cite{Jin2013,Thalmeier2010} or finding new mechanism to couple the SAW devices with other systems, such that more applications of the SAW devices can be found~\cite{Thevenard2013,Kulin2020} or the SAW devices can more easily act as a platform for lab-on-a-chip systems. For example, the interaction between the acoustic wave and the fluid films~\cite{Rambach2016} could result in many potential applications in biomedicine and biotechnology. The acousto-optic interaction provides wide possibilities for light control and probe of acoustic waves~\cite{Sun2016}. In particular, bulk wave acousto-optic devices have been found numerous applications in light deflectors, modulators, and tunable filters~\cite{White1972,Shi2013,Euchner2012}.

With the rapid progress of quantum science and technology, nanomechanical resonators are recently proposed to be coupled to quantum systems. This might open up new applications of nanomechanical devices to quantum technology~\cite{Kockum2014,Guo2017,Gustafsson2014,Aref2016,Manenti2017}, including quantum signal processing, quantum memories, time delays, quantum transducers, quantum metrology, and single-phonon detection. Superconducting quantum circuit (or superconducting qubit)~\cite{Clarke2008,You2011,Devoret2013,Gu2017} is one of the most promising scalable solid-state quantum information devices~\cite{Liu2017,Rosenberg2017,Dunsworth2018}. The interaction between superconducting quantum circuits and phonons at the quantum level is called as circuit quantum acoustodynamics (cQAD)~\cite{Noguchi2017,Manenti2017}, which is similar to the circuit (cavity) quantum electrodynamics~\cite{Wallraff2004,Blais2004,Gu2017} for the interaction between superconducting quantum circuits (atoms) and microwave photons (photons).

Experimental progress for cQAD has been made in recent years. Quantum ground state and single-phonon control for local vibration quanta in a nanomechanical resonator with microwave frequency were achieved by virtue of a superconducting qubit~\cite{Connell2010}. The coupling between travelling SAW phonons and superconducting quantum circuits was recently observed~\cite{Gustafsson2014}. Furthermore, a superconducting qubit strongly coupled to bulk acoustic wave resonators~\cite{Chu2017}, single-mode SAW resonators~\cite{Bolgar2018,Satzinger2018}, multi-mode SAW resonators~\cite{Lehnert2018}, a phononic crystal~\cite{Bolgar2020} were also demonstrated, respectively. The single-phonon acoustic Stark shifts were observed by coupling a superconducting qubit to a SAW resonator in either weak~\cite{Manenti2017} or strong dispersive regime~\cite{Sletten2019}. The multi-phonon Fock states were created and controlled in bulk acoustic wave resonators by a superconducting qubit~\cite{Chu2018}. Phonon-mediated quantum state transfer and remote qubit entanglement were realized by a SAW resonator~\cite{Bienfait2019}.

Most of superconducting qubits in these cQAD experiments are coupled to either a phononic resonator or a phononic waveguide. In the experiment~\cite{Manenti2017}, a superconducting qubit is coupled to both a SAW resonator and a coplanar waveguide resonator, which is employed for
independent dispersive qubit readout. So far, all of these experiments have been performed around $10\sim 40$ mK, which are typical base temperature of a dilution refrigerator, thus the effect of thermal fluctuations on the coherence of the cQAD system is negligibly small. In other words, thermal fluctuations effect on cQAD system due to high temperatures has not been studied, especially, the quantum-to-classical transition induced by thermal fluctuations in cQAD system remains an open question.

In this work, differently from the experiment~\cite{Manenti2017}, we develop a cQAD system, which consists of a transmon qubit coupled to both a two-dimensional (2D) SAW resonator and a one-dimensional ($1$D) microwave transmission line. The strong coupling between the transmon qubit and the SAW resonator is shown by the energy splitting of the anticrossing in the transmission spectrum through the SAW resonator. We also find that the propagating microwave photons in the transmission line can be controlled by phonons inside the $2$D SAW resonator. In particular, we show the quantum-to-classical transition with the increase of the environmental temperature by observing the variations of the transmission spectra through the interaction system of the transmon qubit and the $2$D SAW resonator. The effective crossover temperature corresponding to the transition regime is estimated by thermal energy and the energy level spacing of the system.

The remainder of this paper is organized as follows: in Sec.~\ref{Sec.setup}, we present the experimental setup and measurement approach. We also experimentally characterize the SAW resonator and the transmon qubit. In Sec.~\ref{model}, we give the theoretical model to describe the interaction system between the transmon qubit and the SAW resonator. In Sec.~\ref{Sec.experiment}, the acoustic Stark shift in the weak dispersive regime is observed and used to calibrate the average phonon number inside the SAW resonator. We show that the propagation microwave photons in the transmission line can be controlled by phonons inside the SAW resonator. Then, we demonstrate that the increase of the effective environmental temperature leads to the quantum-to-classical transition of the cQAD system. In Sec.~\ref{numberical}, we theoretically analyze the experimental results by numerically solving the master equation with parameters given in our cQAD system. Results are finally summarized in Sec.~\ref{conclusion}.

\section{EXPERIMENTAL SETUP AND CHARACTERISTIC FREQUENCIES}\label{Sec.setup}

\begin{figure}[ptb]
\includegraphics[scale=0.08,clip]{Fig1.pdf}
 \caption{(Color online) Circuit quantum acoustodynamics device. (a) Optical micrograph of the device. A transmon (inset) is strongly coupled to the 2D SAW resonator via the interdigital transducer (IDT), and simultaneously coupled to a coplanar waveguide transmission line via a capacitor. BM is an abbreviation of the Brag Mirror.  (b) An equivalent electrical circuit of the device. A SAW resonator is excited and probed via two interdigital transducers (IDTs) connected to ports $1$ and $2$.  The coplanar waveguide transmission line is used to measure the transmon via port $3$. The Brag gratings act as the Brag mirrors.}

\label{Figure1}
\end{figure}

Our device, as shown in Fig.~\ref{Figure1}(a), is fabricated on a quartz substrate. The equivalent electrical circuit in Fig.~\ref{Figure1}(a) is schematically shown in  Fig.~\ref{Figure1}(b).  The cQAD device is formed by a tunable SQUID-based transmon qubit strongly coupled to a 2D SAW resonator via a shunted IDT~\cite{Bolgar2018}. The electrodes of IDT are design to position at the antinodes of the coupled acoustic wave in the resonator. Simultaneously, the transmon qubit is also capacitively coupled to $1$D coplanar waveguide transmission line, which is used to measure the qubit and the cQAD system via the port $3$. A SAW is excited and detected by two identical IDTs connected to ports $1$ and $2$. Two Brag gratings, acting as two Brag mirrors, are used to form a SAW resonator. All IDTs are formed by a periodic array of alternating stripe electrodes and have the same period $p$. The SAW travelling speed in quartz at low temperatures is about $3.16$ km/s. In our device, the periodicity of the IDT stripes is $p=980$ nm and that of the Bragg mirrors is $p/2$, which defines the SAW wavelength $\lambda=p$ and the frequency about $3.2$ GHz.  Our device is mounted at the mixing chamber of a dilution refrigerator with the base temperature $T\approx16.5$ mK. A superconducting Al shield and two $\mu-$metal shields are used to reduce the ambient magnetic field.

Our device is characterized as it is discussed below. We measure the transmission coefficient $t$ through the SAW resonator from the port $1$ to the port $2$ at the base temperature, and then obtain frequencies of the SAW resonator. Figure~\ref{Figure2}(a) shows the normalized transmission coefficient $|t|^2$ with the probing power lower than $-134$ dBm, which corresponds to the average phonon number around $1$ inside the SAW resonator. There are three resonances at the frequencies $\omega_{1}=2\pi\times3.144$ GHz, $\omega_r=2\pi\times3.162$ GHz, and $\omega_2=2\pi\times3.184$ GHz, respectively. In the experiment, we use acoustic Stark shifts, in contrast to the ac-Stark shift in the system of the circuit QED~\cite{Schuster2005,Ong2011},  to calibrate the phonon number inside the resonator, which will be further discussed later in the paper. The decay rate of the resonator mode at the frequency $\omega_r=2\pi\times3.162$ GHz is $\kappa=2\pi\times1.56$ MHz by extracting the linewidth from $|t|^2$.

The frequency of the qubit can be characterized by measuring the reflection coefficient $r$ of the microwave signal through the waveguide via the port $3$ using a vector network analyzer (VNA) at the base temperature. Figure~\ref{Figure2}(b) shows a $2$D plot of the normalized reflection coefficient $|r|$ in the frequency range $3.0\sim5.4$ GHz with the magnetic flux bias $\delta\Phi$ from $-400 \,\text m\Phi_0$ to $400\,\text m\Phi_0$, where $\Phi_0$ is the flux quantum. The transition from the ground state to the first excited state of the qubit is revealed as dips. The transition frequency $\omega_a$ at $\delta\Phi=0$ reaches the maximum value $\omega_a=2\pi\times5.260$ GHz and the extracted relaxation rate at this point is $\Gamma=2\pi\times10.48$ MHz, which is similar to that in Ref.~\cite{Bolgar2018}. By fitting the spectroscopy of the qubit with the formula $\omega_a=\sqrt{8E_{\rm C}E_{\rm J}(\Phi)}-E_{\rm C}$, we obtain the qubit charging energy $E_{\rm C}\approx2\pi\times160$ MHz, and the Josephson energy $E_{\rm J}(\Phi)=E_{\rm J,max}|\cos(\pi\Phi)/\Phi_0|$, in which $E_{\rm J,max}\approx2\pi\times22.88$ GHz. In experiment, we can adjust the transition frequency of the qubit by changing the bias flux $\Phi$, injected to the SQUID loop of the qubit.

 \begin{figure}[ptb]
	\includegraphics[scale=0.25,clip]{Fig2.pdf}
		\caption{(Color online) (a) Normalized transmission through the SAW resonator. The transmon is coupled to the mode $\omega_r=2\pi\times3.162\,$GHz and decoupled to another two modes $\omega_{1}=2\pi\times3.144\,$GHz and $\omega_{2}=2\pi\times3.184\,$GHz. (b) Reflection spectrum of the transmon qubit vs biased flux plotted as a reflection coefficient $|r|$.}
	\label{Figure2}
\end{figure}

\section{THEORETICAL MODEL}\label{model}

Our cQAD system is formed by a transmon qubit interacting with a single-mode SAW resonator, thus this system can be modeled as a two-level system interacting with a single-mode cavity, given in the cavity quantum electrodynamics system.  Below, the transmon qubit is called as the qubit for convenience. The single-mode Hamiltonian of the SAW resonator~\cite{Noguchi2017,Manenti2016} can be written as $H/\hbar=\omega_r({\hat{a}^\dagger}\hat{a}+\frac{1}{2})$, where $\hat{a}^\dagger(\hat{a})$ is the bosonic creation (annihilation) operator of the single SAW resonator mode and $\hbar$ is the reduced Plank constant. In our experiment, the mode of the SAW resonator with the frequency $\omega_r=2\pi\times3.162$ GHz  is coupled to the qubit. The Hamiltonian of the qubit can be written as $H/\hbar=\omega_a\sigma_z/2$ with the transition frequency $\omega_a$ of the qubit, $\sigma_z=|e\rangle\langle e|-|g\rangle\langle g|$ is the Pauli $Z$ operator with the ground state $|g\rangle$ and the excited state $|e\rangle$ of the qubit. The interaction Hamiltonian between the SAW resonator and the qubit is $H_{\rm int}/\hbar=g\sigma_x({\hat{a}^\dagger}+\hat{a})$ with the coupling strength $g$. Under the rotating-wave approximation, the total Hamiltonian is given as below
\begin{equation}
H/\hbar=\omega_r{\hat{a}^\dagger}\hat{a}+\frac{1}{2}\omega_a\sigma_z+g(\sigma_{+}\hat{a}+\sigma_{-}\hat{a}^\dagger),\label{Eq:1}
\end{equation}
where $\sigma_{+}=|e\rangle\langle g| \,(\sigma_{-}=|g\rangle\langle e|)$ is raising (lowering) operator of the transmon qubit. Hereafter, the state $|e,n\rangle$ ($|g,n\rangle$) denote that the qubit is in the excited (ground) state and there are $n$ phonons in the SAW resonator. The ground state $|g,0\rangle$ of the system means that the qubit is in the ground state and the SAW resonator is in the vacuum state, this state can be prepared in our experiment by cooling the system down to the base temperature of the dilution refrigerator.

The coupling between the qubit and the SAW resonator mixes the states $|g,n+1\rangle$ and $|e,n\rangle$ as two dressed states
\begin{eqnarray}
   	|n,+\rangle&=\cos\theta_{n}|g,n\rangle + \sin\theta_{n}|e,n-1\rangle, \\
	|n,-\rangle&=\sin\theta_{n}|g,n\rangle - \cos\theta_{n}|e,n-1\rangle,
\end{eqnarray}
as schematically shown in Fig.~\ref{Figure3}(a) for general case, and the corresponding eigenvalues write
\begin{equation}
E_{n,\pm}=\omega_r\left(n-\frac{1}{2}\right)\pm \frac{1}{2}\sqrt{\Delta^2+4g^2n},
\end{equation}
with $n\geq 1$, $\Delta=\omega_a-\omega_r$ and $\theta_{n}$ determined by
\begin{equation}
	\tan \left(2\theta_{n}\right)=\frac{2g\sqrt{n}}{\Delta}.
\end{equation}

\begin{figure}[ptb]
	\includegraphics[scale=0.12,clip]{Fig3.pdf}\
	\caption{(Color online) Energy spectrum of the uncoupled (left) and coupled  qubit-resonator system (right) under the nonresonant condition (a) and resonant condition (b). The degeneracy of the states $|n,e\rangle$ and $|n+1,g\rangle$ in (b) is lifted $2g\sqrt{n+1}$ by the coupling. (c) Energy diagram of the qubit (dashed blue curve $|e,0\rangle$), the SAW resonator (dashed orange curve $|g,1\rangle$), the coupled system (solid purple and yellow curve $|1,\pm\rangle$), parameters we choose here are $E_{\rm C}=2\pi \times 160$MHz, $E_{\rm J,max}=2\pi \times 22.88$GHz, $g=2\pi \times 100$MHz and $\omega_{r}=2\pi \times 4$GHz. (d) Schematic diagram for transitions corresponding to different frequencies as shown in (b).}
	\label{Figure3}
\end{figure}

When the qubit and the SAW resonator are not resonant, and in the larger detuning regime, i.e., $|\Delta=\omega_a-\omega_r|\gg g$, the Hamiltonian in Eq.~(\ref{Eq:1}) can be equivalently written
\begin{equation}
H_{\rm eff}/\hbar=\left(\omega_r+\frac{g^2}{\Delta}\sigma_z\right)\hat{a}^\dagger \hat{a}+\frac{1}{2}\left(\omega_a+\frac{g^2}{\Delta}\right)\sigma_z.
\end{equation}
In this nonresonant case, the dressed qubit level separation is given by $\widetilde{\omega_a}=\omega_a+2ng^2/\Delta+g^2/\Delta$, which depends on the number of phonons $n=\langle \hat{a}^\dagger \hat{a} \rangle$ inside the SAW resonator. We call the terms $2ng^2/\Delta$ and $g^2/\Delta$ as the acoustic Stark and Lamb shifts of the qubit frequency, in contrast to the ac Stark and Lamb shifts in the system of the cavity quantum electrodynamics, and  use Stark shift to calibrate the phonon number inside the SAW resonator induced by driving field in the experiment.

When the transmon qubit resonantly interacts with the SAW resonator, i.e., $\omega_{a}=\omega_{r}=\omega$ and $\Delta=0$ , the energy spacing becomes $\Delta_{n}=2g\sqrt{n}$, corresponding to two dressed states
\begin{equation}
|n,\pm\rangle=\frac{1}{\sqrt{2}}(|g,n\rangle\pm|e,n-1\rangle),
\end{equation}
as schematically shown in Fig.~\ref{Figure3}(b).

If a probe field is applied to the SAW resonator, in which the interaction Hamiltonian between the resonator and the driving field is proportional to $a^{\dagger}+a$, then transitions from states $|1,\pm\rangle$ to the state $|g, 0\rangle$ as well as from the states $|n+1,\pm\rangle$ to the states $|n,\pm\rangle$ can be induced. The transitions between the state $|g,0\rangle$ and the states $|1,\pm\rangle$ with frequencies $\omega+g$ and $\omega-g$ lead to the standard vacuum Rabi splitting, which can be observed with two peaks, separated by $2g$, in the absorption spectrum. This vacuum Rabi splitting corresponds to an anticrossing in the energy diagram of the coupled qubit-resonator system as shown in Fig.~\ref{Figure3}(c).

For the upper levels, there are four transitions from the states $|n+1,\pm\rangle$ to the states $|n,\pm \rangle$ for $n\geq 1$ corresponding to four frequencies
\begin{eqnarray}
\omega_{1,n}&=&\omega+g\left(\sqrt{n+1}-\sqrt{n}\right),\label{Eq:7}\\
\omega_{2,n}&=&\omega-g\left(\sqrt{n+1}-\sqrt{n}\right),\label{Eq:8}\\
\omega_{3,n}&=&\omega+g\left(\sqrt{n+1}+\sqrt{n}\right),\label{Eq:9}\\
\omega_{4,n}&=&\omega-g\left(\sqrt{n+1}+\sqrt{n}\right),\label{Eq:10}
\end{eqnarray}
here, the subscript $``i, n"$ means the $i$th transition from $|n+1,\pm\rangle$ to $|n,\pm \rangle$. $\omega_{1,n}$ and $\omega_{3,n}$ are frequencies from the state $|n+1,+\rangle$ to the states $|n,+\rangle$ and $|n,-\rangle$, respectively. However  $\omega_{2,n}$ and $\omega_{4,n}$ are frequencies from the state $|n+1,-\rangle$ to the states $|n,-\rangle$ and $|n,+\rangle$, respectively.
Thus, as schematically shown in Fig.~\ref{Figure3}(d), in the quantum case, many peaks appear at the positions $\omega\pm g$, $\omega\pm g(\sqrt{n+1}-\sqrt{n})$, and $\omega\pm g(\sqrt{n+1}+\sqrt{n})$ for all $n$ with $n\geq 1$ in the spectrum. If $n$ is a very large number such that
\begin{equation}
\omega_{2,n}-\omega_{1,n}=2g(\sqrt{n+1}-\sqrt{n})\approx \frac{g}{\sqrt{n}}\approx 0,
\end{equation}
then $\omega_{1,n}\approx \omega_{2,n}=\omega$, that is, the two peaks for frequencies $\omega_{1,n}$ and $\omega_{2,n}$ in the spectrum merge into one in the limit of the large $n$. It is not difficult to find that the separation between two successive peaks, e.g., corresponding to frequencies $\omega_{i, n+1}$ and $\omega_{i,n}$ with $i=1,\,2,\,3,\,4$, is of the order of $g/\sqrt{n}$. Thus, many peaks are not resolved in the limit of the large phonon number $n$ inside the resonator. Moreover, we can also find that the transition elements from the state $|n+1,+ \rangle$ to the state $|n, +\rangle$ or from the state $|n+1,- \rangle$ to the state $|n, -\rangle$ are proportional to $(\sqrt{n+1}+\sqrt{n})$, but the transition elements from the state $|n+1,- \rangle$ to the state $|n, +\rangle$ or from the state $|n+1,+\rangle$ to the state $|n, -\rangle$ are proportional to $(\sqrt{n+1}-\sqrt{n})$. Roughly speaking, the heights of the peaks corresponding to the frequencies $\omega_{1,n}$ and $\omega_{2,n}$ are higher than those of peaks corresponding to the frequencies $\omega_{3,n}$ and $\omega_{4,n}$ in the spectrum. Thus, in the limit of the large $n$, peaks corresponding to the frequencies $\omega_{3,n}$ and $\omega_{4,n}$ in the spectrum are not observable, and peaks corresponding to the frequencies $\omega_{1,n}$ and $\omega_{2,n}$ in the spectrum emerge into one, the system is changed from the quantum to the classical regime.

In practice, the environmental temperature of the system can be used to induce quantum-to-classical transition. For example, in our experiments, the mean thermal phonon numbers $\langle n_{\rm th}\rangle=[\exp(\hbar\omega_{r}/k_{B}T)-1]^{-1}$ for the resonator and $\langle n_{\rm th}^{\prime}\rangle=[\exp(\hbar\omega_{a}/k_{B}T)-1]^{-1}$ for the qubit are  about $\langle n_{\rm th}\rangle\approx \langle n_{\rm th}^{\prime}\rangle \approx 0$ with $\omega_r=2\pi\times3.162$ GHz  and $\omega_a=2\pi\times3\sim2\pi\times5$ GHz at the base temperature $T=16.5$ mK, and thus the cQAD system can be prepared to the ground state $|g,0 \rangle$,  the transitions mostly occur between the state $|g, 0 \rangle$ and the states $|1,\pm\rangle$ at the base temperature. In this case, transition spectra of dressed states $|1,\pm\rangle$ to the state $|g, 0 \rangle$, with two peaks separated by $2g$, can be observed. However, with the increase of the temperature, the mean thermal phonon number $\langle n_{\rm th}\rangle$ and $\langle n_{\rm th}^{\prime}\rangle$  are also increased, then more high energy levels are excited, and the spectra with multiple peak structure should be observed at the intermediate temperatures. If the temperature is further increased such that the number of the thermal phonons is very large, multiple peaks merge into two peaks, and finally one peak. This is because: (i) the lower energy levels are all saturated at the high temperature and transition can only occur from large quantum number, whose transition frequencies are densely close to the resonator frequency, as discussed above, (ii) the spectrum width of each peak becomes broader and broader with the increase of the temperature,  thus many peaks overlapped into one. The temperature effect on our experimental results will be carefully analyzed and discussed in following sections.

\section{EXPERIMENTAL RESULTS}\label{Sec.experiment}

At the base temperature $T\approx 16.5$ mK of the dilution refrigerator, the mean thermal excitation phonon numbers for both the acoustic mode frequency $\omega_{r}$  and the qubit frequency are about zero, as discussed above. Thus, the system works in the quantum regime because the effect of thermal fluctuations on the coupled system between the SAW resonator and the qubit is negligibly small and the quantum fluctuations are dominant. In our experiments, the qubit is designed to be coupled to the SAW mode with the frequency $\omega_r=2\pi\times3.162$ GHz by matching $p$ with the wavelength of the frequency $\omega_{r}$,  and decoupled from other two modes with frequencies $\omega_1=2\pi\times 3.144$ GHz and $\omega_2=2\pi\times3.184$ GHz. Note that the coupling between the SAW resonator and the qubit is usually characterized in two ways. (i) When the qubit and the SAW resonator are in the regime of large detuning, the coupling is measured by the acoustic Stark shift of the qubit frequency or the dispersive shift of the SAW resonator frequency induced by the qubit. (ii) When the qubit resonantly interacts with the SAW resonator, the coupling strength can be found by measurement of the anticrossing.

\subsection{Non-resonant case}
We first study the large detuning case. In Refs.~\cite{Manenti2017} and~\cite{Sletten2019}, the single-phonon acoustic Stark shift in the weak and the strong dispersive regimes are observed, respectively. However, the phonon number dependence of the acoustic Stark shift from a few to large number of phonons is still missing, because it was not easy to clearly detect large qubit frequency shift with limited measurement bandwidth in the SAW resonator in previous cQAD systems. In our experiment, the qubit is coupled to the transmission line and the limitation of the measurement bandwidth is overcome. Thus, it is not difficult to detect the phonon number dependence of the acoustic Stark shift in the large phonon number range. We apply the magnetic flux bias such that the qubit works at the frequency $\omega_a=2\pi\times4.442$ GHz. In this case, the detuning between the qubit and the SAW resonator is $\Delta=2\pi\times1.282$ GHz, which is much larger than $g$. Then we measure the phonon-dependent acoustic Stark shift via the reflection spectrum of the qubit through the port $3$ using VNA. The acoustic phonon number is changed by adjusting the driving power on the SAW resonator. As shown in Fig.~\ref{Figure4}(a), we measure the phonon dependent acoustic Stark shift in the range from $-143$ dBm to $-103$ dBm, which corresponds to the average phonon number from $\langle n\rangle\sim0.1$ to $\langle n\rangle\sim1000$ inside the SAW resonator. Here, we define the phonon number operator $n=a^{\dagger} a$.  The acoustic Stark shift induced by one phonon is extracted to be $2g^2/\Delta\approx510$ KHz via the experimental data in Fig.~\ref{Figure4}(b), which means the device works in the weak dispersive regime. With the power calibration, we find that the microwave to phonon conversion efficiency of each IDT is $-8\pm3$ dB. To keep low phonon number excitation, the probing power used in following experiment is set lower than $-134$ dBm, which corresponds to the average phonon number $\langle n\rangle\sim1$ inside the SAW resonator. We note that although the qubit is weakly coupled to the $1$D transmission line, we can clearly resolve the reflection spectrum of the propagation microwave field for the average phonon number inside the SAW resonator down to $\langle n\rangle\sim 21$ because of $21\times2g^2/\Delta\ >\Gamma,\kappa$. This means that the reflection of the microwave field could be controlled by the phonons inside the SAW resonator, which is similar with the demonstration in Ref.\cite{Peng2018} that reflection of the microwave field controlled by the microwave photons inside a microwave transmission line resonator. Here, our result is the first step towards the microwave photon control by phonons in a quantum hybrid system on-chip.

 \begin{figure}[ptb]
	\includegraphics[scale=0.25,clip]{Fig4.pdf}

	\caption{(Color online) (a) Spectrum of the transmon vs driving power of the SAW resonator as a reflection coefficient $|r|$. $-134\,$dBm corresponds to average phonon number $\langle a^{\dagger}a\rangle=\langle n\rangle=1$ inside the SAW resonator. (b) The spectrum with theoretically calculated transition frequency (dashed line) from ac-Stark shift.}
	\label{Figure4}
\end{figure}

 \begin{figure*}[ptb]
	\includegraphics[scale=0.13,clip]{Fig5.pdf}\
	\caption{(Color online) (a), (b), and (c) Normalized transmission amplitude around the SAW mode($\omega_r=2\pi\times3.162\,$GHz) as a function of $\delta\phi$ under different temperatures. The temperature of the system are $16.5\,$mK, $50.5\,$mk, $98.5\,$mk, respectively, from (a) to (c). (d), (e), and (f) Phases of the transmission spectra under different temperatures from the experimental data, the temperatures of the system are $16.5\,$mK, $50.5\,$mK, $98.5\,$mK, respectively, from (d) to (f). (g) The calculated spectra around $\omega_a=\omega_r$ with the model of the lowest five energy levels of transmon (the model of the lowest two levels of transmon), plotted with solid curves (dashed curves) for A(A$'$), B(B$'$), and C(C$'$) with experimental parameters $\kappa=2\pi\times1.56\,$MHz and $\gamma=2\pi\times7.90\,$MHz, and corresponding temperatures are  $16.5\,$mK, $50.5\,$mK, $98.5\,$mK, respectively. (h) The calculated spectra around $\omega_a=\omega_r$ with the model of the lowest five energy levels of transmon(the model of the lowest two levels of transmon) plotted with solid (dashed) curves A(A$'$), B(B$'$), and C(C$'$) with theoretically chosen parameters $\kappa=2\pi\times0.4\,$MHz and $\gamma=2\pi\times6\,$MHz, and the corresponding temperatures are $16.5\,$mK, $50.5\,$mK, $98.5\,$mK, respectively.}
	\label{Figure5}
\end{figure*}

\subsection{Resonant case}

In the resonant case, the interaction between the qubit and the SAW resonator can be studied via the transmission coefficient $t$ through the SAW resonator around the frequency $\omega_r=2\pi\times 3.162$ GHz using the VNA by varying magnetic flux bias $\delta\Phi$. As shown in Fig.~\ref{Figure5}(a) for the amplitudes of normalized transmission spectra measured at the base temperature $T\approx 16.5$ mK, in which the thermal phonon numbers for the qubit and resonator frequecies are $\langle n_{\rm th}\rangle \approx \langle n_{th}^{\prime}\rangle \approx0$, two hybridized modes of the qubit and the resonator are repelled by each other, and the anticrossing due to the vacuum Rabi splitting is clearly observed when $\omega_a=\omega_r$.   The coupling strength $g=2\pi\times18.5$ MHz is extracted by fitting the anticrossing, which is closed to the value obtained in Ref.~\cite{Bolgar2018}. $2g$ is just a frequency difference between two transition frequencies from the dressed states $|1,\pm\rangle$ to the ground state $|0,g\rangle$, and is the distance between two peaks in the transmission spectrum as shown in the curve A (or A$^{\prime}$ calculated using the model of the lowest two energy levels) calculated using the model of the lowest five energy levels of Fig.~\ref{Figure5}(g) for $\omega_{a}=\omega_{r}$, which is simulated using experimental parameters and the master equation  (see, Subsec. B of Sec.~\ref{numberical}).  The corresponding phase spectra from experimental data are presented in Fig.~\ref{Figure5}(d), in which the vacuum Rabi splitting is also observed. This result indicates that the strong coupling between the qubit and the SAW resonator is realized, i.e., $g>\Gamma, \,\kappa$.

The vacuum Rabi splitting can also be considered as the normal mode splitting for two coupled systems using the semiclassical theory. To further explore the quantum nature for the coupled system of the qubit and the SAW resonator, the phonon number dependent coupling should be observed, thus we slightly increase the temperature to, e.g., $T=50.5$  mK and $T=98.5$ mK, which correspond to the thermal phonon number $\langle n_{\rm th}\rangle=\langle n_{\rm th}^{\prime}\rangle\approx0.04 $ and  $\langle n_{\rm th}\rangle=\langle n_{\rm th}^{\prime}\rangle\approx 0.17 $ around the resonant frequency, respectively. In this case, the thermal phonon number is very small and the system still works in the quantum regime, but the temperature induced transitions from high energy levels are also involved, e.g., four lower transitions from the states $|2,\pm \rangle$ to the states $|1,\pm \rangle$. In Figs.~\ref{Figure5}(b), (c), (e) and (f), the amplitudes and the corresponding phases of the transmission spectra are measured at two temperatures.

Using experimental parameters in Figs.~\ref{Figure5}(b) and (c), the spectra corresponding to $\omega_{a}=\omega_{r}$ are simulated by the master equation (see, Subsec. B of Sec.~\ref{numberical}) for the  transmon model with the lowest five energy levels because of the weak anharmonicity and shown in solid curves B and C of Fig.~\ref{Figure5}(g) for these two temperatures. As comparison, the spectra corresponding to the transmon model with the lowest two energy levels are plotted as dashed curves B$'$ and C$'$ in Fig.~\ref{Figure5}(g). In curve C$'$, we find two more peaks appear in the middle besides the two peaks corresponding to the vacuum Rabi splitting. These two peaks correspond to the transitions from the state $|2,+\rangle$ to the state $|1,+\rangle$  and from the state $|2,-\rangle$ to the state  $|1,-\rangle$. These additional peaks show a quantum nature of the coupled system of the qubit and the SAW resonator. However, these additional peaks are difficult to distinguish in curve C because the quantum nature of the coupled system is covered by the weak anharmonicity of the transmon at high temperature. In order to further show the dissipation effect on the spectra, we present the calculated spectra in Fig.~\ref{Figure5}(h) around $\omega_{a}=\omega_{r}$ by theoretically choosing smaller dissipation parameters of the system, e.g., $\kappa=2\pi\times0.4\,$MHz and $\gamma=2\pi\times6\,$MHz. As shown in curve $C^{\prime}$ in Fig.~\ref{Figure5}(h), the multi-peak structure becomes more clear in the case of the small dissipation and large anharmonicity (two energy level model).  In our experiments, the conversion efficiency from microwave photon to phonon through the IDT is low and the transmon is weakly coupled to the transmission line. It is difficult to monitor the additional peaks either through SAW resonator or transmission line in experiments even that the dissipation rates become small. We expect the these additional peaks can be experimentally observed in an improved architecture with small dissipations,  a strong anharmonicity (e.g., a flux qubit\cite{Peng2018}), high conversion efficiency for IDT.

\section{Temperature induced quantum-to-classical transitions: Experimental and numerical results}\label{numberical}

In this section, we will systematically study the temperature effect on the cQAD system when the qubit resonantly interacts with the SAW resonator. In particular, we will study the quantum-to-classical transition induced by the temperature. To show a clear comparison for the effect of different temperatures on the quantum nature of the system, we regulate temperatures of the mixing chamber of the refrigerator such that different thermal phonon occupations for the SAW resonator and the qubit can be induced.

\subsection{Experimental results and qualitative analysis}

 In our experiments, we measure the transmission spectrum through the SAW resonator around the frequency $\omega_r$ by varying the temperature from $16.5$ mK to $349$ mK. Figures~\ref{Figure5}(a)-(c) show the experimental results for the temperatures $16.5$ mK, $50.5$ mK and $98.5$ mK.   Figures~\ref{Figure6}(a), (c), (e) and (g) show the experimental results when the coupled system is in different temperatures, i.e., $149$ mK, $200$ mK, $254$ mK, and $349$ mK. Figures~\ref{Figure6}(b), (d), (f) and (h) are simulation results using experimental parameters. Figure~\ref{Figure6}(i) is transmission spectra corresponding to Figs.~\ref{Figure6}(b), (d), (f) and (h) when
the qubit resonantly interacts with the resonator, i.e., $\omega_{a}=\omega_{r}$. The resonant transmission spectra corresponding to Figs.~\ref{Figure5}(a)-(c) are also plotted in Fig.~\ref{Figure6}(i) for comparison. As shown in Fig.~\ref{Figure5} and Fig.~\ref{Figure6}, with the increase of the temperature $T$, the thermal fluctuations start to play the role, two hybridized modes bend toward each other, the level repelling to the level attraction becomes obvious. When the temperature is close to $T_{c}\sim\hbar\omega_{a}/k_{B}\approx 149$ mK, with the Boltzmann constant $k_{B}$, as shown in Fig.~\ref{Figure6}(a), the two hybridized modes almost coalesce at the point of the resonant frequency $\omega_a=\omega_r=2\pi\times3.162$ GHz, the Rabi splitting almost disappears. When $T\gg T_{c}$, the thermal fluctuation of the system is dominant. The lower energy levels of the system are saturated and the energy level splitting completely disappears, and the transmission spectrum, as shown in Fig.~\ref{Figure6}(g) and the curve corresponding to the temperature $349$ mK in Fig.~\ref{Figure6}(i),  resembles that of the harmonic oscillator, and has only one peak located at the frequency $\omega_r=2\pi\times3.162$ GHz. We also measure the phases corresponding to the amplitudes of transmission spectra through the SAW resonator, the crossover from the anticrossing to the level attraction can also be observed.

In fact, our experiments show a temperature induced transition from the quantum to the classical regime in the cQAD system. In the lower temperature limit, e.g., the base temperature of the dilution refrigerator,  the thermal phonon number is near zero. As shown in Fig.~\ref{Figure5}(a), we observe a clear vacuum Rabi splitting, which corresponds to the distance between two peaks in curve A calculated using the transmon model of the lowest five energy levels or curve A$^{\prime}$ calculated using the transmon model of the lowest two energy levels of Fig.~\ref{Figure5}(g). With the increase of the temperature, the thermal excitations in the system are increased and more peaks should appear in the spectrum.  Although we can not observe multiple peaks in current experiment due to the limitation of sample parameters and measurement techniques, but as shown in the curve C$'$ in  either Fig.~\ref{Figure5}(g) or Fig.~\ref{Figure5} (h), two more peaks for the transitions from the state $|2,+\rangle$ to the state $|1,+\rangle$  and from  the state $|2,-\rangle$ to the state  $|1,-\rangle$ could appear. This shows a quantum nature of the system. If the temperature is further increased, and the thermal excitation number is far larger than one, all of the lower energy levels are occupied. The transitions occur between higher energy levels with the large phonon number $n$. In this case, the frequency differences between $\omega_{1,n}$ in Eq.~(\ref{Eq:7}) and $\omega_{2,n}$ in Eq.~(\ref{Eq:8}), which are proportional to $2g(\sqrt{n+1}-\sqrt{n})$,  becomes smaller than the decay rates of the acoustic resonator and qubit, i.e., $2g(\sqrt{n+1}-\sqrt{n}) << \Gamma,\,\kappa$, thus the peaks corresponding to these frequencies are not resolved and merge into one around the resonant frequency $\omega$, the height of this central peak is proportional to $\sqrt{n}$. Because the heights of the peaks corresponding to the frequencies $\omega_{3,n}$ and $\omega_{4,n}$ in the spectrum are proportional to $2g(\sqrt{n+1}-\sqrt{n})$, thus they are not observable in the limit of the large $n$. That is, all quantum nature disappear in the high temperature, the classical behavior is observed.  The numerical simulations in Fig.~\ref{Figure6}(i) clearly show the spectrum variations when the coupled system is changed from the quantum to classical regime by varying the environmental temperature. The change from the level splitting to the level attraction around $T_{c}=149$ mK indicates that the cross-over temperature of quantum-to-classical transition in our cQAD is about $T_{c}=149$ mK. When the temperature $T\gg149\,$mK, we observe the classical behavior of the system.

 \begin{figure*}[ptb]
	\includegraphics[scale=0.17,clip]{Fig6.pdf}\

	\caption{(Color online) Normalized transmission amplitude around the SAW mode($\omega_r=2\pi\times3.162\,$GHz) as a function of $\delta\phi$ under different temperatures. (a),(c),(e) and (g) are the experimental results at $149$ mK, $200$ mk, $254$ mk, $349$ mk, respectively. The corresponding numerical simulations are shown in (b),(d),(f) and (h), respectively. (i) The calculated spectra when the SAW resonator resonantly interacts with the qubit, i.e. $\omega_{r}=\omega_{a}$, for the temperatures $16.5\,$mK, $50.5\,$mK, $98.5\,$mK, $149\,$mK, $200\,$mK, $254\,$mK and $349\,$mK, respectively.}
	\label{Figure6}
\end{figure*}

\subsection{Theoretical method and numerical simulations}
To further understand the quantum-to-classical transition induced by the temperature, we theoretically analyze and simulate the transmission spectra in Figs.~\ref{Figure6} (a), (c), (e) and (g) via the master equation
\begin{equation}
\begin{split}
\frac{\partial\rho}{\partial t}=&-i[\mathcal{H}+H_p,\rho]+\frac{\kappa}{2}\left[n_{\rm th}(\omega_r)+1\right]D[\hat{a}]\rho \\
&+\frac{\kappa}{2}n_{\rm th}(\omega_r)D[\hat{a}^\dagger]\rho \\
&+\frac{\Gamma}{2}\left[n_{\rm th}^{\prime}(\omega_a)+1\right]D\bigg[\sum_{l=e,f,\cdots}\frac{g_{l-1,l}}{g_{g,e}}\sigma_{l-1,l}\bigg]\rho \\
&+\frac{\Gamma}{2}n_{\rm th}^{\prime}(\omega_a)D\bigg[\sum_{l=e,f,\cdots}\frac{g_{l-1,l}}{g_{g,e}}\sigma_{l-1,l}^\dagger\bigg]\rho,
\end{split} \label{eq12}
\end{equation}
here $\mathcal{H}$ is ladder-model Hamiltonian with $l$ transmon levels~\cite{Koch2007}. In our theoretical study,  we use the Hamiltonian $\mathcal{H}$ with the lowest five energy levels of transmon for numerical simulation and the lowest two energy levels for theoretical analysis and comparison with the five energy levels, $\sigma_{l-1,l}=|l-1\rangle\langle l|$ and $\rho$  is the density matrix operator of the qubit-resonator system, $H_p=\varepsilon(\hat{a}^{\dagger}+\hat{a})/2$ is the interaction Hamiltonian between the probe field and the SAW resonator with the interaction strength $\varepsilon$, and $D[\hat{o}]\rho=2\hat{o}\rho{\hat{o}^\dagger}-{\hat{o}^\dagger}\hat{o}\rho-\rho{\hat{o}^\dagger}\hat{o}$ denotes the Lindblad type of the dissipation. $\kappa$ and $\Gamma$ are damping rates corresponding to the resonator and the transmon qubit, respectively. Noting that the temperature dependence of $\Gamma$ on piezoelectric substrate is ignored because it makes small difference to the results up to $349\,$mK in our experiment. In fact, the loss channels of qubits on piezoelectricity substrate are expected to be more complicate than qubits on Sapphire or Si~\cite{Lisenfeld2007,Sun2012,Diniz2020,Singh2020}. The related discussions are beyond the topic of this paper. The mean thermal phonon numbers for the resonator and the qubit are $n_{\rm th}(\omega_r)$ and $n_{\rm th}^{\prime}(\omega_a)$, respectively, which are calculated using formula of the Bose-Einstein distribution. The transmission rate, calculated analytically and numerically, is defined as
\begin{equation}
	t=\frac{i\kappa\langle{\hat{a}}\rangle_{ss}}{\varepsilon}.
\end{equation}

Before going to the numerical simulations, we first give a simple theoretical analysis. Using Hamiltonian $\mathcal{H}$ with the lowest two levels of transmon is enough for analysing characteristics of spectra. Under the secular approximation~\cite{Tian1992}, and with the detailed calculation steps given in the supplemental material~\cite{supplemental_material}, the rescaled transmission rate is given as
\begin{eqnarray}
t\cdot\bigg(\frac{i \kappa}{\varepsilon}\bigg)^{-1}&=&\langle \hat{a} \rangle_{ss}=\frac{1}{\sqrt{2}}\sum_{\eta=\pm}\rho_{1,\eta;0}^{(1)}\\
&+&\sum_{\eta=\pm,\xi=\pm}\sum_{n=1}^{\infty} \frac{C^{n}_{\eta,\xi}}{2} \rho_{n+1,\eta;n,\xi}^{(1)}, \nonumber
\end{eqnarray}
with $C_{\eta,\xi}^n=\sqrt{n+1}+\epsilon_{\eta}\epsilon_{\xi}\sqrt{n}$ ($\epsilon_{\eta}=\pm1,\epsilon_{\xi}=\pm1$). The detailed expressions of $\rho_{1,\eta;0}^{(1)}$ and $\rho_{n+1,\eta;n,\xi}^{(1)}$, which are temperature-dependent, are given in supplemental material~\cite{supplemental_material}. The term of $\rho_{1,\eta;0}^{(1)}$ corresponds to the vacuum Rabi splitting. For given $n$, the relative strength of $C_{\pm,\pm}^n$ shows that amplitudes corresponding to transition frequencies $\omega_{1,n}$ and $\omega_{2,n}$ are larger than those corresponding to frequencies $\omega_{3,n}$ and $\omega_{4,n}$, and this results in that the peaks corresponding to $\omega_{3,n}$ and $\omega_{4,n}$ are difficult to be observed in the spectrum. In particular, the peaks corresponding to $\omega_{3,n}$ and $\omega_{4,n}$ are not observable in the limit of large $n$.  This result coincides with the conclusion discussed in Section \uppercase\expandafter{\romannumeral3}. For given $\eta$ and $\xi$, the relative strengths of $\rho_{1,\eta;0}^{(1)}$ and $C^{n}_{\eta,\xi} \rho_{n+1,\eta;n,\xi}^{(1)}$ corresponding to the heights of peaks in the spectrum varies when the environmental temperature $T$ is changed. In the lower temperature with the very smaller $n_{\rm th}$ and $n_{\rm th}^{\prime}$, the most obvious peaks correspond to those of vacuum Rabi splitting. While for the higher temperature with the larger $n_{\rm th}$ and $n_{\rm th}^{\prime}$, the most dominated peaks correspond to those of both vacuum Rabi splitting and the transitions from higher energy levels. When $n_{\rm th}$ and $n_{\rm th}^{\prime}$ become very large, all peaks merge into one peak, as shown in experimental result. We note that the result derived from the secular approximation does not agree well with the experiment, this is because all oscillatory terms are neglected as shown in supplemental material~\cite{supplemental_material}. These neglected terms have small effect on the result under the condition of the lower temperature. However, they take great effect on the result under the condition of the high temperature. This is because the effective dissipation becomes comparable or larger than the coupling strength between the SAW resonator and the qubit with the increase of the environmental temperature. Nevertheless, it is helpful enough to clearly show the variation trend of peaks in the spectrum with the increase of the temperature. More accurate result is shown below by numerical simulation.

Simulation on the transmission spectra is completed by numerically solving master equation in Eq.~(\ref{eq12}) under different environmental temperature~\cite{Rau2004}. The simulation results are shown in curves in Fig.~\ref{Figure5}(g), Figs.~\ref{Figure6}(b), (d), (f), (h) and Fig.~\ref{Figure6}(i), respectively. Figure~\ref{Figure5}(g) and Fig.~\ref{Figure6}(i) are
resonant transmission spectra at the different temperature.  We find that the experimental results at different temperatures can be well explained by numerically simulations. At the base temperature around $16.5$ mK, the anticrossing in the transmission spectrum, which is so called vacuum Rabi splitting in Fig.~\ref{Figure5}(a), demonstrates that the qubit is strongly coupled to SAW resonator. Thus two peaks appear in the curve A in Fig.~\ref{Figure5}(g). With the increase of the temperature, high energy levels are excited, the anticrossing becomes not very clear in Fig.~\ref{Figure5}(c), two more peaks occurs as shown in curve C$'$ of Fig.~\ref{Figure5}(g) or Fig.~\ref{Figure5} (h). When the temperature is increased to $T_{c}\sim\hbar\omega_{r}\approx149$ mK, around this temperature, anticrossing for the vacuum Rabi splitting gradually disappears as shown in Fig.~\ref{Figure6}(a) and the dominant peaks become wider. Above this temperature, all peaks in the transmission spectra start to merge into each other with the increase of the temperature as shown in Fig.~\ref{Figure6}(i). When the temperature becomes higher and higher, e.g., at the temperature $T\approx349$ mK which means $T\gg T_{c}$, two hybridized modes coalesce as shown in Fig.~\ref{Figure6}(g), thus there is only one peak at the frequency $\omega_r=2\pi\times3.162$ GHz as shown in curve corresponding to $T=349$ mK of Fig.~\ref{Figure6}(i).

\section{CONCLUSIONS}\label{conclusion}
In this work, we experimentally realize a cQAD system, which consists of a transmon qubit strongly coupled to both a 2D SAW resonator and a $1$D microwave transmission line. With the architecture, we demonstrate AC-Stark shift in weak dispersive system with the average phonon number inside SAW resonator in a wide range, which is from $\langle n\rangle\sim0.1$ to $\langle n\rangle\sim1000$. We show the propagating microwave field in the microwave transmission line can be controlled by phonons inside the SAW resonator with the average phonon number down to $\langle n\rangle\sim21$. This is the first step towards control of microwave photons by phonons in a quantum hybrid system on-chip. We also systematically study the quantum-to-classical crossover of the cQAD system. With the increase of the effective temperature of environment, we observe the crossover from level anticrossing to level attraction in transmission spectra. We expect that we could observe clear evidence of a quantum nature of the coupled system by two more peaks in resonant spectra when the environmental temperature is lower than the crossover temperature with improved architecture, e.g., a flux qubit is strongly coupled to a SAW resonator and a microwave transmission line simultaneously in future. Our study lays a solid foundation for further experiments in the context of quantum physics and information using cQAD system.

\begin{acknowledgments}
This work is supported by NSFC under Grant No.~61833010, No.~12074117 and No.~12061131011, Hunan Province Science and Technology Innovation Platform and Talent Plan (Excellent Talent Award) under Grant No.~2017XK2021. LZ is supported by the NSFC under Grant Nos.~11975095 and 11935006. LMK is supported by the NSFC under Grant No.~11775075 and No.~11434011. YXL is supported by the Key R\&D Program of Guangdong province under Grant No.~2018B030326001 and NSFC under Grant No.~11874037.
\end{acknowledgments}


\section{\uppercase{Introduction}}
	The fast evaluation of kernel density estimators has been well-studied including approaches based on the fast Gauss transform \cite{GreStr91}, hierarchical space decompositions \cite{GreRok87}, locality sensitive hashing \cite{ChaSim17,BacChaInd18,SimRonBai19,BacIndWag19}, and binning \cite{Sco85}, as well as interpolation \cite{Jon98,Kog98}, our main technique in this work. Typically these techniques carefully leverage the structure of the kernel under consideration, and many of them operate in a worst-case framework over the dataset. In this work, we consider the problem of fast evaluation of a density estimator $\hat f$ in a statistical setting where $\hat f$ gives a good pointwise approximation to an unknown density $f$  that lies in a H\"{o}lder class of smooth functions. We show that a pointwise approximation guarantee alone, without assuming any specific structure of the estimator $\hat f$, is enough to construct a new estimator $\tilde f$ that can be stored and queried cheaply, and whose approximation error is similar to that of the original estimator. Our approach is based on a multivariate polynomial interpolation scheme of \citet{Nic72} \citep[see also][]{ChuYao77} and provides an explicit formula for $\tilde f$ in terms of some judiciously chosen queries of the original estimator.

	\subsection{Background and related work}

    Density estimation is the task of estimating an unknown density $f$ given an i.i.d. sample $X_1, \ldots, X_n \sim \mb{P}_f$, where $\mb{P}_f$ is the probability distribution associated to $f$. A popular choice of density estimator is the kernel density estimator (KDE)

    \begin{equation}
        \label{EQ:KDE}
        \hat f(y):=\frac{1}{nh^d}\sum_{j=1}^n K\left(\frac{X_i-y}{h}\right).
    \end{equation}

    With proper setting of the bandwidth parameter $h$ and choice of kernel $K$, the KDE $\hat f$ is a minimax optimal estimator over the $L$-H\"{o}lder smooth densities $\mc{P}_{\mc{H}}(\beta, L)$ of order $\beta$ \citep[see e.g.][Theorem~1.2]{Tsy09}:

    \begin{equation}
            \label{eqn:classical_risk}
             \inf_{ \hat{f} } \sup_{ f \in \mc{P}_{\mc{H}}(\beta, L)  }\mb{E}_f \, \norm{ \hat{f} - f }_2    = \Theta_{\beta, d, L}( n^{-\frac{\beta}{2 \beta + d} } )\,.
    \end{equation}

    Despite its statistical utility, the KDE \eqref{EQ:KDE} has the computational drawback that it naively requires $\Omega(n)$ time to evaluate a query. The problem of improving on these computational aspects has thus received a lot of attention.

    Motivated by multi-body problems, \cite{GreStr91} developed the fast Gauss transform to rapidly evaluate sums of the form \eqref{EQ:KDE} when $K(x) = \exp(-\abs{x}_2^2)$ is the Gaussian kernel. Their work is posed a worst-case batch setting where $\hat f$ is to be evaluated at $m$ points $y_1, \ldots, y_m$ specified in advance and the locations $X_1, \ldots, X_n$ lie in a box. Their techniques use hierarchical space decompositions and series expansions to show that \eqref{EQ:KDE} may be evaluated at $y_1, \ldots, y_m$ with precision $\eps$ in time $h^d (\log \frac{1}{\eps})^d (n + m)$. These results apply to any kernel that has a rapidly converging Hermite expansion \citep[see also][]{GreRok87}. There are also follow up works on the improved fast Gauss transform and tree-based methods that use related ideas \citep{YanDurGum03,LeeMooGra06}.

    More recently, several works \citep{ChaSim17,BacChaInd18,SimRonBai19,BacIndWag19} are devoted to the problem of fast evaluation of \eqref{EQ:KDE} in high dimension using locality sensitive hashing. In these works, the dataset is carefully reweighted for importance sampling such that a randomly drawn datapoint $X_r$'s corresponding kernel value $K(X_r - y)$ gives a good approximation to $\hat f (y)$. This sampling procedure can be executed efficiently using hashing-based methods. For example, \citet{BacIndWag19} show that for the Laplace and Exponential kernels with bandwidth $h = 1$, e.g., the value $\hat f(y)$ can be computed with multiplicative $1 \pm \eps$ error in time $O(\frac{d}{\sqrt{\tau}\eps^2})$ even in worst case over the dataset, where $\tau$ is a uniform lower bound on the KDE.

    Another effective approach to this problem in high dimensions is through coresets \citep{AgaHarVar05,Cla10,PhiTai18,PhiTai19}. A coreset is a representative subset $\{X_i\}_{i \in S}$ of a dataset such that
    \[
    \hat f(y) \approx \frac{1}{nh^d} \sum_{i \in S} K\left(\frac{X_i-y}{h}\right).
    \]
    When $h = O(1)$, for example, the results of \citet{PhiTai19} give a polynomial time algorithm in $n, d$ such that the coreset KDE yields an additive $\eps$ approximation to $\hat f$ using a coreset of size $\tilde{O}(\frac{\sqrt{d}}{\eps})$. Their results hold in worst case over the dataset and for a variety of popular kernels. The methods of \citet{PhiTai19} are powered by state-of-the-art algorithms from discrepancy theory \citep{BanDadGarLov18} \cite[see][for a comprehensive exposition on discrepancy]{Mat99,Cha00}.

    Our approach is most closely related to prior work on the interpolation of kernel density estimators due to \citet{Jon98} and \citet{Kog98}. Motivated by visualization and computational aspects, \citet{Jon98} studies binned and piecewise linearly interpolated univariate kernel density estimators and provides precise bounds on the mean-integrated squared error. \citet{Kog98} extends this work and constructs higher order piecewise polynomial interpolants of multivariate kernel density estimators, and shows that for very smooth densities, this procedure improves the mean-integrated squared error. We also note the recent work of \cite{BelRakTsy19,LiaRak20} demonstrating the perhaps surprising effectiveness of interpolation in nonparametric regression.

    Our work differs from \cite{Kog98} in a few important respects. We do not assume $\hat f$ to be a KDE in the first place, but rather give a general method for effectively interpolating a minimax density estimator. Also, our results hold for the entire range of smoothness parameters $\beta$ and dimension $d$, while \cite{Kog98} requires the density to be at least $qd$ times differentiable when interpolating KDEs with kernels of order $q$ \citep[Definition 1.3]{Tsy09}. On the other hand, our method increases the mean squared by a multiplicative factor $\tilde{O}(c_{\beta, d}),$ while Kogure's approach improves the mean squared error (though our focus here is the sup norm). Finally, we use a different interpolation scheme as detailed in Section \ref{sec:principal_lattice}.

    \subsection{Results}
    We seek to impose minimal requirements on a density estimator $\hat f$ of an unknown density $f$ so that it can be converted to a new estimator $\tilde{f}$ that performs well on the following criteria.

    \begin{enumerate}[itemsep=3pt,parsep=3pt]
    \item (\textbf{Minimax}) $\tilde{f}$ is a minimax estimator for $f$
    \item  (\textbf{Space-efficient}) $\tilde{f}$ can be stored efficiently
    \item (\textbf{Fast querying}) $\tilde{f}$ can be evaluated efficiently
    \item (\textbf{Fast preprocessing}) $\tilde{f}$ can be constructed efficiently
    \end{enumerate}
    In the statistical setup where typically $\beta, d = O(1)$, by \textit{efficient} we mean requiring only polynomial time or space in the sample size $n$. In particular for fixed $\beta$, by \eqref{eqn:classical_risk} consistent estimation is only possible when $d \ll \log n$. In what follows we indicate dependencies on the parameters $\beta$ and $d$.

    The requirement that we place on the estimator $\hat f$ to be converted is the following assumption.

    \begin{assumption}
         \label{ass:tail_bound}
         For all $y \in [0, 1]^d$ and $1 \geq t \geq  \eps$, we have
             \begin{equation*}
             	\label{eqn:pointwise_tail_bound}
             	\sup_{f \in \hd(\beta, L)} \mb{P}_f\left[ \abs{ \hat f(y) - f(y)} > t \right] \leq
             	2 \exp\left( - \frac{t^2}{\eps^2} \right),
             \end{equation*}
         where $\eps := c^* \, n^{-\beta/(2\beta + d)}$ is the minimax rate over all density estimators and $c^* = c_{\beta, d, L} > 0$.
    \end{assumption}

    In particular this is satisfied if the pointwise error is a sub-Gaussian random variable with parameter $\eps$ that captures the minimax rate of estimation. For the KDE built from a kernel $K$ of order $\ell:= \flo{\beta}$ \citep[][Definition 1.3]{Tsy09} and bandwidth $h = n^{-1/(2\beta + d)}$, \eqref{eqn:pointwise_tail_bound} follows from a standard bias-variance trade-off and an application of Bernstein's inequality for bounded random variables \citep[see e.g.][]{Ver18}. Under Assumption \ref{ass:tail_bound} we have our main result.

\begin{theorem}
\label{thm:main}
    Let $\hat f$ denote a density estimator satisfying Assumption \ref{ass:tail_bound} for some $\beta > 0$ and $d \geq 1$. Let $Q$ denote the amount of time it takes to query $\hat f$. Set $\ell = \flo{\beta}$. Then there exists an estimator $\ti f$ that can be constructed in time $O(Q \binom{\ell +d}{d} n^{\frac{d}{2\beta +d}})$, that requires $O(d^2 \ell^2 \binom{\ell +d}{d} n^{\frac{d}{2\beta +d}} \log n)$ bits to store, that can be queried in time $O(d \ell \binom{\ell + d}{d}\log n)$, and that satisfies
    \begin{equation*}
    \mb{E}_f \norm{\ti f - f}_\infty < \ti c (\log n)^{1/2} n^{-\frac{\beta}{2\beta +d}} ,
    \end{equation*}
    where we may take
    \[
    \tilde{c} = 10 \binom{\ell + d}{\ell}^{3} (2 \ell)^{2 \ell} \left\{ c^* + \frac{Ld^{\beta}}{\ell!} \right\} + 2Ld^{\frac12 (3\ell + 1)}\,.
    \]
\end{theorem}

    In particular, for $\beta, d = O(1)$, we can evaluate queries to $\ti f$ in nearly constant time, and the estimator $\ti f$ can be stored using sublinear space. Moreover, $\ti f$ can be preprocessed in subquadratic time, assuming that the evaluation time of the original estimator $\hat f$ is $O_d(n)$, which holds for the KDE \eqref{EQ:KDE}. Finally, $\ti f$ is a minimax estimator in the sup norm up to logarithmic factors and thus minimax up to logarithmic factors in the $L_p$ norms. In addition, our construction of $\ti f$ in Section \ref{sec:principal_lattice} yields an explicit formula for $\ti f$ in terms of a sublinear number of queries of $\hat f$ at a set of judiciously chosen query points. Specifically, the estimator $\ti f$ is a piecewise multivariate interpolation of the estimator $\hat f$ on this lattice of query points.

	\subsection{Setup and notation}

    Fix an integer $d \ge 1$. For any  multi-index $s = (s_1,\ldots, s_d) \in \mathbb{Z}_{\geq 0}^d$, let $\abs{s} = s_1 + \cdots + s_d$ and for $x =(x_1, \ldots, x_d) \in \mathbb{R}^d$, define $s!=s_1!\cdots s_d!$ and $x^s=x_1^{s_1}\cdots x_d^{s_d}$. Let $D^s$ denote the differential operator
    $$
    D^s=\frac{\partial^{|s|}}{\partial x_1^{s_1} \cdots \partial x_d^{s_d}}\,.
    $$
    Fix a positive real number $\beta,$ and let $\lfloor \beta \rfloor$ denote the maximal integer \textit{strictly} less than $\beta$. We reserve the notation $\norm{ \cdot }_p$ for the $L_p$ norm

    and $\abs{\cdot}_p$ for the $\ell_p$ norm.

    Given $L>0$ we let  $\mc{H}(\beta, L)$ denote the space of H\"{o}lder functions $f:\mb{R}^d \to \mb{R}$ that are supported on the cube $[0,1]^d$, are $\lfloor \beta \rfloor$ times differentiable, and satisfy
    \begin{equation*}
    |D^s f(x) - D^s f(y)| \leq L \abs{x - y}_2^{\beta - \lfloor \beta \rfloor}\,,\quad
    \end{equation*}
    for all $x, y \in \mb{R}^d$ and for all multi-indices $s$ such that $\abs{s} = \flo{\beta}$.

    Let $\mc{P}_{\mc{H}}(\beta, L)$ denote the set of probability density functions contained in $\mc{H}(\beta,L)$. For $f \in \mc{P}_{\mc{H}}(\beta,L)$, let $\mb{P}_f$ (resp. $\mb{E}_f$) denote the probability distribution (resp. expectation) associated to $f$.

    The parameter $L$ will be fixed in what follows, so typically we write $\mc{P}_{\mc{H}}(\beta) := \mc{P}_{\mc{H}}(\beta,L)$. The constants $c, c_{\beta, d}, c_{L},$ etc. vary from line to line and their subscripts indicate parameter dependences.

\section{\uppercase{Efficient interpolation of density estimators}}
\label{sec:interpolation}

The important implication of Assumption \ref{ass:tail_bound} is that we can query $\hat f$ at a polynomial number of data points such that for each query $y$, $\hat f(y) \approx f(y)$, where $f$ is the unknown density.

\begin{lemma}
\label{lem:union_bound} Let $A > 0$ and set $N = \Delta n^A$ with $\Delta >1$. Let $y_1, \ldots, y_N \subset [0,1]^d$ denote a set of points. Then with probability at least $1 - 2n^{-2}$,
\[
\abs{ \hat f(y_i) - f(y_i) } \leq \sqrt{ \log( 2 \Delta n^{A+2} ) } \, \, \eps
\]
for all $1 \leq i \leq N$, where $\eps = c_{\beta, d, L} \, n^{-\beta/(2\beta + d)}$ is the minimax rate.
\end{lemma}

\begin{proof}
Set $t = \sqrt{ \log( 2 \Delta n^{A+2} ) } \,  \eps \geq  \eps$ and apply Assumption \ref{ass:tail_bound} to $y_i$. Then by the union bound,
\[
\mb{P}\left[ \exists y_i\,:\, \abs{ \hat f(y_i) - f(y_i) } > t \right] \leq 2 \Delta n^A  e^{-\frac{t^2}{\eps^2}} \leq n^{-2}.
\]
\end{proof}

We now describe our construction of $\ti f$. Define $\ell := {\lfloor \beta \rfloor}$ and $M = \binom{\ell + d}{\ell}$.

\medskip

\noindent \textbf{Construction of} $\tilde{f}$ (\textit{informal}):

\begin{enumerate}[itemsep=3pt,parsep=3pt]
	\item {\sc Partition}: Divide $[0,1]^d$ into $h^{-d}$ sub-cubes $\{I_{\vec{j}}\} \subset [0,1]^d$ of side-length $h = n^{-1/(2\beta + d)}$ where $\vec{j} \in \mb{Z}_{\geq 0}^d$ and $I_{\vec{j}} := [0, h]^d + h \vec{j}$.
	\item {\sc Mesh}: For each $\vec{j}$, construct a mesh consisting of $M = \binom{\ell + d}{\ell}$ points $U^{\vec{j}}_1, \ldots, U^{\vec{j}}_M \in I_{\vec{j}}$.
	\item {\sc Interpolate}: In each sub-cube $I_{\vec{j}}$, construct a multivariate polynomial interpolant $\hat q_{\vec{j}}$ on the $M$ points $( U_{1}\V{j} , \hat f(U^{\vec{j}}_1  )), \ldots, ( U_{M}\V{j} , \hat f( U^{\vec{j}}_M) )$.
\end{enumerate}

\textbf{Return}: $\tilde{f}: [0,1]^d \to \mb{R}$ defined by
\[
\tilde{f}(y) = \sum_{\vec{j}} \hat q_{\vec{j}}(y) {\rm 1}\kern-0.24em{\rm I}(y \in I_{\vec{j}}).
\]

We first give some intuition for why $\tilde{f}$ is an accurate estimator. On each sub-cube $I_{\vec{j}}$, the true density $f \in \mc{P}_{\mc{H}}(\beta,L)$ is approximated up to the minimax error by a polynomial $q_{\vec{j}}$ of degree at most $\ell$ by the properties of H\"{o}lder functions. Upon setting $\Delta = M$ and $A = d/(2\beta+d)$ in Lemma \ref{lem:union_bound}, this guarantees that for all queries $U\V{j}_k$ in the mesh, $\hat f(U\V{j}_k) \approx f(U\V{j}_k) \approx q_{\vec{j}}(U\V{j}_k)$ with high probability. By studying the stability of the resulting polynomial system of equations, we can show that this construction yields a good approximation to the `true' interpolation polynomial $q_{\vec{j}}$ on the sub-cube $I_{\vec{j}}$. This argument, carried out formally in the remainder of this section, yields the estimation bound of Theorem \ref{thm:main}.

Next, we comment on the remaining guarantees of Theorem \ref{thm:main}. As we show later, there is an explicit formula for $\hat q_{\vec{j}}$ so the main preprocessing bottleneck is evaluation $\hat f$ on the $M n^{d/(2\beta + d)}$ queries in the mesh, which naively takes $Q M n^{d/(2\beta + d)}$ time. For the space requirement, it suffices to store the queried values $\{ \hat f(  U\V{j}_k) \}$ up to polynomial precision as well as the elements of the mesh. Hence $O(d^2 \ell^2 \binom{\ell +d}{d} n^{\frac{d}{2\beta +d}} \log n)$ bits suffice, by the uniform boundedness of H\"{o}lder functions (see Lemma \ref{lem:bounded}) and Assumption \ref{ass:tail_bound}. Finally, to query $\hat f$ at a point $y \in [0,1]^d$ requires checking which sub-cube $y$ belongs to by scanning its $d$ coordinates and then evaluating $\hat q_{\vec{j}}(y)$, which is a $d$-variate polynomial of degree $\flo{\beta}$. These considerations lead to the guarantees of Theorem \ref{thm:main}.

\subsection{Interpolation on the principal lattice}
  \label{sec:principal_lattice}

    To construct our interpolant, we refer to the following definition and theorem which are classical in finite element analysis \citep{Nic72,ChuYao77}. The lattice $\mc{P}_\ell$, dubbed the $\ell$-th principle lattice, has the special property that every function defined on $\mc{P}_\ell$ admits a unique polynomial interpolant of degree at most $\ell$. This translates to the invertibility of the associated Vandermonde-type matrix and, equivalently, means that $\mc{P}_\ell$ is not a subset of any algebraic hypersurface of degree at most $\ell$.

    \begin{definition}[$\ell$-th principal lattice of $\Delta_d$]
    	\label{def:principal_lattice}
	Let $\Delta_d \subset [0,1]^d$ denote the simplex on the points $\{0\} \cup \{ e_i \}_{i = 1}^d \subset \mb{R}^d$, where $e_i$ denotes the $i$-th standard basis vector in $\mb{R}^d$. Label the vertices of $\Delta_d$ to be $v_0 = 0, v_i = e_i$ for $0 \leq i \leq d$. For all $x \in \mb{R}^d$, there exists a unique vector $(\lambda_0(x),\ldots, \lambda_d(x))$ with entries summing to one and such that
	\[
	x = \sum_{i = 0}^d \lambda_i(x) v_i\,.
	\]
	Let $\Lambda:\mb{R}^d \to \mb{R}^{d+1}$ denote the function such that $\Lambda(x)=(\lambda_0(x),\ldots,  \lambda_d(x))$.
	The \emph{$\ell$-th principal lattice} $\mc{P}_\ell$ of $\Delta_d$ is defined to be
	\begin{equation}
	\label{eqn:principal_lattice}
	\mc{P}_\ell = \left\{ x \in \Delta_d: \, \ell\Lambda(x) \in \mb{Z}^{d+1}_{\geq 0} \right\}\,.
	\end{equation}
\end{definition}

The following classical result gives an explicit form for the interpolant.

	\begin{theorem}[\cite{Nic72},\cite{ChuYao77}]
		\label{thm:classical_interpolate}
Write $\mathcal{P}_\ell =\{U_1, \ldots, U_M\} \subset \Delta_d$ and let $g: \mc{P}_\ell \to \mb{R}$ denote a function defined on this lattice.
Define the polynomial
		\begin{equation}
		\label{eqn:basis_fn}
			p_i(x) = \prod_{\substack{t = 0 \\ \lambda_t(U_i) > 0 }}^{d} \prod_{r = 0}^{\ell  \lambda_t(U_i) - 1} \frac{ \lambda_t(x) - \frac{r}{\ell}}{ \lambda_t(U_i) - \frac{r}{\ell} }\,,
		\end{equation}
where we recall that $\lambda_t(x)$ is from Definition~\ref{def:principal_lattice}. Then
		\[
			p(x) := \sum_{i = 1}^M p_i(x) g(U_i)
		\]
		satisfies $p(U_i) = g(U_i)$ for all $U_i \in \mc{P}_\ell$. Moreover, this is the unique polynomial of degree at most $\ell$ with this property.
	\end{theorem}

	Since $\lambda_t(x)$ is linear in $x \in \mb{R}^d$, it is easy to see that $p_i(x)$ is a polynomial of degree $\ell$, and moreover $p_i(U_j)=1$ if $i=j$ and zero otherwise.

	We are now ready to give a precise description of the construction of $\tilde{f}$. The idea is to generate the mesh for interpolation using a shifted and rescaled version of the $\ell$-th principal lattice on $\Delta_d \subset [0,1]^d$. Recall that $\hat f$ is a density estimator that satisfies \eqref{eqn:pointwise_tail_bound}.

	\medskip

\noindent \textbf{Construction of} $\tilde{f}$ (\textit{formal version}):

\begin{enumerate}[itemsep=3pt,parsep=3pt]
	\item \textsc{Partition}: Divide $[0,1]^d$ into $h^{-d}$ sub-cubes $\{I_{\vec{j}}\} \subset [0,1]^d$ of side-length $h = n^{-1/(2\beta + d)}$ where $\vec{j} \in \mb{Z}_{\geq 0}^d$ and $I_{\vec{j}} := [0, h]^d + h \vec{j}$.
	\item \textsc{Mesh}: For each $\vec{j}$, construct a mesh on $I_{\vec{j}}$ consisting of $M = \binom{\ell + d}{\ell}$ points given by the shifted and rescaled principal lattice $\mc{P}_\ell^{\vec{j}} := \{ h(x + \vec{j}) :\, x \in \mc{P}_\ell  \} \subset I_{\vec{j}}$. Let $U^{\vec{j}}_1, \ldots, U^{\vec{j}}_M$ denote the points in $\mc{P}_\ell^{\vec{j}}$.

	\item \textsc{Interpolate}: In each sub-cube $I_{\vec{j}}$, construct a multivariate polynomial interpolant $\hat{q}_{\vec{j}}$ through the $M$ points $( U_{1}\V{j} , \hat f(U_{1}\V{j}), \ldots, ( U_{M}\V{j} , \hat f(U_{M}\V{j}) )$ given by $\hat{q}_{\vec{j}}(y) = p(y/h - \vec{j})$, where $p$ is the polynomial interpolant from Theorem \ref{thm:classical_interpolate} given by
	\[
		p(x) = \sum_{k = 1}^M p_k(x) \hat f( U^{\vec{j}}_k ) .
	\]

\end{enumerate}

\textbf{Return}: $\tilde{f}: [0,1]^d \to \mb{R}$ defined by
\[
\tilde{f}(y) = \sum_{\vec{j}} \hat{q}_{\vec{j}}(y) {\rm 1}\kern-0.24em{\rm I}(y \in I_{\vec{j}}).
\]

The interpolant constructed in Step (4) is unique by Theorem \ref{thm:classical_interpolate}. Also recall that in each sub-cube $I_{\vec{j}}$, the true density $f \in \mc{P}_{\mc{H}}(\beta)$ is well-approximated by a polynomial $q_{\vec{j}}$ of degree $\ell$. To prove  Theorem~\ref{thm:main} it remains to show that given the values of this polynomial (up to some small error incurred by the estimation error of $\hat f$) on the points in our mesh, the interpolant $\hat{q}_{\vec{j}}$ from Step (4) gives a good approximation to $q_{\vec{j}}$ on the sub-cube $I_{\vec{j}}$.

\subsection{Proof of Theorem~\ref{thm:main}}
\label{sec:proof_interpolation}
First, we quantify the error in the approximation of the values of $q_{\vec{j}}$ on the mesh points. Let $f_{z, \ell}$ denote the degree $\ell$ polynomial given by the Taylor expansion of $f \in \mc{P}_{\mc{H}}(\beta)$ at $z$. Since $f \in \mc{P}_{\mc{H}}(\beta)$, by a standard fact (see Lemma \ref{lem:taylor_holder}) it holds that
\[
|f(y) - f_{z, \ell}(y)| \leq \frac{L d^{\ell/2} }{\ell !} \abs{y-z}_2^\beta,
\]
where $f_{z, \ell}$ is the degree-$\ell$ Taylor expansion of the function $f$ at $z \in \mb{R}^d$.

For $\vec{j} \in \{0, \ldots, h^{-1} - 1\}^d$, define $q_{\vec{j}} := f_{z_{\vec{j}}, \ell}$, where $z_{\vec{j}}$ is the vertex of $I_{\vec{j}}$ closest to the origin. Then for all $y \in I_{\vec{j}}$, it holds that
\begin{align}
\label{eqn:Taylor_approx}
|f(y) - q_{\vec{j}}(y)| &\leq \left( \frac{L d^{\beta} }{\ell !} \right) \, h^{\beta} \nonumber \\
&= \left(\frac{L d^{\beta} }{\ell !}\right) n^{-\beta/(2\beta + d)} \nonumber \\
&=: \hat c n^{-\beta/(2\beta + d)}
\end{align}
Note that the right-hand side is the minimax rate of estimation in \eqref{eqn:classical_risk} up to constant factors.

Next, by Lemma \ref{lem:union_bound} and \eqref{eqn:Taylor_approx} it holds with probability at least $1 -  n^{-2}$ that
\begin{align}
\label{eqn:noisy_interpolate}
\abs{ q_{\vec{j}}(U^{\vec{j}}_k) - \hat f(U^{\vec{j}}_k) }
&\leq (c^*\log 3M + \hat c) (\log n)^{\frac{1}{2}} n^{-\frac{\beta}{2\beta + d}} \nonumber \\
&=: \breve{c} (\log n)^{\frac{1}{2}} n^{-\frac{\beta}{2\beta + d}}
\end{align}
for all $\vec{j} \in \{0, \ldots, h^{-1} -1\}^d$ and $ k \in [M]$. Using this fact, we can show that the polynomial interpolant built on $\{ (U_k^{\vec{j}}, \hat f( U_k^{\vec{j}} ))\}_{k = 1}^{M}$ provides a good approximation for $q_{\vec{j}}$ on the interval $I_{\vec{j}}$, which is our next task. The following lemma establishes stability of the polynomial approximation.

\begin{lemma}
	\label{lem:interpolation_error}
	Let $\hat{q}_{\vec{j}}$ denote the unique polynomial of degree at most $\ell$ that passes through the points $\{( U_k^{\vec{j}}, \hat f( U_k^{\vec{j}} )) \}_{k = 1}^{M}$. Then with probability at least $1 - 2 n^{-2}$, for all $\vec{j}$ and all $x \in I_{\vec{j}}$,
	\begin{equation}
	\label{eqn:interpolation_error}
	\abs{ q_{\vec{j}}(x) - \hat{q}_{\vec{j}}(x) } \leq c_{\beta, d, L} (\log n)^{\frac{1}{2}} n^{-\frac{\beta}{2\beta + d}}\,.
	\end{equation}
\end{lemma}
\begin{proof}
	Define $g\J(x) = q\J( h(x + \vec{j}) )$ to be a polynomial function $g:[0,1]^d \to \mb{R}$, noting that $h(x + \vec{j}) \in I\J$. Write
	\[g\J(x) = \sum_{\alpha: \abs{\alpha} \leq \ell} s_\alpha x^{\alpha}
	\]
	where $s_\alpha \in \mb{R}$ and $\alpha \in \mb{Z}^d_{\geq 0}$ is a multi-index. Let $V_1, \ldots, V_M$ denote the elements of $\mc{P}_\ell$, and let $s$ denote the vector of coefficients of $g\J$, indexed by the multi-index $\alpha$. Define $y \in \mb{R}^M$ by $y_i = g\J(V_i)$. Consider the Vandermonde-type $M\times M$ matrix $\mc{V}$ defined by $\mc{V}_{i, \alpha} := V_i^\alpha$, where the columns of $\mc{V}$ are indexed by the multi-index $\alpha$ with $|\alpha|\leq \ell$. Then it holds that $\mc{V} s = y$.

	By Theorem \ref{thm:classical_interpolate}, the matrix $\mc{V}$ is invertible. Let $\sigma_0>0$ denote its minimum singular value, noting that $\sigma_0$ is a constant that only depends on $\beta$ and $d$.
Let $\hat{y} \in \mb{R}^M$ denote the vector defined by $\hat{y}_k = \hat f( U_k^{\vec{j}} )$, where $k \in [M]$ indexes the points in the principal lattice $\mc{P}_\ell$, and let $\hat{s} = \mc{V}^{-1} \hat{y}$. By \eqref{eqn:noisy_interpolate},
	\[
		\abs{\hat{y} - y}_2 \leq \sqrt{M} \breve{c} (\log n)^{\frac{1}{2}} n^{-\frac{\beta}{2\beta + d}}.
	\]

	Hence,
	\begin{align*}
		\abs{\hat{s} - s}_2^2 &= \abs{\mc{V}^{-1}(\hat{y} - y)}_2^2  \leq \frac{M \breve{c}^2}{\sigma_0} (\log n)^{\frac{1}{2}} n^{-\frac{\beta}{2\beta + d}}
	\end{align*}

	Define $\hat{q}\J(x) = \hat{g}\J( x/h - \vec{j} ),$
	and note that $\hat{q}\J$ passes through the points $\{ (U_k\V{j}, \hat f( U_k\V{j} )) \}_{k = 1}^M$. Then for all $x \in I_{\vec{j}}$, by the previous display and Cauchy--Schwarz, we have
\begin{align*}
\abs{\hat{q}\J(x) - q\J(x)} &\leq \sum_{\alpha: |\alpha| \leq \ell}\abs{(\hat{s}_\alpha - s_\alpha) \left( x/h - \vec{j} \right)^\alpha } \\
&\leq \frac{M \breve{c}}{\sqrt{\sigma_0}} (\log n)^{\frac{1}{2}} n^{-\frac{\beta}{2\beta + d}}\,.
\end{align*}
\end{proof}

	Define $\tilde{f}(x) = \sum\J \hat{q}\J(x) {\rm 1}\kern-0.24em{\rm I}(x \in I_{\vec{j}})$, and observe that  Theorem~\ref{thm:main} follows from \eqref{eqn:Taylor_approx}, Lemma \ref{lem:interpolation_error}, and the triangle inequality. Though we have derived a high probability bound, the expectation claimed in Theorem \ref{thm:main} follows using the uniform boundedness of H\"{o}lder functions. Tracing constants above and applying Lemma \ref{lem:Vandermonde} below yields the expression for $\tilde{c}$.

\subsection{Stability of the Vandermonde-type matrix}

\label{appendix:singular_values}

We give below an explicit lower bound on the smallest singular value of the Vandermonde-type matrix $\mc{V}$ from Lemma \ref{lem:interpolation_error} associated to the $\ell$-th principal lattice on the simplex. This gives an explicit dependence of our results on the smoothness $\beta$ and dimension $d$.

	\begin{lemma}
	\label{lem:Vandermonde}
		Let $\mc{V}$ denote the Vandermonde-type matrix associated to the $\ell$-th principle lattice $\mc{P}_\ell$ of $\Delta_d$. To be explicit, let $x_1, \ldots, x_M \in [0,1]^d$ denote the elements of $\mc{P}_\ell$, where $M = \binom{\ell + d}{\ell}$. Given a multi-index $\alpha$ with $|\alpha| \leq \ell$, define $\mc{V}_{k, \alpha} := x_k^\alpha$. Let $\sigma_0$ denote the minimum singular value of $\mc{V}$. Then
		\[
		\sigma_0 \geq \binom{\ell + d}{\ell}^{-3} 4^{-\ell} \ell^{-2\ell}.
		\]
	\end{lemma}

	\begin{proof}
		Write
		\[
		p_k(x) = \sum_{\alpha: \abs{\alpha} \leq \ell} u_{k, \alpha} x^{\alpha},
		\]
		where $p_k$ was defined in Theorem  \ref{thm:classical_interpolate} for $1 \leq k \leq M$.

		Then write $u_{k, \alpha} = \omega_k \cdot v_{k, \alpha}$, where
		\begin{equation}
		\label{eqn:normalization}
		\omega_k := \prod_{\substack{t = 0 \\ \lambda_t(x_k) > 0}}^d \prod_{r = 0}^{\ell \lambda_t(x_k) - 1} \frac{1}{\left(\lambda_t(x_k) - \frac{r}{\ell}  \right)}\,.
		\end{equation}
		Since $\Lambda(x_k) = \frac{1}{\ell}(s_0(x_k), \ldots, s_d(x_k))$, where $s_j(x_k) \in \mb{Z}_{\geq 0}$ for all $0 \leq j \leq d$, we conclude that the denominator of each term of \eqref{eqn:normalization} is at least $1/\ell$. Thus
		\begin{equation}
		\label{eqn:normalization_bound}
		\omega_k \leq \ell^{\ell},
		\end{equation}
		since there are $\ell$ terms in the product defining $\omega_k$.

		Note that
		\[
		\lambda_t(x_k) = \begin{cases}
		x_{kt} \quad &\text{if} \, \, 1 \leq t \leq d \\
		1 - \sum_{t = 1}^d x_{kt} \quad &\text{if} \, \, t = 0.
		\end{cases}
		\]
		Therefore, for all $1 \leq k \leq M$, we may write
		\[
		\tilde{p}_k(x) := \frac{p_k(x)}{\omega_k} = \prod_{s = 1}^\ell ( m_{s, k} \cdot x + h_{s, k} ),
		\]
		where $m_{s, k} \in \mb{R}^d$ with $\abs{ m_{s, k} }_{\infty} \leq 1$, and $|h_{s, k}| \leq 1$. It follows that the coefficient $v_{k, \alpha}$ of the monomial $x^\alpha$ in $\tilde{p}_k(x) $ satisfies
		\begin{equation}
		\label{eqn:coefficient_bound}
			| v_{k, \alpha} | \leq \binom{\ell}{|\alpha|}.
		\end{equation}

		Next, observe that $\mc{V}^{-1}_{\alpha, k} = u_{k, \alpha}$. Hence, by \eqref{eqn:normalization_bound} and \eqref{eqn:coefficient_bound}, for all multi-indices $\alpha$ with $|\alpha| \leq \ell$,
		\begin{equation}
		\label{eqn:invert_bound}
			\sum_{k = 1}^M |u_{k, \alpha}| \leq \binom{\ell + d}{\ell} \binom{\ell}{|\alpha|} \ell^\ell \leq \binom{\ell + d}{\ell} 2^\ell \ell^\ell.
		\end{equation}

		Letting $\eta \in \mb{R}^M$ denote a unit vector and applying the previous inequality,
		\begin{align*}
		\eta^T \mc{V}^{-T} \mc{V}^{-1} \eta &\leq M \cdot  \max_{\alpha} \left( \sum_{k = 1}^M |u_{k, \alpha}| \right)^2 \\
		&\leq \binom{\ell + d}{\ell}^3 4^\ell \ell^{2\ell}.
		\end{align*}
		The conclusion of Lemma \ref{lem:Vandermonde} follows.
	\end{proof}

    \section{\uppercase{A result of Kolmogorov and Tikhomirov}}

    	Given a function class $\mathcal{F}$, let $N( \mc{F}, \delta)$ denote the minimal number of cubes of side-length $\delta$ that cover $\mathcal{F}$, and define $H(\mc{F}, \delta) = \log N( \mc{F}, \delta)$ to be the metric entropy. A classical result of \citet{Tik93} shows that
            \begin{equation}
            \label{eqn:covering_number}
            H( \hd(\beta), \delta) \leq c_{\beta, d, L} \, \delta^{-\frac{d}{\beta}}.
            \end{equation}

            Their proof strategy is conceptually similar to our piecewise multivariate polynomial approximation scheme in that they subdivide the cube as we do here, approximate $f$ by its Taylor polynomial in each cube, and then discretize the coefficients. Our techniques imply a slightly weaker version of the bound \eqref{eqn:covering_number}.

            Define a mesh as in steps 1 and 2 of our formal construction of $\ti f$ as in Section \ref{sec:principal_lattice}, but now for a general parameter $h > 0$ to be set later. This mesh has $M h^{-d}$ points that we denote by $\{ U\V{j}_k \}_{\vec{j}, k }$. Let $f, g \in \hd(\beta)$ be such that for all $\vec{j}, k$ it holds that
            \[
            \abs{f(U\V{j}_k) - g(U\V{j}_k)} \leq h^\beta.
            \]
            By the H\"{o}lder condition and Lemma \ref{lem:taylor_holder}, there exists a degree $\ell = \flo{\beta}$ polynomial $q_{\vec{j}}$ approximating $f$ in $I_{\vec{j}}$ and a degree $\ell = \flo{\beta}$ polynomial $r_{\vec{j}}$ approximating $g$ in $I_{\vec{j}}$, each with error $h^\beta$ pointwise. We conclude that
            \[
            \abs{ q_{\vec{j}}(U\V{j}_k) - r_{\vec{j}}(U\V{j}_k)} \leq c_{\beta, d, L} \, h^{\beta}
            \]
            for all $\vec{j}, k$. Following the proof of Lemma \ref{lem:interpolation_error} and using the bounds on the singular values of $\mathcal{V}$, this implies that for all $x \in I_{\vec{j}}$,
            \[
            \abs{ q_{\vec{j}}(x) - r_{\vec{j}}(x) } \leq c_{\beta, d, L} \, h^\beta.
            \]
            Hence we conclude that for all $x \in [0,1]^d$,
            \[\abs{f(x) - g(x)} \leq c_{\beta, d, L} \, h^\beta.\]

            The H\"{o}lder densities are uniformly bounded by some constant $c_{\beta, d, L}$ (see Lemma \ref{lem:bounded}). Hence setting $\delta = h^{\beta}$ and rounding the values of each density at each point $U\V{j}_k$ to multiples of $\delta$, we see that there exists an $\delta$-net of size at most
            \[
            \left(  \frac{c_{\beta, d, L}}{ \delta } \right)^{M \delta^{-d/\beta}}.
            \]
            Therefore
            \[
             H( \hd(\beta), \delta) \leq c_{\beta, d, L} \, \delta^{-\frac{d}{\beta}} \log \frac{1}{\delta},
            \]
            a mildly weaker bound than \eqref{eqn:covering_number}.