First, I'll discuss the case of empty $C$.

Observe that the structures $(\varepsilon_{\alpha};+,\cdot,\mathsf{exp})$ are definitionally equivalent with the structures $HF(\alpha;<)$ (the structure of hereditarily finite sets with ordinals $<\alpha$ as urelements). Namely we consider the bijection $f\colon HF(\alpha)\to \varepsilon_{\alpha}$ define by recursion on $HF(\alpha)$:

- $f\colon \beta\longmapsto \varepsilon_{\beta}$;
- $f\colon \{a_1,\ldots,a_n\}\longmapsto 2^{f(a_1)}+\ldots+2^{f(a_n)}$, where $f(a_1)>\ldots>f(a_n)$.

Using the high expressive power of $HF(\alpha;<)$ it is easy to show that we could define the predicates $f(x)+f(y)=f(z)$, $f(x)f(y)=f(z)$, and $f(x)^{f(y)}=f(z)$. And in $(\varepsilon_{\alpha};+,\cdot,\mathsf{exp})$ we could easily define the predicates

- $\mathsf{Ur}(f^{-1}(x))$ as "$x$ is an $\varepsilon$-number";
- $f^{-1}(x)\in f^{-1}(y)$ as "$y$ isn't an $\varepsilon$-number and $2^x$ is in $2$-base Cantor normal form for $y$" or equivalently as "$y$ isn't an $\varepsilon$-number and $2^{x+1}z+2^{x}\le y< 2^{x+1}z+2^{x+1}$, for some $z$";
- $f^{-1}(x)<f^{-1}(y)$ as "$x$ and $y$ are $\varepsilon$-numbers and $x<y$".

Let us denote as $\mathsf{U}_{\alpha}(C)$ the elementary theory of $HF(\alpha;<,\langle c_{\beta}\mid \beta\in C\rangle)$, where $c_{\beta}$ is the constant for the ordinal $\beta$. The defined equivalence gives us a one-to-one correspondence $h$ between the models of $\mathsf{T}_{\varepsilon_{\alpha}}(\emptyset)$ and models of $\mathsf{U}_{\varepsilon_{\alpha}}(\emptyset)$. There are two important groups of models of $\mathsf{U}_{\alpha}(C)$: those of the form $HF(A)$, for some linear order $A$ and those that have ill-founded membership relation. The models of the second group are akin of non-standard models of arithmetic, but the models of the first group constitute essentially different phenomenon. Namely it is easy to see that models of $\mathsf{T}_{\varepsilon_{\alpha}}(\emptyset)$ of the form $h^{-1}(HF(A))$ would always poses a well-founded initial segment of the order type $\varepsilon_0$ (latter I'll show that provided that $\alpha\ge \omega^{\omega}$ the model $h^{-1}(HF(A))$ always start with the initial fragment of the order type $\varepsilon_{\omega^{\omega}}$). On the other hand the models of $\mathsf{T}_{\varepsilon_{\alpha}}(\emptyset)$ of the form $h^{-1}(\mathfrak{M})$, where $\mathfrak{M}$ have ill-founded membership predicate would always start with the initial segment of the order type $\omega+(\omega^{*}+\omega)A$, where $A$ is a dense linear order without end-points (here $\omega+(\omega^{*}+\omega)A$ is the order type of the non-standard model of arithmetic produced from $h^{-1}(\mathfrak{M})$ by the standard interpretation of arithmetic in hereditarily finite sets). But, I couldn't say much more about the order types of the models of $\mathsf{T}_{\varepsilon_{\alpha}}(\emptyset)$ of the latter group.

However the models of $\mathsf{U}_{\alpha}(\emptyset)$ of the form $HF(A)$ actually are easy to describe. There is a classical result of Ehrenfeucht that structures $HF(\alpha;<)$ and $HF(\beta;<)$ are elementary equivalent iff either

- $\alpha=\beta$ or
- $\alpha,\beta$ and are of the form $\alpha=\omega^\omega(1+\alpha')+\gamma$ and $\beta=\omega^\omega(1+\beta')+\gamma$.

By essentially the same argument as Ehrenfeucht's it is easy to show that for any $\alpha<\omega^{\omega}$ the structures $HF(\omega^{\omega}(1+A)+\alpha)$ and $HF(\omega^{\omega}(1+A')+\alpha)$ are elementary equivalent for any linear orders $A$ and $A'$. That for any $\alpha<\omega^{\omega}$ and any linear order $B$ if $HF(B)\equiv HF(\omega^{\omega}(1+A)+\alpha)$, then $B$ is of the form $\omega^{\omega}(1+A')+\alpha$. And that for $\alpha<\omega^{\omega}$ for any $A$ if $HF(A)\equiv HF(\alpha)$ then $A$ have the order type $\alpha$.

When we go back to models of $T_{\varepsilon_{\alpha}}(\emptyset)$ this shows that all the models of $T_{\varepsilon_{\omega^{\omega}(1+\beta)+\alpha}}(\emptyset)$ from this group are of the order types $\varepsilon_{\omega^{\omega}(1+A)+\alpha}$, where $A$ is some linear order. Here the order $\varepsilon_A$, for a given linear order $A$, is defined as the natural order on nested (formal) Cantor normal forms built from $0$ and constants $\varepsilon_a$, for all $a\in A$.

The addition of constants doesn't change the picture too much.

For each set of ordinals $C\subseteq \varepsilon_{\alpha}$ the structure $(\varepsilon_{\alpha};+,\cdot,\mathsf{exp},\langle c_{\beta}\mid \beta\in C\rangle)$ is definitionally equivalent to the structure $HF(\alpha;<,\langle c_{\beta}\mid \beta \in \mathbb{H}(C)\rangle)$, where $\mathbb{H}(C)$ is the least set of ordinals such that each $\beta\in C$ is the value $t(\varepsilon_{\beta_1},\ldots,\varepsilon_{\beta_n})$, for some $\beta_1,\ldots,\beta_n\in \mathbb{H}(C)$ and $t(x_1,\ldots,x_n)$ is a term build of operations $+,\cdot,\mathsf{exp}$. Further observe that in $HF(\alpha;<)$ (and actually even $(\alpha,<)$) given a constant $c_{\beta}$ for the ordinal $\beta$ with the Cantor normal form $\omega^{\omega+\beta_1}+\ldots+\omega^{\omega +\beta_n}+\omega^{k_1}+\ldots+\omega^{k_n}$ any other ordinal in the interval $[\omega^{\omega+\beta_1}+\ldots+\omega^{\omega +\beta_n},\omega^{\omega+\beta_1}+\ldots+\omega^{\omega +\beta_n}+\omega^{\omega})$ is definable in $HF(\alpha;<,c_{\beta})$. Henceforth, over $HF(\alpha;<)$ the constants for ordinals from a set $C\subseteq \alpha$ are definitionally equivalent to the set of constants for ordinals from $\{\omega^{\omega}(1+\beta)\mid \beta\in \mathbb{B}(C)\}$, where $\mathbb{B}(C)$ is the least set of ordinals such that each $\beta\in C$ either $<\omega^{\omega}$ or there is $\gamma\in \mathbb{B}(C)$ such that $\omega^{\omega}(1+\gamma)\le \beta<\omega^{\omega}(1+\gamma+1)$. Therefore $(\varepsilon_{\alpha};+,\cdot,\mathsf{exp},\langle c_{\beta}\mid \beta\in C\rangle)$ is definitionally equivalent to the structure $HF(\alpha;<,\langle c_{\omega^{\omega}(1+\beta)}\mid \beta\in \mathbb{B}(\mathbb{H}(C))\rangle)$.

It is fairly simple to adopt Ehrenfeucht's result to the case with constants for ordinals of the form $\omega^{\omega}(1+\beta)$. For a given ordinal $\beta$, set $D\subseteq \beta$, and ordinal $\alpha<\omega^{\omega}$ let $\gamma=\omega^{\omega}(1+\beta)+\alpha$ and let us look at orders $A$ such that there is a model of $\mathsf{U}_{\gamma}(\omega^{\omega}(1+D))$ which constant-free part is of the form $HF(A)$. It is easy to show that here the possible order types of $A$ are of the form $$\omega^{\omega}+\sum_{\delta\in D}\omega^{\omega}A_{\delta}+\alpha,$$
for some non-empty linear orders $A_{\delta}$. Thus for each $C\subseteq \varepsilon_{\gamma}$ the theory $\mathsf{T}_{\varepsilon_{\gamma}}(C)$ have models with the order types $\varepsilon_A$, where $A$ is of the form $$\omega^{\omega}+\sum_{\delta\in \mathbb{B}(\mathbb{H}(C))}\omega^{\omega}A_{\delta}+\alpha,$$
for some non-empty linear orders $A_{\delta}$. And as in the case of empty $C$, any model $\mathfrak{M}$ of $\mathsf{T}_{\varepsilon_{\gamma}}(C)$ which order type isn't of this form should correspond to a structure with ill-founded sets and hence in $\mathfrak{M}$ there would be an initial interval with the order type $\omega+(\omega^*+\omega)A$ for some dense linear order $A$ without end-points.

[1] A. Ehrenfeucht. An application of games to the completeness problem
for formalized theories. *Fundamenta Mathematicae*, 49:129–141, 1961.