Solution (source code)

= Solution

Put $d=|\mathbf r_1-\mathbf r_2|$, $g=e^{ikd}/d$ and $s_j=e^{i\mathbf k_i\cdot\mathbf r_j}$. In the range-to-wavelength overlap region around centre $j$, the outgoing wave's singular and constant terms are $c_j/(k\rho_j)+ic_j$. Matching their ratio to the <Bethe-Peierls boundary condition> yields
$$
(1+ika_1)c_1+a_1g c_2=-ka_1s_1,\qquad
(1+ika_2)c_2+a_2g c_1=-ka_2s_2.
$$
Define $A_j=1+ika_j$ and $D=A_1A_2-a_1a_2g^2$. Solving the coupled system gives
$$
\boxed{c_1=-ka_1\frac{A_2s_1-a_2gs_2}{D},\qquad
c_2=-ka_2\frac{A_1s_2-a_1gs_1}{D}}.
$$
This is the <two-centre scattering-length model> result, retaining the outgoing-wave $ik$ term and all repeated exchange between the two centres. If terms of order $ka_j$ are neglected in the local matching, set $A_1=A_2=1$ in these formulae; if also $kd\ll1$, set $g\simeq1/d$ and $s_j\simeq1$. Away from a vanishing denominator this gives $c_1\simeq-ka_1(1-a_2/d)/(1-a_1a_2/d^2)$ and its exchanged partner. The isolated-centre limit is $c_j=-ka_js_j/(1+ika_j)$. Finite-range corrections require additional effective-range data, not supplied here.