(1) The inputs to the function are:
H = [H_1, H_2, ..., H_N]
L
R
psi = [M_1, M_2, ..., M_N]
On the left side, the boundary tensors M_1 and H_1 are assumed to have dangling indices that are shared with tensor L (tensor L is a 3-index tensor, that also has index dag(l') where l is the Index shared between L and M_1). On the right side, those relations must hold for tensors M_N, H_N, and R.
(2) You can think of the boundary tensors as representing some sort of boundary environment, either finite or infinite (but in general, not periodic). If they were trivial (i.e. L and R were product states of decouple degrees of freedom between the MPS links and MPO links), then it would correspond to an open boundary condition. If the environment tensors L and R came from fixed points of some infinite MPS/MPO, it would correspond to infinite boundary conditions. Otherwise, the boundary conditions just represent some finite bath on either side of the "box" you are optimizing.
One could in theory have periodic boundary conditions, but I believe that would require a single boundary tensor connecting the degrees of freedom between tensors M_1, H_1, M_N, and H_N (i.e. a single order-6 boundary tensor, instead of two decoupled boundary tensors L and R). For periodic boundary conditions with a large enough unit cell (and for typical wavefunctions), these two would be approximately equivalent.
One example I know of is performing time evolution of an MPS within a small window using boundary tensors to represent an infinite bath (
https://arxiv.org/abs/1207.0652 ), which would be related to how this might be used in DMRG. Beyond that, I am not sure where this technique has been used, perhaps Miles is aware of other applications.
Cheers,
Matt