## 2 Theoretical Model  
The model of the nonlocal nonlinear medium with complex potentials can be established by the normalized 2D coupled nonlinear Schrodinger equation as follows [24; 26; 36], with the transverse coordinates \(x\) and \(y\) scaled to the input beam width, \(w_{0}\), and the longitudinal coordinate \(z\) scaled to the diffraction length, \(L_{0}=2 kw_{0}^{2}\):  
\[\left\{ \begin{array}{l}iL_{z}+\
abla_{\perp}U+(V+iW)U+\sigma nlU=0\\ d\
abla_{\perp}n-n+|U|^{2}=0\end{array} \right. \tag{1}\]  
where \(\
abla_{\perp}=\partial_{xx}+\partial_{yy}\) indicates the 2D Laplace operator. This model can be realized in the liquid crystal E7 [37]. The dimensionless light field amplitudes \(U_{n}\) and the nonlinear refractive index \(n\) are scaled to \(2E\sqrt{\varepsilon_{0}\Delta\varepsilon_{E}}/\sqrt{\pi}\Delta\varepsilon_{L }k_{0}w_{0}\) and \(\varepsilon_{0}/k_{0}^{2}w_{0}^{2}\Delta\varepsilon_{E}\), respectively, where \(k_{0}\) is the wave numbers in a vacuum, \(\varepsilon_{0}\) is the vacuum permittivity, and \(E\) is the electric field strength. The material parameters of liquid crystals, namely \(\Delta\varepsilon_{E}\), \(\Delta\varepsilon_{L}\), and \(K\), correspond to the anisotropy at the frequency of the quasistatic electric field, the anisotropy at the light frequency, and the elastic constant, respectively. \(d\) is the degree of nonlocality, defined as \(d^{2}=\pi K/2w_{0}^{2}E\Delta\varepsilon_{E}\). \(\sigma=\pm 1\) represents the focusing or defocusingnonlinearity. Here, we fix \(\sigma=-1\), which corresponds to the defocusing type. \((V+iW)\) represents the PT-symmetric complex potential. \(V\) and \(W\) are the real part and the imaginary part, respectively. To simulate the triangular lattices, we set the PT complex as:  
\[\left\{ \begin{array}{l}V=2V_{0}\Big{[}\frac{3}{2}+\cos\!\left(k_{0}x+\frac{k_ {0}}{\sqrt{3}}y\right)+\cos\!\left(k_{0}x-\frac{k_{0}}{\sqrt{3}}y\right)+\cos \!\left(\frac{2k_{0}}{\sqrt{3}}y\right)\!\Big{]}\\ W=\frac{4}{3}W_{0}\Big{[}\sin\!\left(k_{0}x+\frac{k_{0}}{\sqrt{3}}y\right)+\sin \!\left(k_{0}x-\frac{k_{0}}{\sqrt{3}}y\right)+\sin\!\left(\frac{2k_{0}}{\sqrt {3}}y\right)\!\Big{]}\end{array} \right. \tag{2}\]  
where \(V_{0}\) is the depth of the PT-symmetric lattices and \(W_{0}\) is the gain/loss coefficient. These triangular lattices can be fabricated by interfering three ordinarily polarized broad beams [36]. No linear mode can exit via the bandgaps, but gap solitons exist there. Thus, we mainly concentrate on the nonlinear propagation modes in bandgaps. In addition, the threshold of PT symmetry breaking is found to be \(W_{0 th}=3\), at which the first bandgap and the infinite bandgap merge together. We substitute \(U(x,y,z)=q(x,y)e^{i\mu z}\) in Equation (1), and thus it becomes:  
\[\left\{ \begin{array}{l}-\mu q+\
abla_{\perp}q+(V+iW)q+\sigma nq=0\\ d\
abla_{\perp}n-n+\big{|}q\big{|}^{2}=0\end{array} \right. \tag{3}\]  
where \(q(x,y)\) is a complex function and \(\mu\) is the propagation constant. Equation (3) can be numerically solved by the developed modified squared-operator method (MSOM). [38] The power of solitons is defined as \(P=\int\int_{-\infty}^{+\infty}\!\big{|}q\big{|}^{2}dxdy\). We utilize the symmetric split-step Fourier method to simulate the propagation of solitons in the nonlinear triangular lattices. As for the linear stability, we add a small perturbation into the solitary solution:  
\[U(x,y,z)=\Big{[}q(x,y)+F(x,y)e^{\lambda z}+G^{*}(x,y)e^{\lambda^{*}z}\Big{]} \tag{4}\]  
where \(|F|,|G|\ll 1\) is a small perturbation, \(\lambda\) is the complex instability growth rate, and the superscript \"\(\ast\)\" represents the complex conjugation. By substituting Equation (4) into Equation (1) and reducing the equations, we obtain the eigenvalue equations:  
\[\left\{ \begin{array}{l}\lambda F=i[-\mu F+\
abla_{\perp}F+(V+iW)F+\sigma nG +\sigma q\Delta n]\\ \lambda G=i[\mu G-\
abla_{\perp}G-(V-iW)G-\sigma nF-\sigma q^{*}\Delta n] \end{array} \right.