
\section{Feynman rules for the type I seesaw-extended MSSM}\label{app:FR}
The relevant vertices for the lepton flavor changing self energies are triple vertices for the lepton-slepton-gaugino and -higgsino interactions:
\begin{subequations}\label{eq:FRnuMSSM}
\begin{align}
  i\,\Gamma_{\nu_f}^{\tilde\nu_s \tilde\chi^0_k} = \, &
  -\frac{i}{\sqrt{2}} \bigg\lbrace \left[ (g_2 Z_{2k}^\upsh{N} - g_1
      Z_{1k}^\upsh{N}) \mathcal{Z}^{\tilde\nul *}_{is}
      (\Mat{U}_\upsh{PMNS})_{if} \right]
    P_\Ll \\ & \qquad\qquad + \left[(g_2 Z_{2k}^{\upsh{N}*} - g_1
      Z_{1k}^{\upsh{N}*}) \mathcal{Z}^{\tilde\nul}_{is}
      (\Mat{U}^*_\upsh{PMNS})_{if} \right] P_\Rr\bigg\rbrace ,\nonumber \\
  i\,\Gamma_{e_i}^{\tilde e_s \tilde\chi^0_k} =\, & \frac{i}{\sqrt{2}}
  \bigg\lbrace\left((g_2 Z_{2k}^\upsh{N} + g_1 Z_{1k}^\upsh{N})
      W^{\tilde \el}_{is} - y_{ij}^\el Z_{3k}^\upsh{N} W^{\tilde
        \el}_{j+3,s} \right) P_\Ll \\
    & \qquad\qquad- \left(\sqrt{2} g_1 Z_{1k}^\upsh{N} W^{\tilde
        \el}_{i+3,s} + y_{ji}^{\el*} Z_{3k}^\upsh{N} W^{\tilde
        \el}_{js}\right) P_\Rr \bigg\rbrace,\nonumber \\
  i\,\Gamma_{\nu_f}^{\tilde e_s \tilde\chi^+_k} =\, & -i \left[g_2
    Z_{1k}^- W^{\tilde \el*}_{is} - y_{ij}^{\el*} Z_{2k}^- W^{\tilde
      \el*}_{j+3,s} \right] (\Mat{U}_\upsh{PMNS})_{if} P_\Ll
  ,\\
  i\,\Gamma_{e_i}^{\tilde \nu_s \left(\tilde\chi^+_k\right)^c} =\, &
  -i \left[g_2 Z_{1k}^+ \mathcal{Z}^{\tilde \nul *}_{i,s} P_\Ll -
    y_{ij}^{\el} Z_{2k}^{-*} \mathcal{Z}^{\tilde\nul}_{js} P_\Rr
  \right],
\end{align}
\end{subequations}
where summation over double indices is understood.

The vertices of eqs.~\eqref{eq:FRnuMSSM} are given for an incoming
standard model fermion, outgoing chargino or neutralino as well as
sfermion. They generically follow from a interaction Lagrangian like
\[
\mathcal{L}_\mathrm{int} = \bar f_i \, \Gamma^{\tilde f_s
  \tilde\chi_k}_{f_f} \; \tilde\chi_k \, \tilde f_f \, + \mathrm{\ h.\
  c.}
\]

Each vertex comes along with the corresponding chirality projector:
\[
\Gamma^{\tilde f_s \tilde\chi_k}_{f_f} = \Gamma^{\tilde f_s
  \tilde\chi_k}_{\Ll, f_f} P_\Ll + \Gamma^{\tilde f_s
  \tilde\chi_k}_{\Rr, f_f} P_\Rr
\]

The mixing matrices diagonalize the mass matrices in the following manner:
\begin{itemize}
\item Chargino mixing: \(\Mat{Z}^{-\tp} \mathcal{M}_\upsh{C} \Mat{Z}^+ =
  \left(\mathcal{M}_\upsh{C}\right)^\upsh{diag}\),
 \item Neutralino mixing: \(\Mat{Z}^{\upsh{N}\tp} \mathcal{M}_\upsh{N}
   \Mat{Z}^\upsh{N} = \left(\mathcal{M}_\upsh{N}\right)^\upsh{diag}\),
 \item Slepton mixing: \(\Mat{W}^{\tilde \el \dag} \mathcal{M}^2_{\tilde \el}
   \Mat{W}^{\tilde \el} = \left(\mathcal{M}^2_{\tilde
       \el}\right)^\upsh{diag}\),
 \item Sneutrino mixing: \[\mathcal{W}^{\tilde\nul \dag}
   \bar{\mathcal{M}}_{\tilde\nul}^2 \mathcal{W}^{\tilde\nul} =
   \mathcal{W}^{\tilde\nul \dag} \mathcal{P}^\dag
   \mathcal{M}_{\tilde\nul}^2 \mathcal{P} \mathcal{W}^{\tilde\nul} =
   \left(\mathcal{M}^2_{\tilde\nul}\right)^\upsh{diag},\] such that
   \(\mathcal{Z}^{\tilde\nul} = \mathcal{P}\mathcal{W}^{\tilde\nul}\)
   diagonalizes the original mass matrix \(\mathcal{M}^2_{\tilde\nul}\)
   and therefore:
   \begin{align*}
     \mathcal{Z}_{is}^{\tilde \nul} &= \frac{1}{\sqrt{2}}
   \left( \mathcal{W}_{is}^{\tilde \nul} + i \mathcal{W}_{i+3,s}^{\tilde
       \nul}\right) \quad\text{and} \\
   \mathcal{Z}_{i+3,s}^{\tilde \nul} &=
   \frac{1}{\sqrt{2}} \left(\mathcal{W}_{is}^{\tilde \nul} - i
     \mathcal{W}_{i+3,s}^{\tilde \nul}\right)
   \end{align*}
   appear in the vertices of eqs.~\eqref{eq:FRnuMSSM}.
 \item Neutrino mixing: The PMNS mixing matrix can be determined from
   the neutrino mass matrix \(\Mat{m}_\nul^\upsh{diag} = \Mat{U}_\upsh{PMNS}^*
   \Mat{m}_\nul \Mat{U}_\upsh{PMNS}^\dag\), where \(m_\nul\) is the effective light
   neutrino mass matrix and the charged lepton masses can be taken
   diagonal (otherwise there would be a contribution to the PMNS mixing
   similar to the CKM mixing from both up and down sector:
   \(\Mat{U}_\upsh{PMNS} = \Mat{V}_{\el,\Ll}^\dag \Mat{U}_{\nul,\Ll}\), where
   \(\Mat{V}_{\el,\Ll}\) rotates the left-handed electron fields).
\end{itemize}

\paragraph*{Loop functions}
Finally, we give our conventions for the arising \(\MSbar\) subtracted
loop functions for \(p^2=0\) and \(\mu_r\) the renormalization scale:
\begin{subequations}
\begin{align}
  B_0(m_1, m_2) &= - \ln\frac{m_1 m_2}{\mu_r^2} + 1 -\frac{m_1^2 + m_2^2}{m_1^2 - m_2^2} \ln\frac{m_1}{m_2}, \\
        B_1(m_1, m_2) &= \frac{1}{2} \ln \frac{m_1 m_2}{\mu_r^2} - \frac{3}{4} - \frac{m_2^2}{2(m_1^2 - m_2^2)} + \left(\frac{m_1^4}{(m_1^2 - m_2^2)^2} - \frac{1}{2}\right) \frac{m_1}{m_2}.
\end{align}
\end{subequations}
