Multi mode radiation boundary operators for wave problems

Hello all,
I have got a lot of mileage out of the simple waveguide absorbing boundary condition for electromagnetics
\mathbf{n\times E} = \overline{\mathbf{Z}}\{\mathbf{n\times H}\},
which works well if you have only a single mode on the boundary. This is very easy to implement in Fenics. My question is, can one, in Fenics, construct boundary operators of the type
b_{ij} = \sum\limits_m \int\limits_{\partial \Gamma}\int\limits_{\partial\Gamma}\langle v_i\vert \phi_m\rangle\langle\phi_m\vert u_j\rangle \, dS \, dS^\prime,
where u, v are the usual Fenics trial and test functions, respectively and the sum is taken over all the desired modes? The variable \phi_m represents the m^{th} waveguide mode eigenfunction. At a practical level, we need to construct the outer product of the vectors formed by projecting the mode functions \phi onto the finite element space. Is there a way to do this in Fenics when constructing the forms? This is needed if one wants to consider cases where more than one mode is propagating or where evanescent modes have not died out sufficiently before the exit port is encountered.

I am (still) using Fenics 2019.2.0 (which works very well…I am loathe to break during an upgrade it until I have to). If anyone can offer any hints on where to start, I would be very grateful!

I’m not familiar with the application, and the use of bra–ket notation here is unclear to me. (For instance, if each bra–ket product is an L^2 inner product over the domain, then there would be no dependence on spatial coordinate left for the outer integral \int_{\partial\Gamma}(\cdots)\,dS.)

Thanks for your reply @kamensky. I edited the expression above to indicate more precisely what I would like to do. To put this into context of electromagnetic scattering into multiple modes, I would like to augment the Lagrangian with the following boundary intergral(s).
L_t = L + \sum\limits_m\frac{j k_0\eta_0}{2}{\int\!\!\!\int}_{\partial\Gamma_{port}}\mathbf{E}^*\times \mathbf{h}_m {\int\!\!\!\int}_{\partial\Gamma_{port}} \mathbf{E}\times \mathbf{h}_m \, dS \, dS^\prime.
L is the usual Lagrangian form for the Maxwell system. L_t is the Lagrangian augmented with the boundary operator. The electric field \mathbf{E} is represented by our usual test and trial functions. The (known a priori) normalized waveguide transverse magnetic field eigenfunctions \mathbf{h}_m can also be expanded in terms of the trial functions. In short, we need to form a vector of the inner products of the trial functions and the mode eigen function and then form the outer product of this vector with itself. This produces a boundary operator matrix that is densely populated over the affected degrees of freedom. Can something like this be constructed for incorporation in the Fenics variational form?

Perhaps you could introduce auxiliary global unknowns to represent surface integrals, by using a mixed element combining the element for \mathbf{E} with several additional elements of the type VectorElement("Real",cell,0). (See, e.g., this demo for an example of a scalar field coupled to a global scalar unknown.) Then, if I understand the updated equations correctly, you could add terms to the variation of the Lagrangian (\delta L) of the form

jk_0 \eta_0 \left(\iint_{\partial\Gamma} \left(\frac{1}{\vert\partial\Gamma\vert}\mathbf{I}_m - \mathbf{E}\times\mathbf{h}_m\right)\cdot\mathbf{J}_m\,dS + \iint_{\partial\Gamma}\left(\delta\mathbf{E}\times\mathbf{h}_m\right)\cdot\mathbf{I}_m\,dS\right)\text{ ,}

where \mathbf{I}_m is a global vector-valued unknown in \mathbb{R}^3, \mathbf{J}_m\in\mathbb{R}^3 is a corresponding test function, \delta\mathbf{E} is a test function in the electric field’s test space, \vert\partial\Gamma\vert is the measure of \partial\Gamma, i.e., \iint_{\partial\Gamma}1\,dS, and I just gave up and pretended everything was real-valued, to simplify all the complex conjugate bookkeeping, even though it looks like you’re probably using complex numbers. (Is \mathbf{E}^* meant to be the complex conjugate of \mathbf{E}? Should one of the instances of \mathbf{h}_m in your formula also be conjugated?)

The first integral implies that, if the formulation is satisfied for all \mathbf{J}_m\in\mathbb{R}^3, then

\mathbf{I}_m = \iint_{\partial\Gamma}\mathbf{E}\times\mathbf{h}_m\,dS\text{ .}

If this is true, then the second integral adds the desired term to \delta L. Hopefully at least the main idea is clear, and it can be extended to the complex case without too much effort.