I'm implementing several "fundamental" operations defined in arXiv:1405.7786 for matrix product states, such as the Hadamard product. Part of doing this is defining operations on the core tensors, such as the partial Kronecker product (see Definition 2.5 of the aforementioned paper, reproduced below). Define two tensors @@\mathbf{A}@@ and @@\mathbf{B}@@, where

$$

A \in \mathbb{R}^{R_1 \times R_2 \times \cdots \times R_M \times I_1 \times \cdots I_N}

$$

and

$$

B \in \mathbb{R}^{S_1 \times S_2 \times \cdots \times S_M \times I_1 \times \cdots I_N}

$$

with @@R_m@@ an auxiliary mode (type Link) and @@I_n@@ a physical mode (type Site). Then the partial Kronecker product is defined as

$$

\mathbf{C} = \mathbf{A} \boxtimes \mathbf{B} \in \mathbb{R}^{R_1 S_1 \times R_2 S_2 \times \cdots \times R_M S_M \times I_1 \times \cdots \times I_N}

$$

where the subtensors are given by

$$

\mathbf{C}(:, \dots, :, i_1, \dots, i_N) = \mathbf{A}(:, \dots, :, i_1, \dots, i_N) \otimes \mathbf{B}(:, \dots, :, i_1, \dots, i_N)

$$

If possible, I'd like to implement this in a general way, where it could be applied to any ITensor. I'm just getting hung up on how to do the indexing correctly and get the outer product that I need. Any help is greatly appreciated!