It is rather easy to verify numerically, that $\phi(x) = \exp(-c \vert x \vert^{\frac{1}{2(1-H)}} )$ does not satisfy $T_H \phi = \phi$.
To this end it is enough to observe that $(T_H \phi_c)(x) = c^{-2 (1-h)^2} (T_H \phi_1)(x c^{2 h(1-h)})$, which is obtained via changing variables $y \to c^{-2(1-h)} y$. This scaling tells you that if $T_H \phi_c = \phi_c$, then $T_H \phi_1 = c^{-2(1-h)^2} \phi_1( x c^{2 h (1-h)} )$.
Thus it suffices to show that functional forms of $T_H \phi$ and $\phi$ are different for your purported function. So we set $c=1$ and change variables into $y \to t^{2(1-h)}$. This gives
$ (T_H \phi)(x) = 2(1-h) \sqrt{\frac{2}{\pi}} \int_0^\infty t^{2(1-h)^2} \mathrm{e}^{-t} \exp(-\frac{x^2}{2} t^{-4 h (1-h)} ) \frac{\mathrm{d} t}{t} $ The integral can be evaluate in terms of Fox H-function. Specifically, in terms of Mellin-Barnes integral:
$ (T_H \phi)(x) =\sqrt{\frac{2}{\pi}} \frac{z^{2 \rho_h}}{2h} \frac{1}{2 \pi i} \int_{\gamma - i \infty}^{\gamma+i \infty} \left( \Gamma(s) \Gamma( -\rho_h + \frac{s}{\kappa_h} ) z^{-\frac{2s}{\kappa_h }} \right)\mathrm{d}s $ where $\gamma$ is an arbitrary $\gamma > \rho_h \kappa_h$ and $z = \frac{x}{\sqrt{2}}$ and we introduced $\kappa_h = 4 h(1-h)$ and $\rho_h = \frac{1-h}{2 h}$.
The Mellin-Barnes integral can be evaluated as the sum over poles of $\Gamma$ functions at $s= -k$ for $k \in \mathbb{Z}^+$ and at $s = \rho_h \kappa_h - k \, \kappa_h$, giving $ \begin{eqnarray} (T_H \phi)(x) = \frac{z^{2 \rho_h} }{2 h} \sqrt{\frac{2}{\pi}} && \sum_{k=0}^\infty (- z^{\frac{2}{\kappa_h}} )^k \frac{\Gamma(-\rho_h - \kappa_h \, k )}{k!} + \\ \sqrt{\frac{4 \rho_h \kappa_h}{\pi}} && \sum_{k=0}^\infty (- z^{ 2 } )^k \frac{\Gamma(\rho_h \kappa_h - k \, \kappa_h )}{k!} \end{eqnarray} $
When $h = \frac{1}{2}$, $\kappa_{1/2} = 1$ and $\rho_{1/2} = \frac{1}{2}$. Then the sums evaluate exactly:
$ \begin{eqnarray} T_{1/2} \phi &=& \frac{x}{\sqrt{2}} \sqrt{ \frac{2}{\pi} } \sum_{k=0}^\infty (-\frac{x^2}{2})^k \frac{\Gamma(-\frac{1}{2}-k)}{k!} + \sqrt{ \frac{2}{\pi} } \sum_{k=0}^\infty (-\frac{x^2}{2})^k \frac{\Gamma(1/2-k)}{k!} \\ &=& -\sqrt{2} \sinh( \sqrt{2} x) + \sqrt{2} \cosh( \sqrt{2} x) = \sqrt{2} \, \exp( -\sqrt{2} x) \end{eqnarray} $
It is also possible to evaluate these sums for rational values of $h$, but expressions involve hypergeometric functions and stop being so simple.
Added: Below we construct the $T_H \psi$, where $\psi(x)$ is a general Fox H-function. The class of Fox H-functions is a natural place to search for solutions of $T_H \psi = \psi$ because operator $T_H$ maps a Fox H-function into another Fox H-function. Specifically, let
$ \begin{eqnarray} \psi(x) &=& H^{m,n}_{p,q}\left( \left. \begin{array}{c} (a_1, \alpha_1), \ldots, (a_p, \alpha_p) \\ (b_1, \beta_1), \ldots, (b_q, \beta_q) \end{array} \right\vert x\right) \\ &=& \frac{1}{2 \pi i} \int_{\gamma-i \infty}^{\gamma + i \infty} \frac{ \prod_{k=1}^m \Gamma(1-a_k - \alpha_k s) \prod_{k=1}^n \Gamma(b_k + \beta_k s) }{ \prod_{k=m+1}^p \Gamma(1-a_k - \alpha_k s) \prod_{k=n+1}^q \Gamma(b_k + \beta_k s) } x^{-s} \mathrm{d} s \end{eqnarray} $
Then
$ (T_H \psi)(x) = \frac{1}{\sqrt{2 \pi}} \frac{1}{h} \left( \frac{x^2}{2} \right)^{\frac{(1-h)}{2h} } H^{m,n+1}_{p,q+1}\left( \left. \begin{array}{c} (a_1, \alpha_1), \ldots, (a_p, \alpha_p) \\ \left( -\frac{1-h}{2h}, \frac{1}{2h}\right), (b_1, \beta_1), \ldots, (b_q, \beta_q) \end{array} \right\vert \left(\frac{x^2}{2}\right)^{\frac{1}{2h}} \right) $
Now let's understand the mechanics of how it works for $h=\frac{1}{2}$. Your function $\phi(x)$ with $h=\frac{1}{2}$ is $\exp(-x) = H^{0,1}_{0,1}(x)$ with $\beta_1=1$ and $b_1 = 0$, so
$ \begin{eqnarray} \left. (T_H \phi)(x) \right\vert_{h=\frac{1}{2}} &=& \frac{2}{\sqrt{2 \pi}} \left( \frac{x^2}{2} \right)^{\frac{1}{2} } H^{0,2}_{0,2}\left( \left. \begin{array}{c} - \\ \left( -\frac{1}{2}, 1\right), (0, 1) \end{array} \right\vert \frac{x^2}{2} \right) \\ &=& \sqrt{\frac{2}{\pi}} \left( \frac{x^2}{2} \right)^{\frac{1}{2} } \frac{1}{2 \pi i} \int_{\gamma-i \infty}^{\gamma+i \infty} \Gamma\left( -\frac{1}{2} + s\right)\Gamma(s) \left( \frac{x^2}{2} \right)^{-s} \mathrm{d} s \\ &=& \sqrt{\frac{2}{\pi}} \left( \frac{x^2}{2} \right)^{\frac{1}{2} } \frac{1}{2 \pi i} \int_{\gamma-i \infty}^{\gamma+i \infty} \sqrt{\pi} \Gamma(2s-1) 4^{1-s} \left( \frac{x^2}{2} \right)^{-s} \mathrm{d} s \left. = \right|_{s \to \frac{s}{2}} \\ &=& 2 x \frac{1}{2 \pi i} \int_{\gamma-i \infty}^{\gamma+i \infty} \Gamma(s-1) \left( \sqrt{2} x \right)^{-s} \mathrm{d} s \\ &=& (2 x) \left( \frac{\exp(- \sqrt{2} x)}{\sqrt{2} x} \right) = \sqrt{2} \exp(- \sqrt{2} x) \end{eqnarray} $
Notice that the use of tuplication identity for $\Gamma(x)$ was essential here, this essentially proves that within the family of Fox H-functions the solution to $(T_H \psi)(x) = c_1 \psi(c_2 x)$ for some constants $c_1$ and $c_2$ is only possible for $h=\frac{1}{2}$.
Now, the family of Fox H-functions is big and encompasses all hyper-geometric elementary functions, i.e. those that satisfy a linear differential equations with coefficients polynomial in $x$.
If you would like to find the eigenfunction of $T_h$ operator, I would suggest playing with quadrature rules. Suppose $(T_H f)(x) = \sum_{i=1}^n \omega_i K(x, y_i) f(y_i)$. Solve the system of equations
$ \sum_{i=1}^n \omega_i K(x_j, y_i) f(y_i) = f(x_j) $ numerically, which would allow you to explore the eigen-function for $h \not= \frac{1}{2}$.