# A Truly Second-Order and Unconditionally Stable Thermal Lattice Boltzmann Method
## 3 Unconditionally Stable Thermal Lattice Boltzmann Method (USTLBM)  
### Modifications on STLBM  
The initial formulations in the corrector step are shown in Equations (27)-(29). Apparently, Equation (27) recovers Equation (20) exactly. In the present work, Equations (28) and (29) are, respectively, replaced by the following formulas  
\[\rho(\mathbf{r},t)\mathbf{u}(\mathbf{r},t) =\rho^{*}u^{*}-\Big{(}1-\tfrac{1}{\tau_{\mathrm{r}}}\Big{)}{ \sum\limits_{\alpha}\mathbf{e}_{\alpha}\Big{[}f_{\alpha}^{neq}(\mathbf{r}+0.5 \mathbf{e}_{\alpha}\delta_{t},t-0.5\delta_{t})-f_{\alpha}^{neq}(\mathbf{r}-0.5 \mathbf{e}_{\alpha}\delta_{t},t-0.5\delta_{t})\Big{]}+\mathbf{F}_{E}\delta_{t}} \tag{32}\] \[\rho(\mathbf{r},t)e(\mathbf{r},t) =\rho^{*}e^{*}-\bigg{(}1-\frac{1}{\tau_{\mathrm{c}}}\bigg{)}{ \sum\limits_{\alpha}\Big{[}}g_{\alpha}^{neq}(\mathbf{r}+0.5\mathbf{e}_{\alpha }\delta_{t},t-0.5\delta_{t})-g_{\alpha}^{neq}(\mathbf{r}-0.5\mathbf{e}_{\alpha }\delta_{t},t-0.5\delta_{t})\Big{]} \tag{33}\]  
To present theoretical analysis on the accuracy, we start from non-equilibrium distribution functions involved in the above formulations. Performing Taylor-series expansion analysis on Equation (30) with third-order of truncated error gives  
\[ \begin{split}-\tau_{\mathrm{v}}\Big{[}f_{\alpha}^{eq}(\mathbf{r},t)-f_{\alpha}^{eq}(\mathbf{r}-\mathbf{e}_{\alpha}\delta_{t},t-\delta_{t}) \Big{]}&=-\tau_{\mathrm{v}}\delta_{t}\mathrm{D}f_{\alpha}^{eq} \Big{(}\mathbf{r}-\tfrac{1}{2}\mathbf{e}_{\alpha}\delta_{t},t-\tfrac{1}{2} \delta_{t}\Big{)}+O\big{(}\delta_{t}^{3}\big{)}\\ &=f_{\alpha}^{neq}\Big{(}\mathbf{r}-\tfrac{1}{2}\mathbf{e}_{ \alpha}\delta_{t},t-\tfrac{1}{2}\delta_{t}\Big{)}+O\big{(}\delta_{t}^{3}\big{)} \end{split}  \tag{34}\]  
Similarly, we have  
\[g_{\alpha}^{neq}\bigg{(}\mathbf{r}-\frac{1}{2}\mathbf{e}_{\alpha}\delta_{t},t -\frac{1}{2}\delta_{t}\bigg{)}=-\tau_{\mathrm{c}}\Big{[}g_{\alpha}^{eq}( \mathbf{r},t)-g_{\alpha}^{eq}(\mathbf{r}-\mathbf{e}_{\alpha}\delta_{t},t- \delta_{t})\Big{]}+O\Big{(}\delta_{t}^{3}\Big{)} \tag{35}\]  
Equations (34) and (35) indicate that, with desired order of accuracy, the approximated non-equilibrium distribution functions in Equations (30) and (31) are at space/time level of \((\mathbf{r}-0.5\mathbf{e}_{\alpha}\delta_{t},t-0.5\delta_{t})\). Such shifting in time and space favors our construction of a central-difference scheme both in space and in time, with which the second-order of accuracy can be guaranteed. Maintaining its location in time level, we further expand Equation (34) with respect to the space/time levels of \((\mathbf{r},t-0.5\delta_{t})\), which gives  
\[ \begin{split} f_{\alpha}^{neq}\Big{(}\mathbf{r}-\tfrac{1}{2} \mathbf{e}_{\alpha}\delta_{t},t-\tfrac{1}{2}\delta_{t}\Big{)}&=f_ {\alpha}^{neq}\Big{(}\mathbf{r},t-\tfrac{1}{2}\delta_{t}\Big{)}-\tfrac{1}{2} \delta_{t}\mathbf{e}_{\alpha}\cdot\
abla f_{\alpha}^{neq}\Big{(}\mathbf{r},t -\tfrac{1}{2}\delta_{t}\Big{)}\\ +\tfrac{1}{2}\Big{(}\tfrac{1}{2}\delta_{t}\Big{)}^{2}\Big{[} \mathbf{e}_{\alpha}\cdot\
abla f_{\alpha}^{neq}\Big{(}\mathbf{r},t-\tfrac{1}{ 2}\delta_{t}\Big{)}\Big{]}^{2}+O\big{(}\delta_{t}^{3}\big{)}\end{split}  \tag{36}\]  
\[ \begin{split} f_{\alpha}^{neq}\Big{(}\mathbf{r}+\tfrac{1}{2} \mathbf{e}_{\alpha}\delta_{t},t-\tfrac{1}{2}\delta_{t}\Big{)}&=f_ {\alpha}^{neq}\Big{(}\mathbf{r},t-\tfrac{1}{2}\delta_{t}\Big{)}+\tfrac{1}{2} \delta_{t}\mathbf{e}_{\alpha}\cdot\
abla f_{\alpha}^{neq}\Big{(}\mathbf{r},t -\tfrac{1}{2}\delta_{t}\Big{)}\\ +\tfrac{1}{2}\Big{(}\tfrac{1}{2}\delta_{t}\Big{)}^{2}\Big{[} \mathbf{e}_{\alpha}\cdot\
abla f_{\alpha}^{neq}\Big{(}\mathbf{r},t-\tfrac{1}{2} \delta_{t}\Big{)}\Big{]}^{2}+O\big{(}\delta_{t}^{3}\big{)}\end{split}  \tag{37}\]  
Substituting Equations (36) and (37) into Equation (32) leads to  
\[