I don't have a reference, but this is easy to show directly: we have
$$ \| ST \|_{HS}^2 = \sum_{n = 1}^{\infty} \| STe_n \|^{2} \leq \| S\|^2 \sum_{n = 1}^{\infty} \| Te_{n} \|^{2} = \| S \|^2 \| T \|_{HS}^2,$$ where $\| \cdot \|_{HS}$ is the Hilbert-Schmidt norm. This, by itself, however is not enough to prove that the set of Hilbert-Schmidt operators form an ideal - assuming you mean two-sided ideal, you also need to show, with $S$ and $T$ as above, that $TS$ is Hilbert-Schmidt, and that the sum of two Hilbert-Schmidt operators is again Hilbert-Schmidt (along with the perfectly obvious statement that the zero-operator is Hilbert-Schmidt). The latter is easy,
$$ \| S + S' \|_{HS}^2 = \sum_{n = 1}^{\infty} \| (S + S')e_n \|^2 \leq
\sum_{n = 1}^{\infty} (\| Se_n \|^2 + \| S'e_n \|^2) = \|S\|_{HS}^2 + \| S'\|_{HS}^2 < \infty $$
if $S$ and $S'$ are Hilbert-Schmidt, but I expect the former to be more difficult.
Edit: The result follows from the fact that the set of Hilbert-Schmidt operators is $^*$-closed, i.e. if $T$ is H-S, then so is $T^*$. The proof is easy: if $S \in B(H)$ and $T$ is H-S, then $(TS)^* = S^* T^*$ is H-S by the above, and hence $TS$ is H-S also.
Another edit: I'll add that the fact that the class of H-S operators is $^*$-closed follows from the fact that
$$ \| S \|_{HS} = \| S^* \|_{HS}. $$
This is a (small) part of exercise 17c) of section XVIII.9 in Real and Functional Analysis by Serge Lang, so it shouldn't be too difficult to show (however, at the moment, I have no idea how to do it; maybe a sign that I need to revisit some of this stuff).
Let me give a somewhat analogous example:
Definition A sequence of real numbers $(x_n)$ is called convergent if there exists a real number $x$ such that for every $\varepsilon>0$, there exists $n_0$ such that $|x_n-x|<\varepsilon$ for all $n>n_0$. The number $x$ is called a limit of the sequence $(x_n)$.
A simple result is that limits are unique, so we can call $x$ the limit of $(x_n)$, but this is not necessary to the definition of convergence. Indeed, there are spaces where limits need not be unique. If all we are interested in is discussing the convergence or non-convergence of a sequence, the particular limit is unimportant - all we care about is that some limit exists.
Now, how does that relate here?
Definition An operator $T\in\mathcal B(H,F)$ is called a Hilbert-Schmidt operator if there exists an orthonormal basis $(e_n)$ of $H$ such that $\sum_n\|Te_n\|^2<\infty$.
We could continue that definition with something along the lines of "The basis $(e_n)$ is called a Schmidt basis of $T$." (Note: I have made up the term Schmidt basis.) However, rather than uniqueness of Schmidt bases, we have the sort-of opposite result: if $T$ is a Hilbert-Schmidt operator, then every
basis $(f_n)$ of $H$ is a Schmidt basis of $T$. (Hence the term Schmidt basis becomes redundant.)
There is no issue with well-definedness; an operator $T$ is a Hilbert-Schmidt operator if it satisfies the given property for some orthonormal basis, but it is a straightforward result that if $T$ satisfies the property for some basis, then it satisfies it for every basis.
EDIT: Let us show that $\sum_{i\in I}\|Te_i\|^2=\sum_{i\in I}\|Tf_i\|^2$ for any two bases $(e_i)$, $(f_i)$. First note that we may write the left hand side as $\sum_i\langle T^*Te_i,e_i\rangle=:\operatorname{Tr}(T^*T)$. Note that if $A,B\in\mathcal B(H)$ are operators such that both $AB$ and $BA$ are self-adjoint and positive definite, then
$$\operatorname{Tr}(AB)=\sum_i\left\langle A\left(\sum_j\langle Be_i,e_j\rangle e_j\right),e_i\right\rangle=\sum_{i,j}\langle Ae_j,e_i\rangle\langle Be_i,e_j\rangle=\operatorname{Tr}(BA)$$
by symmetry. If $U\in\mathcal B(H)$ is unitary, then both $T^*T$ and $U^*T^*TU$ are self-adjoint and positive definite, and thus
$$\operatorname{Tr}(U^*T^*TU)=\operatorname{Tr}(T^*TUU^*)=\operatorname{Tr}(T^*T).$$
In particular, if $U$ is the unitary operator such that $f_i=Ue_i$ for each $i$, this implies
$$\sum_i\|Tf_i\|^2=\sum_i\langle TUe_i,TUe_i\rangle=\operatorname{Tr}(U^*T^*TU)=\operatorname{Tr}(T^*T)$$
as claimed.
Best Answer
Hilbert-Schmidt operators are more than bounded, as in the case of separable Hilbert space they are compact. Indeed, if $\{e_n\}$ is an orthonormal basis, and $P_N$ the projection over $\operatorname{Span}\{e_j,1\leqslant j\leqslant N\}$, then $\{TP_N\}$ converges in norm to $T$ and $TP_N$ is finite ranked.
The Hilbert-Schmidt operators form an ideal of the set of bounded operators.
An interest of the Hilbert-Schmidt operators is that it can be endowed with an inner product, defining $$\langle S,T\rangle_{HS}:=\sum_{j=1}^{+\infty}\langle Se_n,Te_n\rangle.$$
It can be shown with Bessel's equality that this doesn't depend on the choice of the Hilbert basis.
For applications, I don't know too much. In the context of probability measures on the Borel subsets of the topology induced by an inner product space (with a countable Hilbert basis), Hilbert-Schmidt operators help us to characterize the sets of probability measures which are tight. In Araujo and Giné's book Probability measures in Banach spaces, we encounter the following result: