\chapter{Thermometry}

\noindent This chapter discusses the theoretical thermometry schemes regarding optical lattices.
Throughout, I restrict the discussion to bosons.
The chapter is based on one of my publications \cite{thermometry}.


\section{Introduction}

\subsection{Temperature determination}

Temperature is one of the important physical quantity that determines the phase of any particular system.
In the classical world, for instance, water exists as either solid ice, liquid water or gaseous steam at different temperature and pressure.
In this example, the temperature can be indicated by a mercury-in-glass thermometer, only because its temperature dependence of the volume property has been previously known and calibrated.
In general, a thermometer measures a temperature-dependent property of a subsystem in thermal equilibrium with the rest of the system.
Here, the disturbance to the subsystem is assumed to be negligible in the measurement process, such that it gets quickly relaxed back to equilibrium by the bulk.
\\ \\
\noindent In the quantum limit towards atomic scale, temperature determination poses tough challenges to experimentalists \cite{bloch2}.
To be specific, validity and accuracy are the two concerns.
First, quantum measurements can cause a non-negligible disturbance to the subsystem, because in the atomic limit, the latter is not small after all with respect to the bulk.
We say such measurements to be {\it{destructive}}, i.e. destroying the quantum state of the system, thus questioning fundamentally on the validity of such temperature determination.
Second, the accuracy of quantum measurements is largely challenged by the Heisenberg uncertainty principle in the atomic limit.
No longer could a single measurement be sufficient to determine the temperature accurately.
To enhance on its accuracy, either more information regarding the system must be made available, or more repetitive measurements have to be performed.
The latter, unfortunately, makes temperature determination of quantum systems very {\it{expensive}}. 
\\ \\
In particular, this chapter is only devoted to the subject matter of accurate temperature determination for bosons parabolically-trapped in quantum optical lattices.


\subsection{Physical system}

In the tight-binding limit, parabolically-trapped bosons in a 3D optical lattice can be accurately modelled by the single-band boson Hubbard Hamiltonian \cite{jaksch1}
%
\begin{equation}
\hat{H} - \mu N = -t \sum_{\langle i,j \rangle} \hat{b}_i^\dag \hat{b}_j + \frac{U}{2} \sum_i n_i (n_i -1) - \sum_i ( \mu - V_T r_i^2 ) n_i  \tag{\ref{trapped_boson_hubbard_hamiltonian}}
\end{equation}
%
where $\hat{H}$, in this chapter, specifically refers to the canonical Hamiltonian, and a total number of $N$ bosons exists in the optical lattice with chemical potential $\mu$ and external harmonic trapping $V_T$. 
The nearest-neighbour hopping $t$, and the onsite repulsion strength $U$ are derived from the lattice laser potential $V_0$ through band-structure calculations. 
In the Fock representation $\{ | n_i \rangle \}$, $n_i$ represents the number of bosons at site $i$, and $\hat{b}_i$ ($\hat{b}_j^\dag$) is the annihilation (creation) operator at site $i$ ($j$). 
In the last term, ${\bf{r}}_i$ denotes the displacement vector of site $i$ from the centre of the parabolic trap.
\\ \\
\noindent The thermodynamical properties of the system are theoretically summarised in the partition function $Z = \mathrm{Tr} \exp [ -\beta (\hat{H} - \mu N)]$ with inverse temperature $\beta = 1/k_BT$ \cite{statistical_mechanics}. 
The physics of trapped bosons in an optical lattice (equation \ref{trapped_boson_hubbard_hamiltonian}) is summarized in the phase diagram illustrated in figure \ref{boson_hubbard_model_phase_diagram}, where its phase is also determined by the local density %
\begin{equation}
\label{eq_local_density}
\langle n({\bf{r}}) \rangle = \frac{1}{Z} \, \mathrm{Tr} \, n({\bf{r}}) \exp [ -\beta (\hat{H} - \mu N)] \,\, .
\end{equation}
%
In the limit towards lower densities, the system becomes a normal fluid.
%
\begin{figure}[h]
\centering
\includegraphics[width=14cm]{./figures/density-1.jpg}
\caption{Cross sectional density $n(r)$ of a 3D bosonic $^{87}$Rb optical lattice system with $N=125 000$, averaged over 1000 independent measurements obtained from a QMC simulation. Here, the parameters are $U/t = 10$ (left), $50$ (right), $T/t=1$, and $V_T/t = 0.0277$.}
\label{thermometry_figure_1}
\end{figure}
%
In the presence of a confining harmonic trap, the densities become inhomogeneous.
In particular, the densities become sparse in the wings due to huge potential energy barrier.
Therefore, the wings of the bosonic trapped system exhibit always the normal-fluid phase.
\\ \\
Recently, direct comparisons over time-of-flight images \cite{trotsky}, and radial density profiles \cite{shiang1}, of bosonic optical lattices between theory and experiment have been studied.
Excellent quantitative agreements have been observed, and this builds the confidence over the use of ab-initio Quantum Monte Carlo (QMC) (see chapter \ref{chapter_QMC}) simulations of the single-band boson Hubbard Hamiltonian as a mimic of the bosonic optical lattice experiments.  
In this chapter, independent experimental measurements are mimicked by uncorrelated QMC measurements, where the autocorrelation time of $\tau < 0.2$ has been strictly imposed for the uncorrelation criterion.
Figure \ref{thermometry_figure_1} is an example.


\subsection{Single site addressability}

Recent progress over the use of florescence techniques in bosonic optical lattice experiments has enabled single-site resolution/addressability of in-situ density images \cite{bloch1}.
%
\begin{figure}[h]
\centering
\includegraphics[width=14cm]{./figures/single-site-addressability.jpg}
\caption{A single measurement of atom distribution for an ultracold quantum gas held in a two-dimensional optical lattice. 
The bosons, indicated by bright spots, are confined towards the parabolic trap centre.
The interaction strength U/t increases from left to right, and thus transiting the system from being a superfluid into being a Mott insulator.
Recent experimental advancement enables good visualisation of the Mott plateau (middle) around the trap centre up to single-site resolution.
At extremely large interaction strength $U/t$, there exists a high probability to locate 2 bosons per site in the vicinity around the trap centre. 
However, current fluorescence experiments could not distinguish a doublon from a hole, therefore indicating dark spots around the centre (right). 
{\it{This figure is replicated from Nature 467, 68 (2010).}} \cite{bloch1}
}
\label{single_site}
\end{figure}
%
Figure \ref{single_site} illustrates a single measurement of atom distribution for a 2D bosonic optical lattice at different interaction strengths $U/t$ \cite{bloch1}.
The single-site resolution gives the possibility to extract more information from the system.
Having more inputs such as density-density correlations and incompressibility enhance the accuracy of temperature determination. 
\\ \\
\noindent For 3D optical lattices, only the integrated densities
%
\begin{equation}
\langle n({\bf{\rho}}) \rangle = \int  \langle n({\bf{r}}) \rangle \, dz
\end{equation}
%
along the line-of-sight, i.e. z-direction, can be captured by density imaging.


\section{Fluctuation-dissipation thermometry}

\subsection{Theory}

\begin{figure}
\centering
\includegraphics[width=10cm]{./figures/FD-1.jpg}
\caption{Illustration of the quantities entering the FD thermometry formula. Shown from top to bottom are: (1) cross-sectional density $n(r)$, (2) column-integrated density $n(\rho)$, (3) dissipation term $L(\rho)$, (4) fluctuation term ($\xi = 3$) $R_3(\rho)$, and (5) fluctuation term ($\xi = \infty$) $R_\infty(\rho)$.
We take a 3D bosonic $^{87}$Rb optical lattice system with $N=125 000$, and we average over 1000 independent measurements obtained from a QMC simulation. The parameters in the left column are $U/t=10$, $T/t=1$, and the parameters in the right column are $U/t=50$, $T/t=1$. The trapping frequency is $V_T/t=0.0091$ (left), $0.0277$ (right).
}
\label{FD_1}
\end{figure}

The theory behind fluctuation-dissipation thermometry starts with the differentiation of  $\langle n({\bf{\rho}}) \rangle$ w.r.t. $\mu$, i.e.
%
\begin{eqnarray}
\frac{\partial \langle n({\boldsymbol{\rho}}) \rangle}{\partial \mu} &=&  \frac{\partial}{\partial \mu} \int dz \, \frac{1}{Z} \mathrm{Tr} \, n({\bf{r}}) \exp [ -\beta (\hat{H} - \mu N)]    \nonumber \\
 &=& \int dz \, \frac{1}{Z} \mathrm{Tr} \, \beta N \, n({\bf{r}})  \exp [ -\beta (\hat{H} - \mu N)]     \nonumber \\
 &  & - \int dz \, \frac{1}{Z^2} \left( \mathrm{Tr} \, \beta N   \exp [ -\beta (\hat{H} - \mu N)]  \right)  \left( \mathrm{Tr} \, n({\bf{r}})   \exp [ -\beta (\hat{H} - \mu N)]  \right)  \nonumber \\
 &=& \beta \left(  \langle n({\boldsymbol{\rho}}) N \rangle  - \langle n({\boldsymbol{\rho}}) \rangle \langle N \rangle    \right)
 \end{eqnarray}
%
or
%
\begin{equation}
k_B T \frac{\partial \langle n({\boldsymbol{\rho}}) \rangle}{\partial \mu} =  \langle n({\boldsymbol{\rho}}) N \rangle  - \langle n({\boldsymbol{\rho}}) \rangle \langle N \rangle \,\, .
\end{equation}
%
In the validity of the local density approximation (LDA), the column-integrated density becomes
%
\begin{equation}
\label{LDA_density}
\langle n({\boldsymbol{\rho}}) \rangle \stackrel{\mathrm{LDA}}{=} \langle n(\mu ; {\boldsymbol{\rho}}) \rangle
\end{equation}
%
where $\mu$ here refers to the local chemical potential in equation \ref{LDA_density}. 
The dissipation term is therefore
%
\begin{equation}
L({\boldsymbol{\rho}}) =  \frac{\partial \langle n({\boldsymbol{\rho}}) \rangle}{\partial \mu} \stackrel{\mathrm{LDA}}{=} - \frac{1}{2V_T} \frac{1}{ \rho} \frac{\partial \langle n({\boldsymbol{\rho}}) \rangle}{\partial \rho} \,\, .
\end{equation}
%
Defining the fluctuation term as
%
\begin{equation}
R({\boldsymbol{\rho}}) =   \langle n({\boldsymbol{\rho}}) N \rangle  - \langle n({\boldsymbol{\rho}}) \rangle \langle N \rangle \,\, ,
\end{equation}
%
the temperature T can be derived from the slope of the linear equation
%
\begin{equation}
k_B T L({\boldsymbol{\rho}}) = R({\boldsymbol{\rho}})  \,\, .
\end{equation}
%
This is the proposal by Jason and Ho \cite{zhou1}.
\\ \\
%
\noindent For some reasons, we have to generalize their proposal. 
Interested only in fluctuations
%
\begin{equation}
R_\xi ({\boldsymbol{\rho}})  = \int d{\boldsymbol{\rho}}'  \left\{ \langle n({\boldsymbol{\rho}}) n({\boldsymbol{\rho}}') \rangle  - \langle n({\boldsymbol{\rho}}) \rangle \langle n({\boldsymbol{\rho}}') \rangle \right\} \theta ( \xi - \left| {\boldsymbol{\rho}} - {\boldsymbol{\rho}}' \right| )
\end{equation}
%
within a window size $\xi$, the linear equation
%
\begin{equation}
k_B T L({\boldsymbol{\rho}}) = R_\xi({\boldsymbol{\rho}})  \,\, .
\end{equation}
%
remains still valid only if the value of $\xi$ exceeds the density-density correlation length of the system. 
Here, $\theta(\cdot)$ is the Heaviside step function.
For radial symmetric lattices, both the dissipation and fluctuation terms could be averaged over the angular variable $\phi$, i.e.
%
\begin{eqnarray}
L(\rho) &=& \frac{1}{2\pi} \int  L ({\boldsymbol{\rho}})   \,  d\phi   \\
R(\rho) &=& \frac{1}{2\pi} \int R ({\boldsymbol{\rho}})   \,  d\phi
\end{eqnarray}
%
where ${\boldsymbol{\rho}} = {\boldsymbol{\rho}}(\rho,\phi)$.
Figure \ref{FD_1} illustrates the scheme, and highlights the importance of noise control with correct window sizing.
In the limit of infinite window size, i.e. $\xi \rightarrow \infty$, the noise in $R_\xi(\rho)$ grows beyond control which makes the FD thermometry scheme infeasible for accurate temperature determination. 
\\ \\
\noindent Finally, the temperature T can be estimated from a least-squares fit taking the measurement errors into account, i.e.
%
\begin{equation}
\left( \sum_i \frac{R_i^2}{\Delta_{Li}^2}  \right) -  \left( \sum_i \frac{L_i R_i}{\Delta_{Li}^2 } \right) \hat{T} =
- \left( \sum_i \frac{L_i R_i}{\Delta_{Ri}^2 } \right) \hat{T}^3 +  \left( \sum_i \frac{L_i^2}{\Delta_{Ri}^2}  \right)  \hat{T}^4 \,\, .
\end{equation} 
%
\\ \\
\noindent Note that the FD thermometry scheme reduces to the proposal by Jason and Ho in the limit of infinite window sizing, i.e. $\xi \rightarrow \infty$.


\subsection{Window sizing}


\begin{figure}
\centering
\includegraphics[width=10cm]{./figures/FD-2.jpg}
\caption{Illustration of the FD thermometry scheme by showing different $L(\rho)- R_\xi(\rho)$ plots at various window sizes $\xi=0,1,2,3,5,\infty$. This approach illustrates how the density-density correlation length can be found in an experimental system. When the window size is smaller than density-density correlation length, systematic errors set in, and this results in nonlinear $L(\rho)-R(\rho)$ behaviour, while for $\xi$ larger than the denisty-density correlation length, the behaviour of $L(\rho)-R(\rho)$ is linear. However, statistical noise also increases with increasing window size. The parameters are the same as in figure \ref{FD_1}.
}
\label{FD_2}
\end{figure}

Accurate determination of temperature relies heavily on the noise control of $R_\xi(\rho)$ with a correct window size of $\xi$.
Due to a lack of knowledge of the density-density correlation length, the optimal choice for $\xi$ is unknown at first which is to be determined through a systematic approach as illustrated in figure \ref{FD_2}.
In this approach, $\xi$ is gradually increased from 0 until it exceeds the density-density correlation length, indicated by the linearity of the $L(\rho)$-$R(\rho)$ plot.
Increasing $\xi$ further will smear out the linearity due to increasing statistical noise in $R(\rho)$.
This approach illustrates how the density-density correlation length could be found in an experimental system.
In this example, with a temperature of $T/t=1$, the optimal window size of $\xi=3$ is sufficient to capture (almost) all the correlations.
%
\begin{table} [!htp]
\centering
\begin{tabular}{| c | c | c |}
\hline
System & \multicolumn{2}{|c|}{nr of shots} \\
\cline{2-3}
       & $\xi$=3  &  $\xi$=$\infty$ \\
\hline
$U/t=10$ , $T/t=1$  &  20  &  $\mathcal{O}(10^4)$ \\
$U/t=10$ , $T/t=3$  &  14  &  $\mathcal{O}(10^4)$ \\
$U/t=50$ , $T/t=1$  &  21  &  $\mathcal{O}(10^4)$ \\
$U/t=50$ , $T/t=3$  &  12  &  $\mathcal{O}(10^4)$ \\
\hline
\end{tabular}
\caption{Number of uncorrelated measurements needed to determine the temperature accurately within 5\% error for 3D ${}^{87}$Rb optical lattice experiments trapping 125,000 bosons.   The variance reduction through window-sizing leads to orders of magnitudes improvement, and thus making FD thermometry scheme as a feasible tool for accurate temperature determination. The parameters are the same as in figure \ref{FD_1}.
}
\label{table:nr_shot1}
\end{table}
%
Using $\xi=3$ as the standard for $T/t=1$ and higher temperatures, we show the number of independent measurements needed to estimate the temperature within 5\% error in Table \ref{table:nr_shot1}. 
The enormous variance reduction through window sizing turns the FD thermometry scheme into a feasible tool for ultracold bosonic optical lattice experiments.

\begin{figure}[h]
\includegraphics[width=7cm]{./figures/FD-4.jpg}
\includegraphics[width=7cm]{./figures/FD-5.jpg}
\caption{The FD thermometry scheme remains valid over 20\% spread in temperature T (left), and 1\% spread in the particle number N (right). The parameters for the 3D bosonic $^{87}$Rb system are $N=125 000$, $U/t=10$, $T/t=1$, $V_T/t=0.0091$. Here, the optimal window size is $\xi=3$, and 20 independent measurements are obtained from a QMC simulation for accurate temperature determination within 5\% error.}
\label{FD_4}
\end{figure}
%
\begin{figure}[h]
\centering
\includegraphics[width=8cm]{./figures/FD-3.jpg}
\caption{FD thermometry scheme at slightly lower temperature, for a 3D bosonic $^{87}$-Rb optical lattice with parameters $N= 125 000$, $U/t=10$, $T/t=0.5$, $V_T/t = 0.0091$. Here, the optimal window size is $\xi=5$, and the data is obtained by averaging over 100 independent measurements from a QMC simulation.}
\label{FD_3}
\end{figure}
%
\noindent With 20 independent measurements uniformly distributed over 20\% spread in T and 1\% spread in N, the FD thermometry scheme remains applicable. See figure \ref{FD_4}.
At lower temperature, the scheme remains valid, although a larger $\xi$ is needed, and the statistical noise will inevitably grow.
See figure \ref{FD_3}.
In this example, with a temperature of $T/t=0.5$, the optimal window size is $\xi=5$, and 100 independent measurements are required to attain an accurate temperature estimate within 5\% error.





\subsection{Experimental imperfections}

\subsubsection{Density imaging with few-sites resolution}

Although we aim at single-site resolution detection tools in the analysis of the FD thermometry scheme in this chapter, the scheme remains applicable when the resolution is just a few sites, therefore giving some room for possible experimental imperfections.
%
\begin{figure}[h]
\centering
\includegraphics[width=14cm]{./figures/FD-6.jpg}
\caption{The FD thermometry scheme remains applicable to in-situ density experiments, which have a resolution of a few sites. The system consists of a 3D optical lattice with $N=125 000$ bosonic $^{87}$Rb atoms. Increasing the bin width increases the systematic error in the dissipation term $L(\rho)$, but the temperature estimate remains reliable. The parameters are the same as in figure \ref{FD_1}. Here, the optimal window size is $\xi=3$, and 1000 independent measurements are averaged out from a QMC simulation.}
\label{FD_6}
\end{figure}
%
See figure \ref{FD_6}.
Here, we simulate such experiments with different resolutions up to 5 sites, by varying the bin.
The FD thermometry scheme suffers from increasing systematic errors in the dissipation term $L(\rho)$ both for $U/t=10$ and for $U/t=50$ at $T/t=1$.
Yet, a relatively linear $L(\rho)$-$R(\rho)$ relationship could still be observed on average. For these cases alone, the estimated temperatures do not deviate more than 10\% as shown in Table \ref{table:resolution} .
However, when the resolution is worse than 5 sites, uncontrolled systematic errors dominate, and the scheme fails.

%
\begin{table} [!htp]
\centering
\begin{tabular}{| c | c | c |}
\hline
n-site resolution  & \multicolumn{2}{|c|}{estimated temperature}        \\
(binwidth)         & \multicolumn{2}{|c|}{{(units of t)}}  \\
\hline\hline
& ({\it U/t=10 , T/t=1}): & ({\it U/t=50 , T/t=1}): \\
\hline
1  &  0.977 $\pm$ 0.007  & 0.994 $\pm$ 0.008 \\
2  &  0.990 $\pm$ 0.006  & 1.014 $\pm$ 0.007 \\
3  &  0.997 $\pm$ 0.006  & 1.032 $\pm$ 0.007 \\
4  &  1.016 $\pm$ 0.007  & 1.045 $\pm$ 0.008 \\
5  &  1.048 $\pm$ 0.007  & 1.094 $\pm$ 0.008 \\
\hline
\end{tabular}
\caption{Estimated temperatures.}
\label{table:resolution}
\end{table}
%

\newpage 

\subsubsection{Density imaging with doublon-hole indistinguishability}

\begin{figure}
\centering
\includegraphics[width=10cm]{./figures/FD-7.jpg}
\caption{FD thermometry scheme in the presence of double-hole indistinguishability. The parameters are the same as in figure \ref{FD_1}. Blue circles (green squares) show the curve where doublons can (cannot) be distinguished from holes.}
\label{FD_7}
\end{figure}

Current single-site addressability experiments using fluorescence techniques can only measure the parity (even-odd) of the occupation number per site \cite{bloch1}.
This experimental imperfection affects the FD thermometry scheme.
Figure \ref{FD_7} illustrates the effect of double-hole indistinguishability.
\\ \\
\noindent An inspection of the $L(\rho)$-$R_\xi(\rho)$ plot shows a group of data points that fall on the linear slope.
These points correspond to the bosons at the edges.
Deep in the edges, the number of doublons is very low compared to the number of holes due to the low overall density and the high potential-energy cost of creating a doublon.
By selecting those points in the $L(\rho)$-$R_\xi(\rho)$ plot that are on the linear slope, we could still obtain an acceptable estimate of the temperature. 
\\ \\
\noindent The estimated temperatures for this example are tabulated as follows.


%
\begin{table} [h]
\centering
\begin{tabular} {| c | c |}
\hline
System              & estimated temperature \\
                    & (units of t)        \\
\hline
$U/t=10$ , $T/t=1 $ &  0.985 $\pm$ 0.008  \\
$U/t=50$ , $T/t=1 $ &  1.003 $\pm$ 0.012  \\
\hline
\end{tabular}
\label{table:fitted_temp_even_odd}
\end{table}
%

\section{Wing thermometry}

\subsection{Theory}

Deep in the edges, or {\it{wings}}, there will always be a normal region for any temperature $T/t$ and interaction strength $U/t$ where the system is well described by the second-order high temperature expansion (HTE2).
\\ \\
\noindent In this section, we shall formulate the theory. 
Note that a slightly different notation from the paper will be used in this section for the sake of convenience.
\\ \\
\noindent The partition function of the system is
%
\begin{equation}
Z = \mathrm{Tr} \, e^{-\beta(\hat{H} - \mu N)} = \mathrm{Tr} \, e^{-\beta (D + \hat{V})}
\end{equation}
%
where the diagonal/offdiagonal part reads
%
\begin{eqnarray}
D &=& \frac{U}{2} \sum_i  n_i (n_i - 1) - \sum_i \mu_i n_i    \\
\hat{V} &=& - t \sum_{\langle i , j \rangle}  \hat{b}_i^\dag \hat{b}_j \,\, .
\end{eqnarray}
%
Since
%
\begin{eqnarray}
\left[ \sum_k \mu_k n_k \, , \, \sum_{\langle i, j \rangle} \hat{b}_i^\dag \hat{b}_j \right] &=& \sum_k \sum_{\langle i, j \rangle} \mu_k \left[ \hat{b}_k^\dag \hat{b}_k \, , \, \hat{b}_i^\dag \hat{b}_j   \right]    \nonumber \\
 &=& \sum_k \sum_{\langle i, j \rangle} \mu_k \left\{   \hat{b}_i^\dag \left[ \hat{b}_k^\dag \hat{b}_k \, , \, \hat{b}_j   \right]  + \left[ \hat{b}_k^\dag \hat{b}_k \, , \, \hat{b}_i^\dag   \right] \hat{b}_j   \right\}   \nonumber \\
 &=& \sum_k \sum_{\langle i, j \rangle} \mu_k \left\{  - \hat{b}_i^\dag \hat{b}_k \delta_{kj}  +  \hat{b}_k^\dag \hat{b}_j \delta_{ki}  \right\}  \nonumber \\
 &=& \sum_{\langle i, j \rangle} \, (\mu_i - \mu_j) \, \hat{b}_i^\dag \hat{b}_j   
\end{eqnarray} 
%
\begin{eqnarray}
\left[ \sum_k n_k^2 \, , \, \sum_{\langle i, j \rangle} \hat{b}_i^\dag \hat{b}_j \right] 
 &=& \sum_k \sum_{\langle i, j \rangle} \left\{  n_k  \left[ \hat{b}_k^\dag \hat{b}_k \, , \, \hat{b}_i^\dag \hat{b}_j   \right]  +  \left[ \hat{b}_k^\dag \hat{b}_k \, , \, \hat{b}_i^\dag \hat{b}_j   \right] n_k \right\}   \nonumber \\
 &=& \sum_{\langle i, j \rangle} \left\{  (n_i - n_j) \hat{b}_i^\dag \hat{b}_j + \hat{b}_i^\dag \hat{b}_j  (n_i - n_j) \right\}  \nonumber \\
 &=& \sum_{\langle i, j \rangle} \hat{b}_i^\dag \hat{b}_j (n_i - n_j + 2)   \,\, ,
\end{eqnarray}
%
we derive the following important commutator
%
\begin{equation}
\left[ D \, , \hat{V} \right] = - t \sum_{\langle i, j \rangle}  \hat{b}_i^\dag \hat{b}_j \gamma_{ij} 
\end{equation}
%
where
%
\begin{equation}
\gamma_{ij} = U(n_i - n_j) + U - (\mu_i - \mu_j)  \,\, .
\end{equation}
%
The following are immediate consequences.
%
\begin{eqnarray}
\left[ D \, , \, \left[ D \, , \hat{V} \right]  \right]  &=&  - t \sum_{\langle i, j \rangle}  \hat{b}_i^\dag \hat{b}_j \gamma_{ij}^2   \\
\left[ D \, \left[ D \, , \, \left[ D \, , \hat{V} \right]  \right] \right]  &=& - t \sum_{\langle i, j \rangle}  \hat{b}_i^\dag \hat{b}_j \gamma_{ij}^3
\end{eqnarray}
%
Using the relation
%
\begin{equation}
e^{\tau D} \,  \hat{V} \, e^{-\tau D} = \hat{V} + \tau \left[ D \, , \, \hat{V} \right] + \frac{\tau^2}{2!} \left[ D \, , \, \left[ D \, , \hat{V} \right]  \right]   + \cdots   \,\, 
\end{equation}
%
we arrive at
%
\begin{equation}
e^{\tau D} \,  \hat{V} \, e^{-\tau D} =  - t \sum_{\langle i, j \rangle}  \hat{b}_i^\dag \hat{b}_j e^{\tau\gamma_{ij}}  \,\, .
\end{equation}
%
Next, we expand the partition function up to second order in $\hat{V}$, i.e.
%
\begin{equation}
Z = \mathrm{Tr} \, e^{-\beta D}  +  \int_0^\beta d\tau_1\int_0^{\tau_1} d\tau_2 \, \mathrm{Tr} \, e^{-\beta D} \, e^{\tau_1 D} \,  \hat{V} \, e^{-\tau_1 D}  \, e^{\tau_2 D} \,  \hat{V} \, e^{-\tau_2 D} \,\, .
\end{equation}
%


\subsubsection{Zeroth order}

Defining
%
\begin{equation}
Z^{(0)} = \mathrm{Tr} \, e^{-\beta D}   \,\, ,
\end{equation}
%
we have
%
\begin{equation}
Z^{(0)} = \sum_{\{n_i\}} \prod_i  e^{-\beta D_i} = \prod_i \sum_{n_i} e^{-\beta D_i} = \prod_i Z_i^{(0)}
\end{equation}
%
where
%
\begin{eqnarray}
Z_i^{(0)} &=& \sum_{n_i = 0}^\infty e^{-\beta D_i}   \\
D_i &=& \frac{U}{2} n_i (n_i -1) - \mu_i n_i
\end{eqnarray}
%
The zeroth-order density is therefore
%
\begin{equation}
\langle n_i^{(0)} \rangle  = \frac{1}{Z^{(0)}} \mathrm{Tr} \, n_i \, e^{-\beta D} = \frac{1}{Z_i^{(0)}} \sum_{n_i = 0}^\infty n_i \, e^{-\beta D_i}  \,\, .
\end{equation}
%


\subsubsection{Second order}

Expressing
%
\begin{equation}
Z \approx Z^{(0)} ( 1 + Z^{(2)})   
\end{equation}
%
such that
%
\begin{equation}
Z^{(2)} = \frac{1}{Z^{(0)}} \left(  \int_0^\beta d\tau_1\int_0^{\tau_1} d\tau_2 \, \mathrm{Tr} \, e^{-\beta D} \, e^{\tau_1 D} \,  \hat{V} \, e^{-\tau_1 D}  \, e^{\tau_2 D} \,  \hat{V} \, e^{-\tau_2 D} \right)  \,\, ,
\end{equation}
%
we arrive at
%
\begin{eqnarray}
Z^{(2)} &=& (-t)^2 \sum_{\langle i , j \rangle}  \frac{1}{Z_i^{(0)}}  \frac{1}{Z_j^{(0)}} \int_0^\beta d\tau_1\int_0^{\tau_1} d\tau_2   \left[ \sum_{n_i^\delta = n_i - 1} \sum_{n_j^\delta = n_j + 1}  n_i n_j^\delta e^{-\beta D_i} e^{-\beta D_j}   e^{\tau_1 \gamma_{ij}^\delta} e^{\tau_2 \gamma_{ji}}    \right.           \nonumber \\
& & + \left. \sum_{n_j^\delta = n_j - 1} \sum_{n_i^\delta = n_i + 1}  n_j n_i^\delta e^{-\beta D_i} e^{-\beta D_j}   e^{\tau_1 \gamma_{ji}^\delta} e^{\tau_2 \gamma_{ij}}  \right]  \,\, .
\end{eqnarray}
%
Next, we perform the integrals. If $\gamma_{ij}^\delta = \gamma_{ji} = 0$, 
%
\begin{equation}
 \int_0^\beta d\tau_1\int_0^{\tau_1} d\tau_2 \, e^{\tau_1 \gamma_{ij}^\delta} e^{\tau_2 \gamma_{ji}}  = \frac{\beta^2}{2}  \,\, ,
\end{equation}
%
otherwise,
%
\begin{equation}
\int_0^\beta d\tau_1\int_0^{\tau_1} d\tau_2 \, e^{\tau_1 \gamma_{ij}^\delta} e^{\tau_2 \gamma_{ji}}  = \beta^2 \left[  \frac{1 - e^{\beta \gamma_{ij}^\delta}}{ (\beta \gamma_{ij}^\delta )(\beta \gamma_{ji})}   - \frac{1 - e^{\beta \gamma_{ij}^\delta + \beta \gamma_{ji}}   }{ ( \beta \gamma_{ij}^\delta + \beta \gamma_{ji})(\beta \gamma_{ji})}       \right]
\end{equation}
%
Defining 
%
\begin{equation}
\Gamma_{ij}^\delta = \frac{1 - e^{\beta \gamma_{ij}^\delta}}{ (\beta \gamma_{ij}^\delta )(\beta \gamma_{ji})}   - \frac{1 - e^{\beta \gamma_{ij}^\delta + \beta \gamma_{ji}}   }{ ( \beta \gamma_{ij}^\delta + \beta \gamma_{ji})(\beta \gamma_{ji})}    \,\, ,
\end{equation}
%
we have
%
\begin{equation}
Z^{(2)} = \sum_{\langle i , j \rangle}  \frac{(-\beta t)^2}{ Z_i^{(0)} Z_j^{(0)} } \left[  \sum_{n_i n_j}^{(-,+)} n_i n_j^\delta e^{-\beta(D_i + D_j)} \Gamma_{ij}^\delta \, + \,  \sum_{n_i n_j}^{(+,-)} n_i^\delta n_j e^{-\beta(D_i + D_j)} \Gamma_{ji}^\delta    \right]
\end{equation}
%
where ($-,+$) stands for ($n_i^\delta = n_i - 1$, $n_j^\delta = n_j + 1$), and vice-versa for ($+,-$).
\\ \\
\noindent Since 
%
\begin{equation}
\langle n_i \rangle = -\frac{1}{\beta} \frac{\partial}{\partial \mu_i} \log Z \,\, ,
\end{equation}
%
we have
%
\begin{equation}
\langle n_i \rangle \approx -\frac{1}{\beta} \frac{\partial}{\partial \mu_i} \left[ \log Z^{(0)} + \log ( 1+ Z^{(2)} ) \right] \approx \langle n_i^{(0)} \rangle  - \frac{1}{\beta} \frac{\partial Z^{(2)}}{\partial \mu_i} \,\, .
\end{equation}
%
\begin{figure}[h]
\centering
\includegraphics[width=10cm]{./figures/HTE2-1.jpg}
\caption{Illustrating the concept of wing thermometry for a 3D bosonic ${}^{87}$Rb optical lattice system, {\it ie.,} describing the normal state by high temperature series expansions. Blue circles: In-situ density profile obtained from 100 uncorrelated measurements obtained by a QMC simulation with parameters
$U/t = 10$ , $T/t = 3$ , $N = 125,000$. The superfluid-normal phase boundary occurs at the density $\langle n \rangle=0.42$ or chemical potential $\mu / t = -2.75$.
The second order series captures all the physics in the normal regime, whereas the zeroth order has a very small validity range.} 
\label{HTE2_1}
\end{figure}
%
After differentiation, the density, up to the second order, is
%
\begin{eqnarray}
\langle n_i \rangle &=&  \langle n_i^{(0)} \rangle \, + \,  \sum_{\langle i , j \rangle}  \frac{(-\beta t)^2}{ Z_i^{(0)} Z_j^{(0)} }  \nonumber \\
& & \times \left[ \sum_{n_i n_j}^{(-,+)} \left( n_i -   \langle n_i^{(0)} \rangle  + \frac{\chi_{ij}^\delta}{\Gamma_{ij}^\delta}   \right)  n_i n_j^\delta e^{-\beta(D_i + D_j)} \Gamma_{ij}^\delta \right.  \nonumber \\
& & +  \left.  \sum_{n_i n_j}^{(+,-)} \left( n_i -   \langle n_i^{(0)} \rangle  + \frac{\chi_{ji}^\delta}{\Gamma_{ji}^\delta}   \right)  n_i^\delta n_j e^{-\beta(D_i + D_j)} \Gamma_{ji}^\delta    \right]
\end{eqnarray}
%
where
%
\begin{equation}
\chi_{ij}^\delta =  \frac{e^{\beta \gamma_{ij}^\delta}}{ (\beta \gamma_{ij}^\delta )(\beta \gamma_{ji})}   + \frac{1 - e^{\beta \gamma_{ij}^\delta + \beta \gamma_{ji}}   }{ ( \beta \gamma_{ij}^\delta + \beta \gamma_{ji})(\beta \gamma_{ji})^2} - \frac{(1 - e^{\beta \gamma_{ij}^\delta})(\beta \gamma_{ij}^\delta - \beta \gamma_{ji} )}{ (\beta \gamma_{ij}^\delta )^2(\beta \gamma_{ji})^2}  \,\, .
\end{equation}

%


\subsection{Advantages and limitations}


Figure \ref{HTE2_1} illustrates the concept of wing thermometry.
%
\begin{enumerate}
\item The zeroth order expansion does a poor job in describing the densities at the edges of the system, as also indicated in reference \cite{zhou2}.
\item The second order expansion gives a fairly accurate description of the densities over the edges if they exhibit normal behaviour. The second order expansion captures the entire physics in the normal regime.
\end{enumerate}
%
As illustrated in the example, the second order densities match almost exactly with the in-situ density profile, obtained from 100 independent measurements of a QMC simulation, over the entire normal region.
Here, the superfluid-normal phase boundary is determined from LDA calculations.
\\ \\
%
\begin{figure}[h]
\centering
\includegraphics[width=10cm]{./figures/HTE2-2.jpg}
\caption{Second order high temperature series expansion thermometry scheme for a bosonic system that is entirely in the normal phase. No more than a single shot of cross-sectional density is needed to estimate the temperature and chemical potential within 10\% accuracy. We take a 3D optical lattice system with bosonic ${}^{87}$Rb atoms and parameters $N = 125,000$, $\mu/t = 25.97$, $U/t = 50$, $T/t = 3$. The blue circles are obtained from a single measurement in a QMC simulation; and the red line is a least-square fit over the entire normal region where $\mu_{\rm fit}/t = 25.92$ and $T_{\rm fit}/t = 2.824$ nK.}
\label{HTE2_2}
\end{figure}
%
\noindent In the extreme case where the entire region is normal, the HTE2 thermometry scheme becomes an extremely powerful thermometer.
See figure \ref{HTE2_2}.
In fact, only a single measurement is sufficient to accurately estimate both the temperature $T$ and chemical potential $\mu$ up to 10\% error within the HTE2 thermometry scheme.
\\ \\
\noindent In the limit of very low temperatures where the normal wings become too narrow, it may be that the density in the edges is so low that it cannot be measured due to low signal-to-noise ratio. 
When this happens, HTE2 thermometry scheme becomes infeasible.
It is possible to use higher-order thermometry schemes, but the gain is minimal compared to the additional effort.
