How to couple PDEs of different dimensionality

I am trying to solve a problem in kinetic plasma physics (namely the Vlasov-Poisson-system). The problem is roughly (omitting constants and restricting to 1d1v) given by

\frac{\partial f_e(x,v_x,t)}{\partial t} + v_x \frac{\partial f_e(x,v_x,t)}{\partial x} - E_x(x,t) \frac{\partial f_e(x,v_x,t)}{\partial v_x} = 0
\frac{\partial f_i(x,v_x,t)}{\partial t} + v_x \frac{\partial f_i(x,v_x,t)}{\partial x} + E_x(x,t) \frac{\partial f_i(x,v_x,t)}{\partial v_x} = 0
\frac{\partial E_X(x,t)}{\partial x} = n_i(x,t)-n_e(x,t) = \int_{-\infty}^{\infty}f_i(x,v_x,t) - f_e(x,v_x,t) \mathrm{d}v_x

The timestepping I actually want to solve myself (probably just backward Euler for now), but I need help with how to deal with functions that are defined on different dimensional spaces. As you can see f_e(x,v_x) and f_i(x,v_x) are defined on a 2d phasespace spanning x-v_x. But the electric field E_x(x) only depends on x. To obtain it I need to integrate the distribution functions over v_x and then have to have E_x(x) act on all f(x,v_x) at that x irrespective of v_x. Is this something that fenics can help with?

Searching for keywords that seems relevant to this topic I found mention of the “mixed-dimensional” branch, but I found no example there that is close enough for me to make sense of. I am admittedly new to FEM, so I may be missing something.

PS: I managed to integrate f_e and f_i to obtain n_e and n_i and can even compute E_x(x) as a 1d function, but I have had no luck getting that back into the problem on the 1d1v space.