## Appendix A Transfer Matrix Method  
In this Appendix, the method to calculate the stability limit is described [14]. Equations (9) and (10) are modified in a matrix form as follows:  
\[\frac{d}{dx}\left(\begin{array}{c}P(x)\\ U(x)\end{array}\right)=C(x)\left(\begin{array}{c}P(x)\\ U(x)\end{array}\right), \tag{10}\]  
\[C(x)\equiv\left(\begin{array}{cc}0&-\frac{i\omega p_{m}}{1-\lambda x}\\ -\frac{i\omega[1+(\gamma-1)\chi_{\kappa}]}{\gamma T_{m}}&\frac{\chi_{\kappa}- \chi_{\kappa}}{(1-\lambda x)(1-\lambda y_{\Gamma})}\frac{1}{T_{m}}\frac{dT_{m }}{dx}\end{array}\right).\]  
A forward different scheme using the fourth-order Runge-Kutta method to Equation (10) is applied to obtain the following expressions:  
\[\left(\begin{array}{c}P(x+\Delta x,t)\\ U(x+\Delta x,t)\end{array}\right)=\left(E+\Delta xC^{\prime}(x)\right)\left( \begin{array}{c}P(x)\\ U(x)\end{array}\right), \tag{11}\]  
\[C^{\prime}(x) = \frac{1}{6}\left(RK_{A}+2 RK_{B}+2 RK_{C}+RK_{D}\right),\] \[RK_{A} = C(x),\] \[RK_{B} = C(x+\Delta x/2)\left(E+\frac{\Delta x}{2}RK_{A}\right),\] \[RK_{C} = C(x+\Delta x/2)\left(E+\frac{\Delta x}{2}RK_{B}\right),\] \[RK_{D} = C(x+\Delta x)\left(E+\Delta xRK_{C}\right),\]  
where \(E\) denotes a unit matrix. When the value of \(\frac{dT_{m}}{dx}\) is known and a flow channel is uniform, we can numerically integrate Equation (11) and obtain the following expression:  
\[\left(\begin{array}{c}P(x)\\ U(x)\end{array}\right)=M_{\Pi}(x,x_{0})\left(\begin{array}{c}P_{0}(x_{0})\\ U_{0}(x_{0})\end{array}\right), \tag{12}\]\[M_{\Pi}(x,x_{0})\equiv(E+\Delta xC^{\prime}_{n-1})(E+\Delta xC^{\prime}_{n-2})\ldots (E+\Delta xC^{\prime}_{1})(E+\Delta xC^{\prime}_{0}).\]  
Here, \(n\) denotes the number of partitions between \(x_{0}\) and \(x\), \(\Delta x\) is defined as \((x-x_{0})/n\), and \(C^{\prime}_{j}\) represents \(C^{\prime}\) at \(x=x_{0}+j\Delta x\).  
The thermoacoustic cooler shown in Figure 1 was divided into the ten components. Their names and numbers are defined as follows: (1) the engine stack, (2) the hot heat exchanger, (3) the thermal buffer tube A, (4) the waveguide A, (5) the ambient heat exchanger A, (6) the cooler stack, (7) the cold heat exchanger, (8) the thermal buffer tube B, (9) the waveguide B, and (10) the ambient heat exchanger B. The numerical integration along each component was performed and their transfer matrices were obtained. The total flow-path areas in the components are different from each other, and thus the connecting matrix is as follows:  
\[O_{k,l}=\left(\begin{array}{cc}1&0\\ 0&A_{l}/A_{k}\end{array}\right), \tag{10}\]  
where \(A\) denotes the total-path area and the subscripts \(k\) and \(l\) denote the number of the components. The transfer matrices of the components and the connecting matrices are used to express the transfer matrix of the total system as follows:  
\[M_{all}=M_{10}O_{10,9}M_{9}O_{9,8}M_{8}O_{8,7}M_{7}O_{7,6}M_{6}O_{6,5}M_{5}O_{5,4}M_{4}O_{4,3}M_{3}O_{3,2}M_{2}O_{2,1}M_{1}. \tag{11}\]  
Additionally, \(M_{all}\) was used such that the pressure \(P_{e,A}\) and velocity \(U_{e,A}\) at the ambient end of the engine stack corresponds to the same as follows:  
\[M_{all}\left(\begin{array}{c}P_{e,A}\\ U_{e,A}\end{array}\right)=\left(\begin{array}{c}P_{e,A}\\ U_{e,A}\end{array}\right). \tag{12}\]  
The solution \((P_{e,A}\), \(U_{e,A})\) of Equation (12) is nonzero if the determinant of the matrix \((M_{all}-E)\) is zero, i.e., if the following expression is applicable:  
\[(m_{11}-1)(m_{22}-1)-m_{12}m_{21}=0, \tag{13}\]  
where \(E\) denotes the unit matrix and \(m_{ij}\) denotes the element of \(M_{all}\). Therefore, Equation (13) was numerically solved to achieve the condition of the stability limit of the spontaneous gas oscillation induced in the thermoacoustic cooler.