The coupling is brought about by the overlap of the relevant wave functions. To appreciate this, you should consider the expression for the polarization function, or the dielectric function, which (assuming zero temperature for convenience) involves matrix elements corresponding to the N-particle ground state and the N-particle excited states of the system under investigation. In a perturbational approach (notably in the random-phase approximation, RPA), where ultimately the relevant N-particle states are expressed in terms of N-particle Slater determinants composed of single-particle wave functions, one can immediately notice the nature of the coupling mediated by the overlap of the underlying single-particle wave functions.
Surface plasmons are electron density oscillation coupled with electromagnetic field. The interaction of electric charge on the outer and inner surfaces of the shells are due to the electric field. This coupling strongly depends on the geometry. One has to solve the electrostatic problem. For spherical shells, the problem is not difficult.
We are not pioneers in this field but you may look at DOI:10.1016/j.chemphys.2012.03.011
There are two coupled plasmon modes for each pair of inner/outer shell mode. One can see a mode-splitting of the mode energy depending on the shell thickness, which is captured phenomenologically by the hybridization model [10.1126/science.1089171]. Some models for calculation are available [10.1021/nl101335z]
Analytically these mode can be captured by Mie-calculations for spherical shells, but for more complex geometries only numerical methods help.