\section{Introduction}
The recent exciting first detections of gravitational waves~\cite{ligo,neutronstar}, which marked a new era in astrophysics and cosmology, have pushed the scientific community towards the construction of ever more sophisticated ground and space based detectors~\cite{tian,decigo,et,kagra,lisa} to observe waves in a variety of ranges in order to sketch a first map of the stochastic gravitational background. Within the framework of quantum theory, a stochastic gravitational background affects the dynamics of matter propagation~\cite{aha,linet} and, when the quantum state is in a superposition, it leads to decoherence effects, as typical of noisy environments. Different models for the description of this phenomenon have been proposed~\cite{goku,breuer,sanchez,power,blencowe,ana,lamine}. However, they do not agree on the decoherence mechanism (the preferred basis and  rates) at which it takes place. In order to solve such apparent contradictory results, we have derived in \cite{bosons} a novel model for the decoherence effect induced by a stochastic gravitational perturbation on non relativistic scalar bosonic particles. Our model has so far proven to be able to describe more general scenarios than those present in the literature, as it is able to qualitative recover them as appropriate limiting cases, thus solving the decoherence basis puzzle. However, it might not be general enough to describe the outcome of a real experiment. The particles commonly employed in experiments (atoms, neutrons, electrons...) in fact have a charge, a spin and could be coupled to other external fields, like the Maxwell one for instance. For the above reasons, in this paper we will derive an analogous model, this time for spin 1/2 fermions interacting with both a gravitational perturbation and an external electromagnetic field.\\
\\
The paper is organized as follows. In section \ref{ii} we derive the equations of motion in Hamiltonian form for a spin 1/2 fermionic field minimally coupled to a weakly perturbed flat metric. We then specialize such equation to the non relativistic regime in section \ref{iii} and proceed with the canonical quantization of the bosonic field in the single particle sector, obtaining a Schr\"{o}dinger like equation for a test particle interacting with a weakly perturbed gravitational field.\\
In section \ref{iv} we compare the here derived fermionic model with the bosonic one derived in \cite{bosons}.\\
In section \ref{v} we specialize to the case of a stochastic gravitational perturbation and derive the corresponding master equation. We discuss the decoherence effect with explicit reference to the preferred eigenbasis and characteristic decoherence time. We also show under which assumptions our master equation is able to reproduce decoehrence in the position or momenutm eigenbasis only thus recovering the results of the literature \cite{bosons,ana,blencowe,breuer,sanchez,power}.
\section{Equations of motion}\label{ii}
We first derive the equations of motion (EOM) for a spin 1/2 fermionic field minimally coupled to linearized gravity. We start from the action for the Dirac field in curved spacetime~\cite{bd}:
\begin{equation}\label{actionf}
\begin{split}
S = \int d^4x \sqrt{-g}\mathcal{L}
\end{split}
\end{equation}
with the Lagrangian density:
\begin{equation}
\mathcal{L}  = \frac{i\hbar c}{2}[ \bar{\psi}\gamma^{\mu}{e^{\small{A}}}_\mu\mathcal{D}_{\small{A}}\psi -{e^{\small{A}}}_\mu\mathcal{D}_{\small{A}}\bar{\psi}\gamma^{\mu}\psi  ] - mc^2 \bar{\psi}\psi
\end{equation}
where ${e^A}_\mu(x)$ is the so called vierbein field \cite{vierbein}, an auxiliary field used in order to extend the definition of fermions as irreducible spin 1/2 representations of the Poincar\'{e} group to curved spacetimes (see appendix A), and
\begin{equation}
\mathcal{D}_{\mu} \psi= \partial_{\mu}\psi +\frac{1}{8}[\gamma_a,\gamma_b]{\omega_\mu}^{ab} \psi + \frac{ie}{\hbar c}A_{\mu}\psi
\end{equation}
is the covariant derivative with respect to both the spin ($\omega_{\mu}^{ab}$) and the electromagnetic ($A_\mu$) connections.
In this framework the metric tensor $g_{AB}$ and the affine connection $\Gamma_{A}^{BC}$ are describeb by the following pair of equations:
\begin{eqnarray}\label{veir}
\left\lbrace
\begin{array}{cc}
{e_{A}}^\mu\eta_{\mu\nu}{e_{B}}^{\nu} =  g_{AB}\\
{\omega_{A}}^{\mu\nu} = {e_B}^{\mu}\eta^{\nu\rho}\partial_A {e^{B}}_\rho
+ {e_{B}}^{\mu}\eta^{\nu\rho} {e^{C}}_\rho{\Gamma^{B}}_{AC}
\end{array}
\right.
\end{eqnarray}
Note that Eq.~\eqref{veir} holds only for a torsion free, metric compatible connection \cite{vierbein}. \\
In order to describe a weak perturbation of the metric, we now write the metric as the sum of a flat background $\eta_{\mu\nu} = \text{diag}(+---)$, and a perturbation $h_{\mu\nu}$:
\begin{equation}
g_{\mu\nu} = \eta_{\mu\nu}+h_{\mu\nu},
\end{equation}
and since we are interested in studying the dynamics of the Dirac field interacting with a weak gravitational perturbation, perform a Taylor expansion of the fermioni action around the flat background metric and truncate the series at the first perturbative order (See Appendix B for the explicit calculation). Thus, we obtain the effective Lagrangian $\mathcal{L}_{eff}$ acting on flat spacetime:
\begin{equation}\label{effectiveaction}
\begin{split}
S &= \int d^4 x \:\Big(\frac{i\hbar c}{2}[\bar{\psi}\gamma^{\mu}\nabla_\mu\psi - \nabla_\mu(\bar{\psi})\gamma^{\mu}\psi](1+\frac{tr(h)}{2})\\
-&(1+\frac{tr(h)}{2})mc^2\bar{\psi}\psi- \frac{i\hbar c}{4}h_{\mu\nu}[\bar{\psi}\gamma^{\mu}\nabla^{\nu}\psi - \nabla^\nu(\bar{\psi})\gamma^{\mu}\psi ]\Big)\\
+& O(h^2)\\
&\equiv \int d^4 x\: \mathcal{L}_{eff} + O(h^2)
\end{split}
\end{equation}
where $\nabla_{\alpha}$ is the flat covariant derivative with respect to the electromagnetic connection. The EOM for the matter field are obtained (at first order in the perturbation $h_{\mu\nu}$) from the Euler Lagrange equations:
\begin{equation}
\frac{\partial \mathcal{L}_{eff}}{\partial \bar{\psi}} - \nabla_\alpha \frac{\partial \mathcal{L}_{eff}}{\partial \nabla_\alpha \bar{\psi} } = 0
\end{equation}
and in the harmonic gauge they read:
\begin{equation}\label{dirac}
\begin{split}
&i\hbar\partial_t\psi = eA_0\psi +  mc^2(1+\frac{h_{00}}{2})\gamma^0\psi - \frac{mc^2}{2}h_{0j}\gamma^j\psi\\
&-i\hbar c (1+\frac{h_{00}}{2})\gamma^0\gamma^i(\partial_i+\frac{ie}{\hbar c}A_i)\psi + \frac{i\hbar c}{2}h_{0i}(\partial^i+\frac{ie}{\hbar c}A^i)\psi\\
&+ \frac{i\hbar c}{2}h_{ij}\gamma^0\gamma^i(\partial^j+\frac{ie}{\hbar c}A^j)\psi+ \frac{i\hbar c}{2}h_{0i}\gamma^i\gamma^j(\partial_j+\frac{ie}{\hbar c}A_j)\psi\\
&- \frac{i\hbar c}{8}\partial_\alpha(tr(h))\gamma^0\gamma^\alpha\psi + O(h^2)\psi\\
&\quad\quad\: = : \mathcal{H}\psi +O(h^2)\psi
\end{split}
\end{equation}
As in the case of a scalar field  studied in~\cite{bosons} , we are not allowed to give a probabilistic interpretation to the field $\psi$, becuse the conserved charged $Q$ associated to the internal $U(1)$ symmetry ($\psi \rightarrow e^{ie}\psi$ ; $\bar{\psi}\rightarrow e^{-ie}\bar{\psi}$) via Noether's Theorem reads:
\begin{equation}\label{chargef}
\begin{split}
Q &\equiv -ie \int d^3 x \: \Big(\frac{\partial\mathcal{L}_{eff}}{\partial(\nabla_0 \psi)}\psi - \psi^\dagger \frac{\partial\mathcal{L}_{eff}}{\partial(\nabla_0 \psi^\dagger)}\Big)\\
&= \hbar e c \int d^3 x\: \Big( \psi^\dagger (1-tr(h)-\frac{h_{00}}{2})\psi - \psi^\dagger \frac{h_{0i}}{2}\gamma^0\gamma^i\psi  \Big)
\end{split}
\end{equation}
instead of the required:
\begin{equation}\label{eq:stchf}
\rho = \int d^3 x \:
\psi^\dagger \psi
\end{equation}
We therefore apply the transformation:
\begin{eqnarray}\label{eq:chargpres}
\left\lbrace
\begin{array}{ccc}
T &=& (1 -\frac{tr(h)}{2}-\frac{h_{00}}{4} -\frac{h_{0i}}{4}\gamma^0\gamma^i )\\
 \psi &\rightarrow& T \psi \\
\mathcal{H}&\rightarrow& \mathfrak{H} := T \mathcal{H}T^{-1}+ i\hbar T \partial_t(T^{-1})
\end{array}
\right.
\end{eqnarray}
so that, in the new representation,the conserved charge can be expressed by the standard form in Eq.~\eqref{eq:stchf}. \\
After some algebra the EOM~\eqref{dirac} read:
\begin{equation}\label{eqmotf}
i\hbar\partial_t \psi = [mc^2\gamma^0 + \mathfrak{E}+ \mathcal{O}]\psi
\end{equation}
where
\begin{equation}
\begin{split}
\mathfrak{E} = & e A_0 +\frac{mc^2}{2}h_{00}\gamma^0 + i\hbar c \: h_{0i} (\partial^i-\frac{ie}{\hbar c}A^i)+\frac{i\hbar c}{4}\partial_{i}(h_0^i)\\
&+\frac{\hbar c}{4}\epsilon^{ijk}\partial_i(h_{0j})\Sigma_k-\frac{3i\hbar}{8}\partial_t(tr(h))+\frac{i\hbar}{4}\partial_t(h_{00})
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\mathcal{O} =& - i\hbar c (1+\frac{h_{00}}{2}) (\partial_j - \frac{ie}{\hbar c}A_j)\alpha^j +\frac{i\hbar}{4}\partial_t(h_{0i})\alpha^i \\
&+ \frac{i\hbar c}{2}h_{ij}(\partial^j - \frac{ie}{\hbar c}A^j)\alpha^i + \frac{i\hbar c}{4}\partial_i(\frac{tr(h)}{2}-h_{00})\alpha^i
\end{split}
\end{equation}
are respectively the even (diagonal) and odd (off diagonal) parts of the Hamiltonian $\mathfrak{H}$, with $\alpha^{\mu}=\gamma^0\gamma^\mu$ and $\Sigma^i = \textrm{diag}(\sigma^i,\sigma^i)$.\\
We are interested in the description of the dynamics of a positive energy particle system in the non relativistic limit. In such a limit, the particle and antiparticle sectors are non interacting with one another, that is to say, the EOM~\eqref{dirac} can be recast to a system of two decoupled equations respectively for the large ($\psi_L$) and the small ($\psi_s$) component of the bispinor $\psi = \left(\begin{array}{cc} \psi_L \\ \psi_s \end{array}\right)$. While this is evident and straightforward for the free case \cite{greiner}, for an interacting theory decoupling the two components is a very complicated task that can only be achieved perturbatively.\\
\\
In the next section we will provide a standard prescription for the diagonalization of the EOM in the non relativistic limit.
\section{Non relativistic limit and canonical quantization}\label{iii}
We aim to find a representation of the bispinor field $\psi$ in which the EOM \eqref{eqmotf} are diagonal. This representation can be found in non relativistic limit
following the Foldy-Wouthuysen Method~\cite{foldy}, which allows one to write perturbatively (at any order in $\frac{v}{c}$) two decoupled equations, one for each component of the field. The method is operatively characterized by the application of an appropriate unitary transformation $U$:
\begin{equation}
\psi  \rightarrow  \psi^\prime = U\psi
\end{equation}
\begin{equation}
\begin{split}
\mathfrak{H} \rightarrow  \mathfrak{H}^\prime =& \:U(\mathfrak{H}-i\hbar\partial_t)U^{-1}\\
 =& \: mc^2 \gamma^0+ \mathfrak{E}^{\prime} + \mathcal{O}^{\prime} +O(h^2)
\end{split}
\end{equation}
such that, in the new representation, the antidiagonal part $\mathcal{O}^\prime$ is of higher order in $\frac{v}{c}$ than the diagonal $\mathfrak{E}^\prime$. By neglecting $\mathcal{O}^\prime$ one recovers two decoupled equations.
By performing iteratively the transformation, one can always find a representation of the bispinor field for which the EOM are diagonal at any desired order in $\frac{v}{c}$.\\
In our case, the task is easily achieved by applying the subsequent transformations:
\begin{eqnarray}\label{transf}
\left\{
\begin{array}{ccc}
U =& e^{- i \gamma^0 \mathcal{O} /(2mc^2)} \\
U^\prime =& e^{-i\gamma^0 \mathcal{O^\prime} /(2mc^2)}\\
U^{\prime\prime} =& e^{-i\gamma^0 \mathcal{O^{\prime\prime}} /(2mc^2)}
\end{array}
\right.
\end{eqnarray}
after which, with some algebra (see Appendix C) and by neglecting the terms containing the derivatives of the gravitational perturbation of order $\frac{v^3}{c^3}$ or higher\footnote{For the sake of compactness we relegate such terms to Appendix D. Note also that in most experimental situations such contributions are negligible in the case of gravity, and they wouldn't add any further informative content to the analysis in any case.}, the Hamiltonian density to order $\frac{v^4}{c^4}$ reads:
\begin{widetext}
\begin{equation}\label{hamf}
\begin{split}
    H=&  eA_0+\gamma^{0}\Bigg[mc^2(1+\frac{h_{00}}{2})-\frac{\hbar^2}{2m}(1+\frac{h_{00}}{2})(\mathbf{\nabla}-\frac{ie}{\hbar c}\mathbf{A})^2 -\frac{\hbar e}{2m c}(1+\frac{h_{00}}{2})B^k\Sigma_k -\frac{\hbar^2}{2m}h_{ij}(\partial^i-\frac{ie}{\hbar c}A^i)(\partial^j-\frac{ie}{\hbar c}A^j)\\
    &+\frac{\hbar e}{4 mc}\epsilon^{ijl}h_{jk}{F_i}^k\Sigma_l \Bigg] +\frac{i\hbar^2 e}{4m^2c^2}(1+\frac{h_{00}}{2})\Big( \frac{\mathbf{\nabla}}{2}\times\mathbf{E} - \mathbf{E}\times\mathbf{\nabla}\Big)\cdot\boldsymbol{\Sigma} - (1+h_{00})\frac{\hbar^2e}{8m^2c^2}\mathbf{\nabla}\cdot\mathbf{E}\\
    &-\frac{i\hbar^2 e}{16m^2c^2}\epsilon^{ikl}h_{ij}\partial^{j}(E_k)\Sigma_l -\frac{i\hbar^2e}{8m^2c^2}\epsilon^{ikl}h_{ij}E_k(\partial^j-\frac{ie}{\hbar c}A^j)\Sigma_l + \frac{i\hbar^2 e}{4m^2c^2}\epsilon^{ijl}h_{0k}{F_j}^k(\partial_i-\frac{ie}{\hbar c}A_i)\Sigma_l\\
    &-\frac{\hbar^2 e}{8m^2c^2}h_{0j}\partial_i(F^{ij})+\frac{i\hbar^2 e}{8m^2c^2}\epsilon^{ijl}h_{0k}\partial_i({F_j}^k)\Sigma_l-\frac{\gamma^0}{8m^3c^6}\Bigg[\hbar^4c^4(1+2h_{00})(\mathbf{\nabla}-\frac{ie}{\hbar c}\mathbf{A})^4 + \hbar^2 e c^2 (1+2h_{00})B^2 \\
&+ 2 \hbar^4 c^4 h_{ij}(\mathbf{\nabla}-\frac{ie}{\hbar c}\mathbf{A})^2(\partial^i-\frac{ie}{\hbar c }A^i)(\partial^j-\frac{ie}{\hbar c }A^j) - \frac{\hbar^3 e c^3}{2}\epsilon^{ijl}h_{jm}{F_i}^mB^k\lbrace \Sigma_k ,\Sigma_l\rbrace\\
& +\frac{\hbar^3 e c^3}{2}\epsilon^{ijl}\lbrace (\mathbf{\nabla}-\frac{ie}{\hbar c}\mathbf{A})^2 , h_{jk}{F_i}^k \rbrace \Sigma_l-\hbar^3 e c^3 (1+2h_{00})\lbrace(\mathbf{\nabla}-\frac{ie}{\hbar c}\mathbf{A})^2,B^k \rbrace\Sigma_k \Bigg]\\
&+H_d+O(h^2)+O(\partial h) +O(\frac{v^5}{c^5})
\end{split}
\end{equation}
\end{widetext}
where $B$ and $E$ are the magnetic and electric field, and in terms of the four-potential they read:
\begin{eqnarray}
\left\lbrace
\begin{array}{cc}
\mathbf{E} =& - \mathbf{\nabla} A_0 - \frac{1}{c}\dot{\textbf{A}}\\
\mathbf{B} =& \mathbf{\nabla}\times\mathbf{A}\\
B^k =& -\frac{1}{2} \epsilon^{ijk}F_{ij}\\
F_{ij} =& - \epsilon_{ijk}B^{k}
\end{array}
\right.
\end{eqnarray}
$\epsilon_{ijk}$ represent the Levi-Civita symbol, and
\begin{equation}
\begin{split}
H_d= &-\frac{\hbar^2}{8m}\partial_i(h_{00})(\partial^i-\frac{ie}{\hbar c}A^i)\gamma^0 -\frac{\hbar^2}{16m}\partial^i\partial_i(h_{00})\gamma^0 \\
&+\frac{i\hbar c}{4}\partial_{i}(h_0^i)+\frac{\hbar c}{4}\epsilon^{ijk}\partial_i(h_{0j})\Sigma_k -\frac{3i\hbar}{8}\partial_t(tr(h))\\
&+\frac{i\hbar}{4}\partial_t(h_{00}) +\gamma^0 \Bigg[ \frac{\hbar^2}{2m}\partial^i(h_{00})\nabla_i -\frac{\hbar^2}{4m}\partial^i(h_{ij})\nabla^j\\ &-\frac{\hbar^2}{2m}\partial_i(\frac{tr(h)}{2}-h_{00})\nabla^i-\frac{\hbar^2}{4m}\partial^i\partial_i (\frac{tr(h)}{2}-h_{00})\\
&-\frac{i\hbar^2}{4m}\epsilon^{ijk}\Big(\partial_i(h_{00})\nabla_j-\partial_i(h_{jl})\nabla^l \Big)\Sigma_k \Bigg]
\end{split}
\end{equation}
Note that as the transformations \eqref{transf} are unitary~\cite{armin}, they preserve the conserved charge in \eqref{chargef}, i.e. the probability density in the non relativistic limit.\\
In the non relativistic limit the EOM \eqref{hamf} do not mix the two components $\psi_L$ and $\psi_s$  of the field (up to a very small correction). As we are interested in the dynamics of particles only, we restrict the analysis to the first field component $\psi_L$, that we rename as $\psi$ in what follows.\\
Since the dynamics preserves the probability density, we are allowed to apply the canonical quantization prescription and impose the equal time commutation relations:
\begin{equation}
\begin{split}
[\hat{\psi}(t,\mathbf{x}),\hat{\psi}(t,\mathbf{x'})] =  & [\hat{\psi}^{\dagger}(t,\mathbf{x}),\hat{\psi}^{\dagger}(t,\mathbf{x'})] =0\\
[\hat{\psi}(t,\mathbf{x}),\hat{\psi}^{\dagger}(t,\mathbf{x'})] =& \: \delta^{3}(\mathbf{x}-\mathbf{x'})
\end{split}
\end{equation}
to obtain the EOM for the quantum field. The equation thus obtained does not allow for the creation or annihilation of particles. We can thus safely project it onto a single particle sector to obtain the single particle Schr\"{o}dinger like equation:
\begin{equation}\label{schrf}
i\hbar\partial_t\vert\phi(t)\rangle = (\hat{H}_0+\hat{H}_r+\hat{H}_p+\hat{H}_{rp}+\hat{H}_d)\vert\phi(t)\rangle
\end{equation}
with:
\begin{equation}
\hat{H}_{0} = mc^2 +\frac{1}{2m}\Big(\hat{\mathbf{p}}-\frac{e}{c}\mathbf{A}(\hat{\mathbf{x}})\Big)^2+eA_0(\hat{\mathbf{x}})- \frac{\hbar e}{2mc}\mathbf{B}(\hat{\mathbf{x}})\cdot\boldsymbol{\sigma}
\end{equation}
\begin{equation}
\begin{split}
\hat{H}_r =& \frac{\hbar e}{4m^2c^2}\Big( \frac{\hat{\mathbf{p}}}{2}\times\mathbf{E}(\hat{\mathbf{x}}) - \mathbf{E}(\hat{\mathbf{x}})\times\hat{\mathbf{p}}\Big)\cdot\boldsymbol{\sigma} \\
&- \frac{\hbar^2e}{8m^2c^2}\mathbf{\nabla}\cdot\mathbf{E}(\hat{\mathbf{x}})-\frac{\gamma^0}{8m^3c^6}\Bigg[c^4(\hat{\mathbf{p}}-\frac{e}{ c}\mathbf{A}(\hat{\mathbf{x}}))^4 \\
&+ \hbar^2 e c^2 B^2(\hat{\mathbf{x}})-\hbar e c^3 \lbrace(\hat{\mathbf{p}}-\frac{e}{ c}\mathbf{A}(\hat{\mathbf{x}}))^2,B^k(\hat{\mathbf{x}}) \rbrace\sigma_k \Bigg]
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\hat{H}_p =& +\frac{mc^2}{2}h^{00}(t,\hat{\mathbf{x}})- \frac{1}{8m}\lbrace h^{00}(t,\hat{\mathbf{x}})\:,\:\Big(\hat{\mathbf{p}}-\frac{e}{c}\mathbf{A}(\hat{\mathbf{x}})\Big)^2\rbrace \\
&+ \frac{c}{2} \lbrace h_{0i}(t,\hat{\mathbf{x}})\: ,\: \hat{p}^i\rbrace -\frac{\hbar e }{4mc}\epsilon^{ikl}h_{ij}(\hat{\mathbf{x}},t){F^j}_k(\hat{\mathbf{x}})\sigma_l \\
&-\frac{1}{4m}\lbrace h^{ij}(t,\hat{\mathbf{x}})\: ,\:\Big(\hat{p}_i-\frac{e}{c}A_i(t,\hat{\mathbf{x}})\Big)\Big(\hat{p}_j-\frac{e}{c}A_j(t,\hat{\mathbf{x}})\Big)\rbrace\\
&- \frac{\hbar e}{4mc}h_{00}(\hat{\mathbf{x}},t)\mathbf{B}(\hat{\mathbf{x}})\cdot\boldsymbol{\sigma}
\end{split}
\end{equation}
\begin{widetext}
\begin{equation}
 \begin{split}
\hat{H}_{rp}=& \frac{\hbar e}{16m^2c^2}\lbrace h_{00}(\hat{\mathbf{x}},t),\Big( \frac{\hat{\mathbf{p}}}{2}\times\mathbf{E}(\hat{\mathbf{x}}) - \mathbf{E}(\hat{\mathbf{x}})\times\hat{\mathbf{p}}\Big)\cdot\boldsymbol{\sigma}\rbrace -\frac{i\hbar^2 e}{16m^2c^2}\epsilon^{ikl} h_{ij}(\hat{\mathbf{x}})\partial^{j}(E_k(\hat{\mathbf{x}}))\Sigma_l - \frac{\hbar^2e}{8m^2c^2}h_{00}(\hat{\mathbf{x}})(\hat{\mathbf{x}},t)\mathbf{\nabla}\cdot\mathbf{E}(\hat{\mathbf{x}})\\
&-\frac{\hbar e}{16m^2c^2}\epsilon^{ikl}\lbrace h_{ij}(\hat{\mathbf{x}},t),E_k(\hat{\mathbf{x}})(\hat{p}^j-\frac{e}{ c}A^j(\hat{\mathbf{x}}))\rbrace\Sigma_l+ \frac{\hbar e}{8m^2c^2}\epsilon^{ijl}\lbrace h_{0k}(\hat{\mathbf{x}}),{F_j}^k(\hat{p}_i-\frac{e}{ c}A_i)\rbrace\Sigma_l\\
&-\frac{\hbar^2 e}{8m^2c^2}h_{0j}(\hat{\mathbf{x}})\partial_i(F^{ij}(\hat{\mathbf{x}}))+\frac{i\hbar^2 e}{8m^2c^2}\epsilon^{ijl}h_{0k}(\hat{\mathbf{x}})\partial_i({F_j}^k(\hat{\mathbf{x}}))\Sigma_l-\frac{\gamma^0}{8m^3c^6}\Bigg[c^4\lbrace h_{00}(\hat{\mathbf{x}}),(\mathbf{\nabla}-\frac{e}{ c}\mathbf{A}(\hat{\mathbf{x}}))^4 \\
& +  c^4 \lbrace h_{ij}(\hat{\mathbf{x}}),(\hat{\mathbf{p}}-\frac{e}{ c}\mathbf{A}(\hat{\mathbf{x}}))^2(\hat{p}^i-\frac{e}{ c }A^i(\hat{\mathbf{x}}))(\hat{p}^j-\frac{e}{ c }A^j(\hat{\mathbf{x}}))\rbrace - \frac{\hbar^3 e c^3}{2}\epsilon^{ijl}h_{jm}(\hat{\mathbf{x}}){F_i}^m(\hat{\mathbf{x}})B^k(\hat{\mathbf{x}})\lbrace \Sigma_k ,\Sigma_l\rbrace \\
&+ 2\hbar^2 e c^2 h_{00}(\hat{\mathbf{x}})B^2(\hat{\mathbf{x}})+\frac{\hbar e c^3}{2}\epsilon^{ijl}\lbrace (\hat{\mathbf{p}}-\frac{e}{ c}\mathbf{A})^2 , h_{jk}(\hat{\mathbf{x}}){F_i}^k(\hat{\mathbf{x}}) \rbrace \Sigma_l-\hbar e c^3 \lbrace h_{00}(\hat{\mathbf{x}})\lbrace(\hat{\mathbf{p}}-\frac{e}{c}\mathbf{A})^2,B^k \rbrace\rbrace\Sigma_k \Bigg]\\
\end{split}
\end{equation}
\begin{equation}
\begin{split}
\hat{H}_d =&  -\frac{\hbar}{16m}\lbrace\partial_i(h_{00}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}^i-\frac{e}{ c}A^i(\hat{\mathbf{x}}))\rbrace\gamma^0 -\frac{\hbar^2}{16m}\partial^i\partial_i(h_{00}(\hat{\mathbf{x}}))\gamma^0 +\frac{i\hbar c}{4}\partial_{i}(h_0^i(\hat{\mathbf{x}}))+\frac{\hbar c}{4}\epsilon^{ijk}\partial_i(h_{0j}(\hat{\mathbf{x}}))\sigma_k\\ &-\frac{3i\hbar}{8}\partial_t(tr(h(\hat{\mathbf{x}})))+\frac{i\hbar}{4}\partial_t(h_{00}(\hat{\mathbf{x}})) +\gamma^0 \Bigg[ \frac{\hbar^2}{4m}\lbrace\partial^i(h_{00}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}_i-\frac{e}{c}A_i(\hat{\mathbf{x}}))\rbrace -\frac{\hbar^2}{8m}\lbrace\partial^i(h_{ij}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}^j-\frac{e}{c}A^j(\hat{\mathbf{x}}))\rbrace \\ &-\frac{\hbar^2}{4m}\lbrace\partial_i(\frac{tr(h(\hat{\mathbf{x}}))}{2}-h_{00}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}^i-\frac{e}{c}A^i(\hat{\mathbf{x}}))\rbrace-\frac{\hbar^2}{4m}\partial^i\partial_i (\frac{tr(h(\hat{\mathbf{x}}))}{2}-h_{00}(\hat{\mathbf{x}}))-\frac{e}{c}A^l(\hat{\mathbf{x}}))\rbrace \Big)\sigma_k\\
&-\frac{i\hbar^2}{8m}\epsilon^{ijk}\Big(\lbrace\partial_i(h_{00}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}_j-\frac{e}{c}A_j(\hat{\mathbf{x}}))\rbrace-\lbrace\partial_i(h_{jl}(\hat{\mathbf{x}})),\hat{\mathbf{p}}^l-\frac{e}{c}A^l(\hat{\mathbf{x}})\rbrace \Bigg]
\end{split}
\end{equation}
\end{widetext}
where $\hat{\mathbf{x}},\hat{\mathbf{p}}$ are respectively the single particle  position and the momentum operator. The term $\hat{H}_0$ is the usual Pauli Hamiltonian \cite{greiner} plus an irrelevant global phase $mc^2$ that can be reabsorbed with the transformation:
\begin{equation}
\vert\phi (t)\rangle \rightarrow e^{imc^2 t/\hbar}\vert\phi (t)\rangle
\end{equation}
The term $\hat{H}_r$ encodes the standard relativistic corrections~\cite{armin} due to the presence of the electromagnetic field up to order $\frac{v^4}{c^4}$. Finally, $\hat{H}_p$, $\hat{H}_d$ and $\hat{H}_{rp}$ account for the corrections due to the presence of the weak gravitational field respectively to $\hat{H}_0$ and $\hat{H}_r$.\\
\section{Differences with the bosonic model}\label{iv}
Equation~\eqref{hamf} is rather instructive as it extends the non relativistic Hamiltonian obtained in \cite{bosons} for the scalar field (namely Eq.(21) of \cite{bosons}) which we recall below~\footnote{We have added the superscripts (B) and (F) for respectively bosonic and fermionic in order to avoid confusion}:
\begin{equation}\label{hamB}
\begin{split}
{\hat{H}_0}^{(B)} =& mc^2 +\frac{\hat{\mathbf{p}}^2}{2m}\\
{\hat{H}_p}^{(B)} =& \frac{mc^2}{2}h^{00}(t,\hat{\mathbf{x}}) - \frac{\hbar^2}{8m}\lbrace h^{00}(t,\hat{\mathbf{x}}),\hat{\mathbf{p}}^2\rbrace + \frac{c}{2}\lbrace h^{0i}(t,\hat{\mathbf{x}}),\hat{p}_i\rbrace\\ &-\frac{1}{4m}\lbrace h^{ij}(t,\hat{\mathbf{x}}),\hat{p}_i\hat{p}_j\rbrace \\
{\hat{H}_d}^{(B)}=&\frac{\hbar^2}{8m}\mathbf{\nabla}^2(tr[h^{\mu\nu}(t,\hat{\mathbf{x}})])+\frac{i\hbar}{2}\partial_t(h^{00}(t,\hat{\mathbf{x}}))\\
&-\frac{i\hbar}{4}\partial_t(tr[h^{\mu\nu}(t,\hat{\mathbf{x}})])
\end{split}
\end{equation}
to describe the dynamics of a charged quantum particle (with spin 1/2) subject to a gravitational perturbation and to an external electromagnetic field. It is however as instructive to consider the electromagnetic free case , i.e. the limit $A(t,\hat{\mathbf{x}})\rightarrow 0$, in order to directly compare the fermionic and bosonic Hamiltonian. By taking the limit $A(t,\hat{\mathbf{x}})\rightarrow 0$ of Eq.~\eqref{hamf}, we obtain:
\begin{equation}\label{hamilF}
\begin{split}
{\hat{H}_0}^{(F)} =& mc^2 +\frac{\hat{\mathbf{p}}^2}{2m}\\
{\hat{H}_p}^{(F)} =& \frac{mc^2}{2}h^{00}(t,\hat{\mathbf{x}}) - \frac{\hbar^2}{8m}\lbrace h^{00}(t,\hat{\mathbf{x}}),\hat{\mathbf{p}}^2\rbrace + \frac{c}{2}\lbrace h^{0i}(t,\hat{\mathbf{x}}),\hat{p}_i\rbrace\\ &-\frac{1}{4m}\lbrace h^{ij}(t,\hat{\mathbf{x}}),\hat{p}_i\hat{p}_j\rbrace
\end{split}
\end{equation}
\begin{equation}
\begin{split}
&{\hat{H}_d}^{(F)}= -\frac{\hbar}{16m}\lbrace\partial_i(h_{00}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}^i-\frac{e}{ c}A^i(\hat{\mathbf{x}}))\rbrace\gamma^0\\ &-\frac{\hbar^2}{16m}\partial^i\partial_i(h_{00}(\hat{\mathbf{x}}))\gamma^0 +\frac{i\hbar c}{4}\partial_{i}(h_0^i(\hat{\mathbf{x}}))-\frac{3i\hbar}{8}\partial_t(tr(h(\hat{\mathbf{x}})))\\
&+\frac{\hbar c}{4}\epsilon^{ijk}\partial_i(h_{0j}(\hat{\mathbf{x}}))\sigma_k +\frac{i\hbar}{4}\partial_t(h_{00}(\hat{\mathbf{x}}))\\
&+\gamma^0 \Bigg[ \frac{\hbar^2}{4m}\lbrace\partial^i(h_{00}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}_i-\frac{e}{c}A_i(\hat{\mathbf{x}}))\rbrace \\
&-\frac{\hbar^2}{8m}\lbrace\partial^i(h_{ij}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}^j-\frac{e}{c}A^j(\hat{\mathbf{x}}))\rbrace \\ &-\frac{\hbar^2}{4m}\lbrace\partial_i(\frac{tr(h(\hat{\mathbf{x}}))}{2}-h_{00}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}^i-\frac{e}{c}A^i(\hat{\mathbf{x}}))\rbrace\\
&-\frac{\hbar^2}{4m}\partial^i\partial_i (\frac{tr(h(\hat{\mathbf{x}}))}{2}-h_{00}(\hat{\mathbf{x}}))\\
&-\frac{i\hbar^2}{8m}\epsilon^{ijk}\Big(\lbrace\partial_i(h_{00}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}_j-\frac{e}{c}A_j(\hat{\mathbf{x}}))\rbrace\\
&-\lbrace\partial_i(h_{jl}(\hat{\mathbf{x}})),(\hat{\mathbf{p}}^l-\frac{e}{c}A^l(\hat{\mathbf{x}}))\rbrace \Big)\sigma_k \Bigg]
\end{split}
\end{equation}
As expected, the bosonic and the fermionc description match for the gravity free case (${\hat{H}_0}^{(B)}={\hat{H}_0}^{(F)}$). They also match for the terms proportional to the gravitational perturbation $h_{\mu\nu}$. This is also to be expected: suppose in fact that there were actually a difference in the terms $h_{00}\hat{\mathbf{p}}^2/2m$ or $mc^2h_{00}$. This would imply that e.g. a simple change form Cartesian to Rindler \cite{mtw} coordinates would predict for a boson and a fermion to fall with the same acceleration in the first (Cartesian) but not the second (Rindler) reference frame, which would violate the weak equivalence principle.
The same line of reasoning can be applied to the other terms containing $h_{ij}$and $h_{0i}$. It is interesting however to notice that some differences arise for the terms containing the derivatives of the gravitational perturbation $\partial h_{\mu\nu}$. Such differences originated when we required the matter field to allow for a probabilistic interpretation in order to canonically quantize the system (see Eq. (\ref{eq:chargpres}) and Eq. (x) of \cite{bosons}).
\section{Master equation with electromagnetic field}\label{v}
In this section we derive a master equation to describe the decoherence effect induced by a weak stochastic gravitational perturbation on a spin 1/2 fermionic particle, as done for the scalar case in the previous chapter. For the sake of simplicity and compactness of the result, we will restrict our analysis to the Pauli Hamiltonian $\hat{H}_0$ and its gravitational corrections $\hat{H}_p$, as the terms $\hat{H}_r$ and $\hat{H}_{rp}$ are of higher order in the non relativistic expansion \footnote{Note however that one needs to be careful when applying the results of this section to a real experimental situations, as the term $\hat{H}_r$ might dominate over $\hat{H}_p$ (depending on the size of $\mathbf{E}$,$\mathbf{B}$, and $h_{\mu\nu}$) and should therefore be taken in consideration.}, and the term $\hat{H}_d$ contains derivatives of the gravitational perturbations, as in typical experimental situations~\cite{ligo,tian,decigo,et,kagra} they are negligible and would not add any further informative content to the analysis in any case.\\
This means that we approximate Eq.~\eqref{schrf} to:
\begin{equation}\label{schrf2}
i\hbar\partial_t\vert\phi(t)\rangle = (\hat{H}_0+\hat{H}_p)\vert\phi(t)\rangle
\end{equation}
If the metric is random, Eq.~\eqref{schrf2} becomes a stochastic differential equation. As a consequence the predictions are given by taking the stochastic average over the random gravitational field. We then need to specify its stochastic properties.\\
As done for the bosonic particle case, we assume the noise to be gaussian and with zero mean. For the sake of simplicity, we also assume the different components of the metric fluctuation to be uncorrelated. This means that the noise is fully characterized by:
\begin{equation}\label{stocpropf}
\begin{split}
\mathbb{E} [ h_{\mu\nu}(\mathbf{x},t)] =& 0\\
\mathbb{E} [h_{\mu\nu}(\mathbf{x},t)h_{\mu\nu}(\mathbf{y},s)] =& \alpha^2f_{\mu\nu}(\mathbf{x},\mathbf{y};t,s)
\end{split}
\end{equation}
where we recall that $\mathbb{E}[\, \cdot \,]$ denotes the stochastic average, $\alpha$ represents the strength of the gravitational fluctuations, and $f(\mathbf{x},\mathbf{y};t,s)$ is the two point correlation function.\\
We move to the density operator formalism, and write the von Neumann equation for the averaged density matrix :
\begin{equation}
\begin{split}\label{vonf}
\partial_{t}\hat{\rho}(t)=& -\frac{i}{\hbar}\Big[ \hat{H}_0(t),\hat{\rho}(t)\Big]-\frac{i}{\hbar}\mathbb{E}\Big[ [\hat{H}_p(t),\hat{\Omega}(t)]\Big]\\
 \equiv & \mathbb{E}\Big[ \mathfrak{L}[\hat{\Omega(t)}]\Big]
\end{split}
\end{equation}
We solve the above equation perturbatively exploiting the cumulant expansion~\cite{vank} (see Appendix C of \cite{bosons}). With the further help of the gaussianity, zero mean, uncorrelation of different components, we can rewrite Eq.~\eqref{vonf}  in Fourier space~\footnote{Recall that our choice for the Fourier transform is: $$f(\mathbf{x}) = \frac{1}{(\sqrt{2\pi}\hbar)^{3}}\int d^3q\:\tilde{f}(\mathbf{q}) e^{i\mathbf{q}\cdot\mathbf{x}/\hbar} $$} as:
\begin{widetext}
\begin{equation}\label{nonmarkovf}
\begin{split}
&\partial_t \hat{\rho} =  -\frac{i}{\hbar}[\hat{H}_0,\hat{\rho}(t)]+\\
&-\frac{\alpha^2}{\hbar^8}\int \frac{d^3q\:d^3q^\prime}{(2\pi)^{3}}\int_{0}^{t}dt_{1}\:\frac{\tilde{f}^{00}(\mathbf{q},\mathbf{q^\prime};t,t_1)}{4}\Big[ \Big\lbrace e^{i\mathbf{q}\cdot\hat{\mathbf{x}}/ \hbar},\Xi_{00}(\hat{\mathbf{x}},\hat{\mathbf{p}})\Big\rbrace,\Big[\Big\lbrace e^{i\mathbf{q^\prime}\cdot\hat{\mathbf{x}}_{t_1}/\hbar},\Xi_{00}(\hat{\mathbf{x}}_{t_{1}},\hat{\mathbf{p}})\Big\rbrace ,\hat{\rho}(t)\Big]\Big]  \\
&  -\frac{\alpha^2c^2}{\hbar^8}\int\frac{d^3q\:d^3q^\prime}{(2\pi)^{3}}\int_{0}^{t}dt_{1}\:\frac{\tilde{f}^{0i}(\mathbf{q},\mathbf{q^\prime};t,t_1)}{4}\Big[\Big\lbrace e^{i\mathbf{q}\cdot\hat{\mathbf{x}}/ \hbar},\hat{p}_i\Big\rbrace,\Big[\Big\lbrace e^{i\mathbf{q^\prime}\cdot\hat{\mathbf{x}}_{t_1}/\hbar},\hat{p}_i\Big\rbrace,\hat{\rho}(t)\Big]\Big]\\
& -\frac{\alpha^2}{\hbar^8}\int \frac{d^3q\:d^3q^\prime}{(2\pi)^{3}}\int_{0}^{t}dt_{1} \:\frac{\tilde{f}^{ij}(\mathbf{q},\mathbf{q^\prime};t,t_1)}{4}\Big[ \Big\lbrace e^{i\mathbf{q}\cdot\hat{\mathbf{x}}/ \hbar},\Xi_{ij}(\hat{\mathbf{x}},\hat{\mathbf{p}})\Big\rbrace,\Big[\Big\lbrace e^{i\mathbf{q^\prime}\cdot\hat{\mathbf{x}}_{t_1}/\hbar},\Xi_{ij}(\hat{\mathbf{x}}_{t_{1}},\hat{\mathbf{p}})\Big\rbrace ,\hat{\rho}(t)\Big]\Big]\\
&+O(t\alpha^3\tau_c^2)\end{split}
\end{equation}
\end{widetext}
where we have introduced:
\begin{equation}
\begin{split}
\Xi_{00} (\hat{\mathbf{x}},\hat{\mathbf{p}}) &= \frac{(\hat{\mathbf{p}}-\frac{e}{c}\mathbf{A}(t,\hat{\mathbf{x}}))^2}{4m}+\frac{m c^2}{2}-\frac{\hbar e}{2m c}\mathbf{B}(t,\hat{\mathbf{x}})\cdot\mathbf{\sigma}\\
\Xi_{ij}(\hat{\mathbf{x}},\hat{\mathbf{p}}) &= \frac{(\hat{p}_i-\frac{e}{c}A_i(t,\hat{\mathbf{x}}))(\hat{p}_j-\frac{e}{c}A_j(t,\hat{\mathbf{x}}))}{4m}+\frac{m c^2}{2}\\
&\quad+\frac{\hbar e}{2m c}\epsilon_{kil}{F^k}_j(t,\hat{\mathbf{x}})\sigma^l
\end{split}
\end{equation}
for the sake of compactness, and $\hat{x}_{t_1} = e^{i\hat{H}_0t_{1}}\hat{x}e^{-i\hat{H}_0t_{1}}$.\\
The above equation describes the dynamics of a pointlike spin 1/2 fermionic particle in presence of an external weak, stochastic gravitational field (with the further assumptions made in this section) and an external electromagnetic field.\\
\\
We specialize Eq.~\eqref{nonmarkovf} to the Markovian limit, i.e. we assume the noise to be delta correlated in time, with the further assumptions of isotropy and homogeneity of the noise, so that its correlation function reads:
\begin{equation}
f^{\mu\nu}(\mathbf{x},\mathbf{y};t,s) = \lambda u^{\mu\nu}(\mathbf{x}-\mathbf{y})\delta(t-s)
\end{equation}
where the factor $\lambda$ is in principle a generic coefficient with the dimension of a time. Note that the white noise assumption makes physical sense only if the correlation time ($\tau_{c}$) of the gravitational fluctuations is much smaller than the free dynamics' characteristic time  ($\tau_{free}$), or in the case where the contribution to the dynamics due to the gravitational perturbation is not affected by the free evolution dynamics, i.e. the operators describing the perturbation commute with the free dynamics operator $\hat{H}_0$. In such cases, as a first approximation, we can take $\lambda$ to be:
\begin{equation}\label{eq:lambda}
\lambda = \min(\tau_c \: ,\: t)
\end{equation}
Note that this choice does not affect the generality of the analysis as we leave $u^{\mu\nu}(\mathbf{x}-\mathbf{y})$ unspecified.\\
In such a regime Eq.~\eqref{nonmarkovf} is exact and it is easy to show that it reduces to:
\begin{widetext}
\begin{equation}\label{markovf}
\begin{split}
\partial_t \hat{\rho} = & -\frac{i}{\hbar}[\hat{H},\hat{\rho}(t)]\\
& -\frac{\alpha^2 \lambda}{(2\pi)^{3/2}\hbar^5 }\int d^3q \:\tilde{u}^{00}(\mathbf{q})\:\Big[\Big\lbrace e^{i\mathbf{q}\cdot\hat{\mathbf{X}}/ \hbar},\Xi_{00}(\hat{\mathbf{x}},\hat{\mathbf{p}})\Big\rbrace ,\Big[\Big\lbrace e^{-i\mathbf{q}\cdot\hat{\mathbf{X}}/\hbar},\Xi_{00}(\hat{\mathbf{x}},\hat{\mathbf{p}})\Big\rbrace ,\hat{\rho}(t)\Big]\Big]\\
&  -\frac{\alpha^2\lambda c^2}{(2\pi)^{3/2}\hbar^5 }\int d^3q \:\tilde{u}^{0i}(\mathbf{q})\:\Big[\Big\lbrace e^{i\mathbf{q}\cdot\hat{\mathbf{X}}/ \hbar},\hat{P}_i\big\rbrace,\Big[\Big\lbrace e^{-i\mathbf{q}\cdot\hat{\mathbf{X}}/\hbar},\hat{P}_i\Big\rbrace ,\hat{\rho}(t)\Big]\Big]\\
& -\frac{\alpha^2\lambda}{(2\pi)^{3/2}\hbar^5 }\int d^3q \:\tilde{u}^{ij}(\mathbf{q})\:\Big[\Big\lbrace e^{i\mathbf{q}\cdot\hat{\mathbf{X}}/ \hbar},\Xi_{ij}(\hat{\mathbf{x}},\hat{\mathbf{p}})\Big\rbrace,\Big[\Big\lbrace e^{-i\mathbf{q}\cdot\hat{\mathbf{X}}/\hbar},\Xi_{ij}(\hat{\mathbf{x}},\hat{\mathbf{p}})\Big\rbrace,\hat{\rho}(t)\Big]\Big
\end{split}
\end{equation}
\end{widetext}
Eq.~\eqref{markovf} describes decoherence in a complex combination of position momentum and energy bases, as it contains double commutators of functions of the position, momentum and free kinetic energy operators with the averaged density matrix.\\
\\
In what follows we will specialize Eq.~\eqref{markovf} to determine under which approximations it recovers decoherence in the position or momentum eigenbasis only.\\
As for the bosonic case \cite{bosons}, the conditions:
\begin{eqnarray}
\left\{
\begin{array}{ll}
h^{00}\gtrsim h^{0i}\\
h^{00}\gtrsim h^{ij} \\
\Delta E \ll Mc^2\:(1-u^{00}(\Delta \mathbf{x}))
\end{array}
\right.
\end{eqnarray}
are sufficient for our master equation to describe decoherence in the position eigenbasis only, where in this case the energy coherence needs to be modified to take into account the presence of the electromagnetic field, as $E= \frac{(\mathbf{p}-\frac{e}{c}\mathbf{A})^2}{2m}$. \\
Under the above assumptions, Eq.~\eqref{markovf} reads:
\begin{widetext}
\begin{equation}
\begin{split}
&\partial_t \hat{\rho} =  -\frac{i}{\hbar}[\hat{H},\hat{\rho}(t)] -\frac{\alpha^2 \lambda}{(2\pi)^{3/2}\hbar^5 }\int d^3q \:\tilde{u}^{00}(\mathbf{q})\Big[ e^{i\mathbf{q}\cdot\hat{\mathbf{X}}/ \hbar}\Big(\frac{mc^2}{2}-\frac{\hbar e \mathbf{\sigma}}{2m c}\cdot\mathbf{B}(t,\hat{\mathbf{x}})\Big),\Big[\Big\lbrace e^{-i\mathbf{q}\cdot\hat{\mathbf{X}}/\hbar}\Big(\frac{mc^2}{2}-\frac{\hbar e \mathbf{\sigma}}{2m c}\cdot\mathbf{B}(t,\hat{\mathbf{x}})\Big) ,\hat{\rho}(t)\Big]\Big]\\
\end{split}
\end{equation}
\end{widetext}
Contrary to the bosonic case, the condition of low momentum transfer
\begin{equation}
e^{i\mathbf{q}\cdot\hat{\mathbf{X}}/\hbar}\sim \hat{\mathbb{1}}
\end{equation}
is necessary but not sufficient to recover deocherence in the momentum or energy eigenbasis starting from Eq.~\eqref{markovf}. One in fact needs the further condition:
\begin{eqnarray}
\left\{
\begin{array}{ll}
|\mathbf{p}|\gg |\frac{e}{c}\mathbf{A}|\\
\frac{p^2}{2m}\gg|\frac{\hbar e \mathbf{\sigma}}{2mc}\cdot\mathbf{B}| \\
\end{array}
\right.
\end{eqnarray}
In this regime, Eq.~\eqref{markovf} can be approximated as:
\begin{equation}
\begin{split}
\partial_t \hat{\rho} = & -\frac{i}{\hbar}[\hat{H},\hat{\rho}(t)]\\
-& \frac{\alpha^2 \lambda}{(2\pi)^{3/2}\hbar^5 }\int d^3q \:\tilde{u}^{00}(\mathbf{q})\:\Big[\frac{\hat{\mathbf{p}}^2}{2m},\Big[\frac{\hat{\mathbf{p}}^2}{2m},\hat{\rho}(t)\Big]\Big]\\
-&  \frac{\alpha^2\lambda c^2}{(2\pi)^{3/2}\hbar^5 }\int d^3q \:\tilde{u}^{0i}(\mathbf{q})\:\Big[\hat{p}_i,\Big[\hat{p}_i,\hat{\rho}(t)\Big]\Big]\\
-& \frac{\alpha^2\lambda}{(2\pi)^{3/2}\hbar^5 }\int d^3q \:\tilde{u}^{ij}(\mathbf{q})\:\Big[\frac{\hat{p}_i\hat{p}_j}{2m},\Big[\frac{\hat{p}_i\hat{p}_j}{2m},\hat{\rho}(t)\Big]\Big]
\end{split}
\end{equation}
which indeed describes decoherence in the momentum eigenbasis.\\
\\
\section{Conclusions}
In this paper we have extended the results of our previous paper \cite{bosons} to the ferminonic case. We have  derived a
model of decoherence for non relativistic spin 1/2  particles interacting with both a weak stochastic gravitational perturbation and an external electromagnetic field. The resulting Hamiltonian and master equation correctly reproduce the results of the bosonic model in \cite{bosons} up to very small corrections. Such corrections account for relativistic effects (different spin of the two kind of particles) and the different quantization scheme employed.\\
The dynamics predicts also in this case decoherence in the position, momentum and energy eigenbasis, though under  different limiting cases than those described in \cite{bosons}.
\section*{Acknowledgments}
LA  thanks S. Ansoldi, L. Curcuraci, J.L. Gaona Reyes, C.I. Jones and S. Liberati for the helpful and inspiring discussions. The authors acknowledge financial support form the EU Horizon 2020 research and innovation program under Grant Agreement No. 766900 [TEQ].
GG acknowledges support from support from the Leverhulme Trust (RPG-2016- 046), the Spanish Agencia Estatal de Investigación, project PID2019-107609GB-I00 and QuantERA grant C'MON-QSENS!, by Spanish MICINN PCI2019-111869-2. LA thanks the University of Trieste and INFN.