I was just about to look into this problem, and suggest using `plussers()`. Glad you came across it yourself! Unfortunately right now ITensor doesn't have an easy built-in way to do this operation, but Miles and I were discussing functions we could add to help with this.
My quick answer is that you should be able to determine the QN structure from the QN structure of the IQTensor you are using to "pad" your MPS tensors (I don't think the procedure works if you only pad all of your MPS tensors with zeros, otherwise effectively the bond dimension is not being increased).
For example, using this paper as a reference:
https://arxiv.org/pdf/1501.05504.pdf
you can see in Equations (14) and (15) that they pad one tensor with zeros, but the other tensor with a tensor `P`, which in that paper they determine using Equation (19).
Then, the expanded index has a QN structure that is the direct sum of the old index and the index from `P`.
I will look into this to see if there are functions that I can add to ITensor to make this procedure easier. In theory, we could add a function like `directSum(Index i, Index j) -> Index k` where `k` is the direct sum of indices `i` and `j`. We can also add a function like `directSum(ITensor A, ITensor B) -> ITensor C`, where ITensor `C` has indices that are common between ITensors `A` and `B`, and direct sums are performed over the indices that are not common between `A` and `B`. I will try to write a prototype within the next few days that I can send to you, but please let me know if the suggestion above helps you with your own implementation.
Cheers,
Matt