# A Truly Second-Order and Unconditionally Stable Thermal Lattice Boltzmann Method
## 3 Unconditionally Stable Thermal Lattice Boltzmann Method (USTLBM)
### Stability Analysis
For illustration purposes, we take the term \([{\partial y_{2}}^{*}/{\partial y_{2}}]^{n}\) to elaborate the analyzing process:  
\[\left[\frac{\partial y_{2}^{*}}{\partial y_{2}}\right]^{n}=\sum_{\alpha}\left[ \frac{w_{a}e_{ax}^{2}}{c_{s}^{2}}+\frac{w_{a}e_{axy}^{3}y_{2\alpha}^{n}}{c_{s} ^{4}y_{1\alpha}^{n}}+\frac{w_{\alpha}e_{ax}^{2}e_{ay}y_{3\alpha}^{n}}{c_{s}^{4 }y_{1\alpha}^{n}}-\frac{w_{\alpha}e_{ax}y_{2\alpha}^{n}}{c_{s}^{2}y_{1\alpha}^ {n}}\right] \tag{53}\]  
In von-Neumann stability analysis, the variables are approximated by the Fourier serial expansions [2]. Consider a random component of the variables encountered here, which can be written as  
\[\mathrm{d}y_{i}^{*}=\overline{\mathrm{d}y}_{i}^{*}e^{i\mathbf{k}\cdot\mathbf{r}} \text{; }\mathrm{d}y_{i}^{n}=\overline{\mathrm{d}y}_{i}^{n}e^{i\mathbf{k}\cdot \mathbf{r}}\text{; }y_{i\alpha}^{n}=y_{i}^{n}e^{i\mathbf{k}\cdot(\mathbf{r}-\mathbf{e}_{ \alpha}\delta_{i})},i=1-4 \tag{54}\]  
In addition, note that the following isotropic properties of the lattice tensors are generally applicable within LBM framework:  
\[ \begin{array}{l}\sum\limits_{\alpha}w_{a}e_{ax}^{2}e_{ay}=\sum\limits_{ \alpha}w_{a}e_{ax}e_{ay}^{2}=\sum\limits_{\alpha}w_{a}e_{ax}^{3}=\sum\limits_ {\alpha}w_{a}e_{ay}^{3}=0\\ \sum\limits_{\alpha}w_{a}e_{ax}^{2}=\sum\limits_{\alpha}w_{a}e_{ay}^{2}=c_{s} ^{2}\\ \sum\limits_{\alpha}w_{a}e_{ax}e_{ay}=0\\ \sum\limits_{\alpha}w_{\alpha}=1\\ \sum\limits_{\alpha}w_{a}e_{ax}=\sum\limits_{\alpha}w_{a}e_{ay}=0\\ \end{array}  \tag{55}\]  
Substituting Equations (54) and (55) into Equation (53), we can derive that  
\[ \begin{array}{l}\left[\frac{\partial y_{2}^{*}}{\partial y_{2}}\right]^{n} =\frac{1}{c_{s}^{2}}\sum\limits_{\alpha}w_{\alpha}e_{ax}^{2}+\frac{y_{2}^{n} }{c_{s}^{2}y_{1\alpha}^{n}}\sum\limits_{\alpha}w_{\alpha}e_{ax}^{3}+\frac{y_{ 3}^{n}}{c_{s}^{3}y_{1\alpha}^{n}}\sum\limits_{\alpha}w_{\alpha}e_{ax}^{2}e_{ay }+\frac{y_{2}^{n}}{c_{s}^{2}y_{1\alpha}^{n}}\sum\limits_{\alpha}w_{\alpha}e_{ ax}\\ =1\end{array}  \tag{56}\]  
Following similar processes shown in Equations (53)-(56), we can find the values of all terms inside the brackets in Equation (52) (i.e., \([{\partial y_{1}}^{*}/{\partial y_{1}}]^{n}\)); and the equation can be simplified as  
\[\overline{\mathrm{d}\mathbf{Y}}^{*}=\overline{\mathrm{d}\mathbf{Y}}^{n} \tag{57}\]  
In the corrector step, it is shown in Equation (42) that no correction is performed on the density. Therefore, we can simply derive that  
\[\overline{\mathrm{d}y}_{1}^{n+1}=\overline{\mathrm{d}y}_{1}^{*}=\overline{ \mathrm{d}y}_{1}^{n} \tag{58}\]  
Next, we still take the variable \(y_{2}\) to illustrate the analyzing process in the corrector step. According to Equations (40), (43), and (45), we have  
\[y_{2}^{n+1}=y_{2}^{*}+(\tau_{v}-1)\sum\limits_{\alpha}e_{ax}f_{\alpha}^{eq}( \mathbf{r}+\mathbf{e}_{\alpha}\delta_{t},t)-(\tau_{v}-1)y_{2}^{*} \tag{59}\]  
Substitution of the equilibrium distribution based on the intermediate macroscopic variables into the second term on RHS of the above equation gives a nonlinear algebraic formula of  
\[\begin{array}{l}S=\sum\limits_{\alpha}e_{ax}f_{\alpha}^{eq}(\mathbf{r}+ \mathbf{e}_{\alpha}\delta_{t},t)\\ =\sum\limits_{\alpha}e_{ax}w_{a}y_{1-\alpha}^{*}+\sum\limits_{\alpha}e_{ax} \frac{w_{a}}{c_{s}^{2}}\left(e_{ax}y_{2-\alpha}^{*}+e_{ay}y_{3-\alpha}^{*}\right) \\ +\sum\limits_{\alpha}e_{ax}\frac{w_{a}}{where the subscript \"\(-\alpha\)\" represents the space level of \((\mathbf{r}+\mathbf{e}_{\mathrm{a}}\delta_{t})\). Similar to the procedures conducted in the corrector step, here we linearize the formula by differentiating it, which gives  
\[\mathrm{d}S=\left(\frac{\partial S}{\partial y_{1}}\right)^{*}\mathrm{d}y_{1}^{* }+\left(\frac{\partial S}{\partial y_{2}}\right)^{*}\mathrm{d}y_{2}^{*}+\left( \frac{\partial S}{\partial y_{3}}\right)^{*}\mathrm{d}y_{3}^{*} \tag{61}\]  
Following similar processes shown in Equations (53)-(56), the values of terms inside brackets can be evaluated.