Neutrino oscillation experiments have unambiguously proved the fact that neutrinos have tiny non-zero mass and they mix among each other. The discovery of neutrino oscillations reveals that charged lepton flavour violating phenomena must occur. In the Standard Model, due to massless neutrinos lepton flavour violating decays like $\mu \to e \gamma$ are forbidden. Even if we go one step further and assume that there exist massive neutrinos in the Standard Model, the rates for lepton flavour violating processes are absurdly small (for example ${\rm BR}(\mu \to e \gamma) \lesssim 10^{-54}$) owing to the GIM mechanism and the fact that neutrino masses are tiny compared to the weak scale. This is clearly undetectable from the experimental point of view. Thus any observable signal to this process will indisputably prove the existence of new physics. In many new physics models, new particles occur at the loop-level and contribute considerably to these lepton flavour violating processes.

The LFV module in $\texttt{HEPfit}$ calculates the decay rates for various lepton flavour violating observables in a given new physics model. The expression for various $\Delta F = 1$ and $\Delta F = 0$ processes are calculated in terms of the Wilson coefficients of the responsible operators. Therefore, given any new physics model, once the user provides the Wilson coefficients generated in the model, $\texttt{HEPfit}$ will calculate the decay rates for all the lepton flavour violating observables generated by those operators. This enables the user to include any new physics model to $\texttt{HEPfit}$ and analyze the lepton flavour violating processes efficiently. For example, the effective weak Hamiltonian for the $\ell_i \to \ell_j\, \gamma$ process is given by
\begin{align}
\mathcal{H}_{\rm eff}^{\Delta F=1} &=
\mathcal{C}_{7} \mathcal{O}_{7} + \mathcal{C}_{7}^{\prime} \mathcal{O}_{7}^{\prime} , \nonumber
\end{align}
where
\begin{align}
\mathcal{O}_{7} &= e\, m_{l_i}
\overline{\ell}_j \sigma_{\mu \nu} P_R \ell_i F^{\mu \nu}\ , \notag \\
\mathcal{O}_{7}^{\prime} &= e\, m_{l_i}
\overline{\ell}_j \sigma_{\mu \nu} P_L \ell_i F^{\mu \nu}\ , \nonumber
\end{align}

$\sigma_{\mu \nu} = \frac{i}{2}[\gamma^{\mu},\gamma^{\nu}]$, $P_{R,L} = \frac{1}{2} (1\pm\gamma_5)$, and $F_{\mu\nu} = (\partial_{\mu} A_{\nu} - \partial_{\nu} A_{\mu})$. $m_{l_i}$ is the mass of the initial state fermion. The decay rate for the process is
\begin{align}
\Gamma (\ell_i \to \ell_j\, \gamma) &= \frac{e^2}{4\pi} m_{l_i}^5
\bigg(|\mathcal{C}_{7}|^2 + |\mathcal{C}_{7}^\prime|^2\bigg) \ . \nonumber
\end{align}
In the present version of the HEPfit, all the Wilson coefficients for the processes $\ell_i \to \ell_j\, \gamma$, $\ell_i \to 3\ell$, $\mu \to e$ conversion as well as $(g-2)_{\mu}$ at two-loop level is included in the general MSSM.