\section{Introduction}

In general spacetimes the notion of particles loses its invariant meaning. One consequence of this observation is the detection of thermal particles in the Minkowski vacuum by a uniformly accelerated observer. This phenomenon is known as the Unruh effect \cite{Unruh:1976db}. It has immediate consequences for quantum states observed by accelerated observers. Their description becomes observer-dependent. In particular, the entanglement present in quantum states cannot be assumed to be independent of the observer \cite{Peres:2002wx, FuentesSchuller:2004xp}.

In the past, several investigations of entanglement in uniformly accelerating systems used nonlocalized Fock states, see e.g. \cite{FuentesSchuller:2004xp, Bruschi:2010mc, Richter:2015wha} and references therein. Therefore, it was not clear how to obtain an operational meaning of entanglement in these states. Furthermore, as pointed out in \cite{Dragan:2012hy, Ahmadi:2016fbd}, it is questionable that these simplified models capture the consequences of acceleration correctly. Another approach that is commonly used to study entanglement in the Minkowski vacuum is employing accelerated Unruh-de Witt detectors \cite{birrell1984quantum}. In these models, a two-level system is linearly coupled to a scalar field and the effect of finite time interactions is studied \cite{benatti2004entanglement}. In particular, it was shown that entanglement can be extracted from vacuum correlations \cite{Salton:2014jaa, Lorek:2014dwa}. The process of extracting entanglement from the vacuum sometimes is referred to as ``entanglement harvesting'' \cite{MartinMartinez:2012sg}. Since these effects are typically small, a direct observation is a challenging task. However, analog systems to study these and related phenomena were proposed \cite{delRey:2011jt, garcia2016entanglement, richter2017}. Another approach to localize quantum states was taken in \cite{Dragan:2012hy, Ahmadi:2016fbd, Kubicki:2016ngz, Dragan:2013dna}. Instead of making use of Unruh-de Witt detectors, the modes constituting the state itself were constructed to be localized. That is also the approach we are following in this work.

The goal of this work is to give a fully analytic treatment of the effect of acceleration on general localized fermionic Gaussian states. These states are of great relevance, as the class of fermionic Gaussian states contains a broad variety of states such as vacuum states of quadratic Hamiltonians, thermal states and Bell states. Our approach is based on a method that was developed to study bosonic two-mode Gaussian states of localized wave packets \cite{Ahmadi:2016fbd}. Similarly to the bosonic case, we are able to formulate the transformation connecting inertial and accelerated observers as the action of a Gaussian quantum channel.

The framework developed in this work is then used to study the entanglement of the vacuum. In particular, we quantify the amount of entanglement that could possibly be harvested from the vacuum by local detectors and show that it increases with increasing acceleration. Furthermore, we are able to obtain the entanglement in Bell states, a problem that attracted a lot of attention \cite{FuentesSchuller:2004xp, Bruschi:2010mc, Richter:2015wha}, and show that this effect is due to an inevitable mismatch of the wave functions in Minkowski and Rindler spacetime.

The structure of this work is the following. In Sec.\ \ref{framework} we introduce the setting that we study in this work. In Sec.\ \ref{fermionicnoninertial} we calculate the transformation of general fermionic Gaussian states for the case of two uniformly accelerated observers with arbitrary accelerations. With these results in hand, we study the vacuum entanglement in Sec.\ \ref{vacuumcase}. In particular, we compare our findings to the case of bosonic Gaussian states. In Sec.\ \ref{bellstates}, we study the entanglement in Bell states. Finally, we give the conclusions in Sec.\ \ref{conclusions}.

\section{Framework}\label{framework}

\subsection{Outline}

As anticipated above, we study the effect of acceleration on spatially localized fermionic modes. These modes are solutions of the equation of motion for the Dirac field $\hat{\Psi}$. The Dirac equation is obtained by varying the action
\begin{equation}\label{actiondirac}
S_D=\int \text{d}^2 x \sqrt{-g}\left(i\bar{\hat{\Psi}}g_{\mu\nu}\gamma^\mu \partial^\nu \hat{\Psi}-m \bar{\hat{\Psi}}\hat{\Psi}\right),
\end{equation}
where $g$ denotes the determinant of the metric $g_{\mu\nu}$, $\gamma^\mu$ are the gamma matrices, $m$ is the mass of the field and $\bar{\hat{\Psi}}=\hat{\Psi}^\dagger \gamma^0$. The field can be expanded in a complete set of solutions of the Dirac equation.  Since we are interested in modes that are sufficiently localized, we choose to expand $\hat{\Psi}$ in terms of localized wavepackets $\phi_k^\pm$. This gives the expansion
\begin{equation}\label{expansionminkowski}
\hat{\Psi}=\sum_k \left(\phi_k^+ \hat{f}_k +\phi_k^- \hat{g}_k^\dagger \right)
\end{equation}
in Minkowski space, where the $\hat{f}_k^\dagger/\hat{g}_k^\dagger$ are the creation operators for the wavepackets $\phi_k^\pm$ of particles/antiparticles. The Dirac equation dictates the anticommutativity of the creation and annihilation operators; $\{\hat{f}_k , \hat{f}_l\}= 0$, $\{\hat{f}_k^\dagger , \hat{f}_l^\dagger\}= 0$  and $\{\hat{f}_k,\hat{f}_l^\dagger\}= \delta_{kl}$, where the same holds for $\hat{g}_k$ and $\hat{g}_k^\dagger$ with vanishing mixed anticommutators. Alternatively, we can carry out an equivalent expansion in Rindler space and obtain in terms of the respective creation and annihilation operators $\hat{d}_k^\dagger/\hat{e}_k^\dagger$ and $\hat{d}_k/\hat{e}_k$
\begin{equation}\label{expansionrindler}
\hat{\Psi}= \sum_k \left(\psi_k^+ \hat{d}_k +\psi_k^- \hat{e}_k^\dagger\right).
\end{equation}
Here the $\psi_k^\pm$ denote the respective wavepackets and the $\hat{d}_k$ and $\hat{e}_k$ satisfy the same algebra as the $\hat{f}_k$ and $\hat{g}_k$.

Let Alice and Bob be two observers that are each in possession of a (localized) mode of the fermion field $\hat{\Psi}$. These modes have negligible overlap and we denote these by $\phi_\text{I}^+$ and $\phi_\text{II}^+$, respectively. It follows that the corresponding operators commute; $[\hat{f}_\text{I}, \hat{f}_\text{II}^\dagger]=0$. Initially, these modes are prepared in a fermionic Gaussian state $\rho$ that is completely characterized by its first and second moments (covariance matrix). Then to describe the state of the shared pair of modes from an accelerated perspective,  Rindler space offers a convenient reference frame. The modes in both spacetimes are related to each other by a Bogolyubov transformation. Hence, the Gaussianity of state $\rho$ is preserved under such a transformation. The resulting state of Alice's and Bob's modes is thus again a Gaussian of the transformed modes $\psi_\text{I}^\pm$ and $\psi_\text{II}^\pm$. As for bosons \cite{Ahmadi:2016fbd, PhysRevA.63.032312}, the map transforming the covariance matrix of the inertial modes, into the state of the accelerated modes, is a trace preserving Gaussian map $\sigma^{(f)}\rightarrow\sigma^{(d)}$, which takes the form
\begin{equation}\label{trafocovmatrix}
\sigma^{(d)}=M\sigma^{(f)} M^T +N,
\end{equation}
where $M$ and $N$ are $8\times 8$ matrices \cite{Bravyi05} and  we use the superscripts $(f)$ and $(d)$ to indicate the modes that are used to calculate the covariance matrix. One of the main results of the present work is the exact analytic calculation of the matrices $M$ and $N$ for the generic case of two uniformly accelerated observers.

\subsection{Fermionic Gaussian states}\label{fergaustates}

While bosonic Gaussian states are intensively studied in the context of quantum information and extensive literature exists \cite{RevModPhys.77.513, RevModPhys.84.621}, fermionic Gaussian states are much less investigated. Therefore, we insert a brief discussion of these states. The most prominent members of the family of fermionic Gaussian states are vacuum states of quadratic Hamiltonians, thermal states and the Bell states. Even so fermionic Gaussian states share some similarities with their bosonic counterparts, there are crucial differences due to the anticommutativity of fermions. Starting from the creation and annihilation operators $\hat{f}_k$ and $\hat{f}_k^\dagger$ that satisfy the CAR algebra, i.e, $\{\hat{f}_k , \hat{f}_l\}= 0$ and $\{\hat{f}_k,\hat{f}_l^\dagger\}= \delta_{kl}$, we define the Majorana fermion operators  $\hat{c}_k$, as
\begin{equation}\label{defofmajorana}
\hat{c}_{2j-1}=\frac{1}{\sqrt{2}}(\hat{f}_j^\dagger+\hat{f}_j), \hspace{5mm} \hat{c}_{2j}=\frac{1}{i\sqrt{2}} (\hat{f}_j^\dagger-\hat{f}_j).
\end{equation}
These form a Clifford algebra
\begin{equation}
\{\hat{c}_k,\hat{c}_l\}= \delta_{kl}.
\end{equation}
By definition, the density matrix of an even fermionic Gaussian state can be written as
\begin{equation}
\rho=C \,e^{-\frac{i}{2}\hat{c}^T A \hat{c}},
\end{equation}
where $A$ is a real, antisymmetric matrix and $C$ is the normalization \cite{Bravyi05}. That is, these states are thermal states of a quadratic Hamiltonian $H=\frac{i}{2}\hat{c}^T A \hat{c}$. Upon an $SO(2n)$ transformation $O$ with $n$ being the number of Majorana fermions $\hat{c}_i$, we can write the density matrix in the form
\begin{equation}
\rho=\frac{1}{2^n} \prod_{k=1}^n \left(1+i\lambda_k \hat{c}_{2k-1}' \hat{c}_{2k}'\right),
\end{equation}
where the $\lambda_i$ live in $[-1,1]$ and $\hat{c}'=O^T\hat{c}$. Further, using the Bloch-Messiah reduction \cite{bloch1962canonical}, we can find a basis of modes ${\hat{h}_k}$  such that every pure fermionic Gaussian state can be written in the form
\begin{equation}
\rho=\prod_{k=1}^n \left(u_k+v_k \hat{h}_{k}^\dagger \hat{h}_{-k}^\dagger\right),
\end{equation}
where $|u_k|^2+|v_k|^2=1$. Even Gaussian states are completely characterized by the corresponding covariance matrix $\sigma_{kl}$ \cite{Bravyi05}.  The real antisymmetric covariance matrix $\sigma_{kl}$ for a Gaussian state $\rho$ is given by
\begin{equation}\label{covmatrix}
\sigma_{kl}= i \,Tr(\rho[\hat{c}_k,\hat{c}_l])= i \langle \hat{c}_k \hat{c}_l - \hat{c}_l \hat{c}_k \rangle.
\end{equation}
The diagonal elements in (\ref{covmatrix}) are vanishing due to the commutator and the covariance matrix is completely characterized by the elements
\begin{equation}
\sigma_{kl}= 2i \langle \hat{c}_k \hat{c}_l \rangle,\hspace{10mm}\text{with} \,\, k>l.
\end{equation}
All higher moments of $\rho$ can be obtained by Wick's theorem \cite{Bravyi05}. So far we concentrated on even states. While these are completely characterized by $\sigma_{kl}$, the description of odd Gaussian states requires, in general, also the knowledge of the first moments
\begin{equation}\label{firstmoments}
 Tr(\rho \hat{c}_k)=  \langle \hat{c}_k \rangle.
\end{equation}
These are naturally vanishing for even states due to the vanishing commutator with the parity operator $P$, $[P, \rho]=0$ \cite{Bravyi05}. More information regarding fermionic Gaussian states and Gaussian linear maps can, for example, be found in \cite{kraus2009pairing, de2013power, bravyi2005classical,  Bradler:2010st, greplova2016degradability}.

\subsection{Dirac field in Rindler spacetime}

We briefly introduce solutions of the Dirac equation that is obtained by varying action (\ref{actiondirac}) with respect to $\bar{\hat{\Psi}}$. For simplicity, we restrict our attention to the $1+1$ dimensional case. However, the generalization to $3+1$ dimensions is straightforward. Furthermore, we use units such that $c=\hbar=1$ throughout the entire work. We start by deriving the solutions of the Dirac equation in Minkowski coordinates $(t,x)$. The Dirac equation takes the form
\begin{equation}
i\partial_t\hat{\Psi}=\left(-i\alpha_3\partial_x+m\beta\right)\hat{\Psi},
\end{equation}
where $m$ denotes the mass of the field and the matrices $\beta$ and $\alpha_3$ are given by
\begin{equation}\label{matrices}
\beta=\left( \begin{array}{cc}
1 & 0  \\
0& -1   \end{array} \right), \hspace{10mm} \alpha_3=\left( \begin{array}{cc}
0 & 1  \\
1& 0   \end{array}\right).
\end{equation}
We find the solutions in Minkowski coordinates to be given by
\begin{align}
u_{k,\pm}=&\frac{1}{\sqrt{4\pi\omega_k}}  \left(
\begin{array}{c}
\sqrt{\omega_k\pm m}\\
\pm\sqrt{\omega_k\mp m}
\end{array}
\right)e^{\mp i\omega_k t+ikx}  ,
\end{align}
where $\omega_k=\sqrt{k^2+m^2}$. These are normalized as given in (\ref{orthonormalization}). For further details, see \cite{birrell1984quantum}. Next, to describe the physics of an accelerated observer \cite{Crispino:2007eb, takagi1986vacuum}, we work in Rindler coordinates
\begin{subequations}\label{coordtrafo}
\begin{align}
t=& \chi \sinh(a\eta),\\
x=&\chi \cosh(a\eta),
\end{align}
\end{subequations}
 where $a$ is a parameter that does not have an immediate physical meaning. However, if we consider the worldline of a particle with a proper acceleration $\cal A$, namely $\chi=\frac{1}{\cal A}$, its proper time $\tau$ can be related
to the coordinate time $\eta$ as ${\cal A} \tau =a \eta$. In Rindler coordinates, the equation of motion obtained from (\ref{actiondirac}) is given by
\begin{equation}\label{diracinrindler}
i \frac{1}{a}\partial_\eta\hat{\Psi}=\left(m\chi \beta-\frac{i}{2} \alpha_3-i\chi\alpha_3\partial_{\chi}\right)\hat{\Psi},
\end{equation}
where $m$ denotes the mass of the field and the matrices $\beta$ and $\alpha_3$ are given by (\ref{matrices}). Solving the Dirac equation using the ansatz $w^\pm_{\text{I}\Omega}(\eta,\chi)= e^{\mp i\Omega\eta}w^\pm(\chi)$ for particles and antiparticles respectively, one finds the normalized solutions to be given by
\begin{widetext}
\begin{align}
w^\pm_{\text{I}\Omega}=&\sqrt{\frac{m \cosh(\frac{\pi \Omega}{a})}{2\pi^2 a}} \left(
\begin{array}{c}
K_{\pm i\frac{\Omega}{a}+\frac{1}{2}}(m\chi)+i K_{\pm i\frac{\Omega}{a}-\frac{1}{2}}(m\chi)\\
- K_{\pm i\frac{\Omega}{a}+\frac{1}{2}}(m\chi)+i K_{\pm i\frac{\Omega}{a}-\frac{1}{2}}(m\chi)
\end{array}
\right)  e^{\mp i\Omega\eta}, \hspace{13mm} \text{in Rindler wedge $I$},\\
w^\pm_{\text{II}\Omega}=&\sqrt{\frac{m \cosh(\frac{\pi \Omega}{a})}{2\pi^2 a}}  \left(
\begin{array}{c}
K_{\pm i\frac{\Omega}{a}+\frac{1}{2}}(-m\chi)+i K_{\pm i\frac{\Omega}{a}-\frac{1}{2}}(-m\chi)\\
- K_{\pm i\frac{\Omega}{a}+\frac{1}{2}}(-m\chi)+i K_{\pm i\frac{\Omega}{a}-\frac{1}{2}}(-m\chi)
\end{array}
\right) e^{\pm i\Omega\eta}, \hspace{6mm} \text{  in Rindler wedge $II$},
\end{align}
\end{widetext}
where $K_{ y}(x)$ is the modified Bessel function of the second kind of order $y$.  The Dirac inner product of the mode functions $\omega_1$ and $\omega_2$ is defined as $\left(\omega_1(x),\omega_2(x)\right)_\Sigma=\int \text{d}\Sigma^\mu \omega_1^\dagger(x)\gamma^0\gamma_\mu \omega_2(x)$, where $\text{d}\Sigma^\mu=n^\mu \text{d}s$ is space-like with normal vector $n^\mu$ and volume element $\text{d}s$. It has the properties $\left(\omega_1(x),\omega_2(x)\right)^*_\Sigma=\left(\omega_1^*(x),\omega_2^*(x)\right)_\Sigma =\left(\omega_2(x),\omega_1(x)\right)_\Sigma$. Here, in Rindler coordinates, we make the particular choice
\begin{equation}
\left(\omega_1,\omega_2\right)=\int \text{d}\chi\, \omega_1^\dagger\omega_2,
\end{equation}
where we suppressed the coordinate dependence of the mode functions. The positive and negative energy solutions satisfy the relations
\begin{equation}\label{orthonormalization}
(w^\pm_{\text{I}\Omega},w^\pm_{\text{I}\Omega'})= \delta( \Omega-\Omega'),\hspace{5mm}(w^\pm_{\text{I}\Omega},w^\mp_{\text{I}\Omega'})=0.
\end{equation}
The decomposition of $\hat{\Psi}$ in terms of localized modes is given by
\begin{equation}\label{fieldexpansion}
\hat{\Psi}=\sum_k \left(\phi_k^+ \hat{f}_k +\phi_k^- \hat{g}_k^\dagger \right)= \sum_k \left(\psi_k^+ \hat{d}_k +\psi_k^- \hat{e}_k^\dagger\right),
\end{equation}
where $\hat{f}_k/\hat{g}_k$ are the annihilation operators of particles/antiparticles in the Minkowski case and $\hat{d}_k/\hat{e}_k$ are the annihilation operators of particles/antiparticles in the Rindler case. These are related to the Rindler particle and antiparticle operators, $\hat{b}_\Omega$ and $\hat{a}_\Omega$,  by
\begin{subequations}\label{relationdb}
\begin{align}
\hat{d}_\text{I}=&\int \text{d}\Omega (\psi_\text{I}^+, w^+_{\text{I}\Omega})\hat{b}_{\text{I}\Omega},\\
\hat{d}_\text{II}=&\int \text{d}\Omega (\psi_\text{II}^+, w^+_{\text{II}\Omega})\hat{b}_{\text{II}\Omega},\\
\hat{e}_\text{I}^\dagger=&\int \text{d}\Omega (\psi_\text{I}^-, w^-_{\text{I}\Omega})\hat{a}^{\dagger}_{\text{I}\Omega},\\
\hat{e}_\text{II}^\dagger=&\int \text{d}\Omega (\psi_\text{II}^-, w^-_{\text{II}\Omega})\hat{a}^{\dagger}_{\text{II}\Omega}.
\end{align}
\end{subequations}
The Minkowski vacuum $|0\rangle_\text{M}$ is related to the Rindler vacuum $|0\rangle_\text{R}$ by a squeezing operator $S$ as
\begin{equation}
|0\rangle_\text{M}=S |0\rangle_\text{R},
\end{equation}
where $S$ acts on $\hat{b}_{\text{I}\Omega}$ and $\hat{a}_{\text{II}\Omega}$  as
\begin{subequations}\label{squeezing}
\begin{align}
S^\dagger \hat{b}_{\text{I}\Omega} S=& \cos(r_\Omega) \hat{b}_{\text{I}\Omega}-\sin(r_\Omega) \hat{a}_{\text{II}\Omega}^{\dagger},\\
S^\dagger \hat{a}_{\text{II}\Omega} S=& \cos(r_\Omega) \hat{a}_{\text{II}\Omega}+\sin(r_\Omega) \hat{b}_{\text{I}\Omega}^{\dagger}
\end{align}
\end{subequations}
with $\tan(r_\Omega)=e^{-\frac{\pi |\Omega|}{a}}$. To obtain the transformations of the operators $\hat{b}_{\text{II}\Omega}$ and $\hat{a}_{\text{I}\Omega}$, one simply has to interchange $\hat{a}\leftrightarrow \hat{b}$ in (\ref{squeezing}).

\subsection{Modes}\label{modes}

  \begin{figure}[t]

\includegraphics[scale=1]{figure1.pdf}
\caption{ Comparison between the real parts of the spatial modes $\phi^+$ and $\psi^+$ for the following choice of parameters: $x_0^{-1} = {\cal A}=0.1$, $L=2$, $\Omega_0=4.71$, $m=0.1$. These modes are localized and, therefore, a single proper acceleration ${\cal A}$ can be assigned.}
\label{picmodes}
\end{figure}

Most of our calculations are general and do not depend on the particular choice of modes. However, when considering concrete examples, namely in Sec.\ \ref{vacuumcase} and  Sec.\ \ref{bellstates}, we have to make a specific choice. Physically, cavity modes are suitable candidates for localized  solutions. In particular, motivated by the solutions obtained in \cite{Friis:2013eva}, for $\Lambda\in \{\text{I},\text{II}\}$, we consider the Minkowski modes
\begin{equation}\label{inercavity}
\phi_\Lambda^+(x,0)= C_\phi \, \xi_\phi(k) \, e^{-2\left(\frac{x_0}{L}\log(\frac{x}{x_0})\right)^2+ik x},
\end{equation}
where   $L$ is the width of the wave packet that is centred at $x_0$ and $\xi_\phi(k)$ is
\begin{equation}
\xi_\phi(k)=\left(
\begin{array}{c}
\cos\left(\frac{\kappa}{2}\right)+\sin\left(\frac{\kappa}{2}\right)\\
\cos\left(\frac{\kappa}{2}\right)-\sin\left(\frac{\kappa}{2}\right)
\end{array}
\right)  e^{-\frac{i}{2}\kappa}e^{-i k x_0}
\end{equation}
with $\kappa=\arctan(\frac{m}{k})$ and $C_\phi$ being a constant of normalization. As in \cite{Ahmadi:2016fbd}, we use a Gaussian profile to sufficiently localize the mode.

For the accelerated case, we choose wave packets of particle solutions of the Dirac equation in Rindler space, (\ref{diracinrindler}). Specifically, the localized modes we work with are
\begin{equation}\label{acccavity}
 \psi_\Lambda^+(\chi,0)= C_\psi \, \xi_\psi(\chi) \, e^{-2\left(\frac{\chi_0}{L}\log(\frac{\chi}{\chi_0})\right)^2},
\end{equation}
where  $\xi_\psi(\chi)$ is given by the spinor
\begin{equation}
\xi_\psi(\chi)=\left(
\begin{array}{c}
I_{i\frac{ \Omega }{a}-\frac{1}{2}}(m \chi )+i I_{i\frac{\Omega }{a}+\frac{1}{2}}(m \chi )\\
I_{i\frac{ \Omega }{a}-\frac{1}{2}}(m \chi )-i I_{i\frac{ \Omega }{a}+\frac{1}{2}}(m \chi )
\end{array}
\right)
\end{equation}
and $C_\psi$ is
\begin{equation}
C_\psi=\tilde{C}_\psi\left( I_{-i\frac{ \Omega }{a}-\frac{1}{2}}\left(m\chi_0\right)-I_{-i\frac{ \Omega }{a}+\frac{1}{2}}\left(m\chi_0\right)\right)
\end{equation}
with the constant of normalization, $\tilde{C}_\psi$ and $I_{i\mu}(x)$ being the modified Bessel functions of the first kind. The modes are plotted in Fig.\ \ref{picmodes}.

The physical motivation for choosing localized cavity modes is that these modes can possibly be detected by a local detector (in our case a cavity). Therefore, the choice of localized modes enables us to assign an operational meaning to quantities like entanglement. Furthermore, the modes (\ref{inercavity}) and (\ref{acccavity}) are purely positive frequency solutions. This guarantees, that the average particle and antiparticle numbers are zero in the Rindler vacuum, as expected. In consequence, the modes do necessarily have noncompact support. However, due to the Gaussian profile they are still sufficiently localized.

\section{Accelerated fermionic Gaussian states}\label{fermionicnoninertial}

In this section, we study the Bogolyubov transformation of fermionic Gaussian states  relating modes in Minkowski space to those of Rindler space. We give the analytical solution for the general case of two uniformly accelerated observers.

\subsection{Transformation of the Minkowski vacuum}

To study the effect of acceleration on the Minkowski vacuum, we restrict our attention to two localized fermion modes with respective creation operators $\hat{f}_1^\dagger$ and $\hat{f}_2^\dagger$. Then, using (\ref{defofmajorana}), we can  write the corresponding Majorana operators $\hat{c}_i$  as
\begin{equation}
\left(
\begin{array}{c}
\hat{c}_1\\
\hat{c}_2\\
\hat{c}_3\\
\hat{c}_4\\
\end{array}
\right)
=\frac{1}{\sqrt{2}}
\left(
\begin{array}{c}
\hat{f}_\text{I}^\dagger+\hat{f}_\text{I}\\
\frac{1}{i} (\hat{f}_\text{I}^\dagger-\hat{f}_\text{I})\\
\hat{f}_\text{II}^\dagger+\hat{f}_\text{II}\\
\frac{1}{i} (\hat{f}_\text{II}^\dagger-\hat{f}_\text{II})\\
\end{array}
\right).
\end{equation}
From (\ref{covmatrix}) it is clear that in the case of two modes, there are six independent entries in the fermionic covariance matrix. The same is true for the corresponding antiparticles. Therefore, we find the covariance matrix  of the Minkowski vacuum
\begin{align}\label{covmatrixinertial}
\sigma^{(f)}_\text{M}=&\left(\begin{array}{cccccccc}
0 & 1& 0 & 0  & 0 & 0& 0 & 0     \\
-1 & 0& 0 & 0  & 0 & 0& 0 & 0    \\
0 & 0& 0 & 1 & 0 & 0& 0 & 0     \\
0 & 0 & -1 & 0 & 0 & 0& 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 1& 0 & 0 \\
0 & 0 & 0 & 0 & -1 & 0& 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0& 0 &1 \\
0 & 0 & 0 & 0 & 0 & 0& -1 & 0 \\
\end{array}\right),
\end{align}
where this expression describes two particle modes and two antiparticle modes. However, it generalizes trivially to an arbitrary number of particles and antiparticle modes; see Appendix \ref{appminkvac} for details. As it is known for pure states \cite{Bravyi05}, $\sigma^{(f)^T}_\text{M} \sigma^{(f)}_\text{M}=1$ and the first moments vanish, as the vacuum is an even state.  To obtain  the covariance matrix $\sigma^{(d)}_\text{M}$ of the Minkowski vacuum in the Rindler frame, we need to transform the operators $\hat{f}$ and $\hat{g}$ (for details see Appendix \ref{appminkvactrafo}). We find the covariance matrix of the Minkowski vacuum in the Rindler basis to be of the form
\begin{equation}\label{covmatrixnoninertial}
\sigma^{(d)}_{}=\left(\begin{array}{cc}
\sigma_{+}^{(d)} & \sigma_{c}^{(d)}    \\
\tilde{\sigma}_{c}^{(d)}  & \sigma_{-}^{(d)}
\end{array}\right),
\end{equation}
where the block $\sigma_{+}^{(d)}$ captures the effect of acceleration on the respective mode and displays the thermal character of the Unruh effect ($\sigma_{-}^{(d)}$ is the analog for antiparticles). $\sigma_{c}^{(d)} $ and $\tilde{\sigma}_{c}^{(d)}$ describe the correlations between particles and antiparticles that arise due to the acceleration. Explicitly, (\ref{covmatrixnoninertial}) reads
\begin{widetext}
\begin{equation}\label{covmatrixnoninertialdntriv}
\sigma^{(d)}_{}=\left(\begin{array}{cccccccc}
0 & N_\text{I}^+& 0 & 0   & 0 & 0& \Im[N_\text{I,II}^+] & \Re[N_\text{I,II}^+]   \\
-N_\text{I}^+ & 0& 0 & 0   & 0 & 0& \Re[N_\text{I,II}^+] & -\Im[N_\text{I,II}^+]  \\
0 & 0 & 0 & N_\text{II}^+  &   -\Im[N_\text{I,II}^-] & -\Re[N_\text{I,II}^-]   & 0 & 0  \\
0& 0& -N_\text{II}^+ & 0  &  -\Re[N_\text{I,II}^-] & \Im[N_\text{I,II}^-]  & 0 & 0 \\
0 & 0& \Im[N_\text{I,II}^-] & \Re[N_\text{I,II}^-]  &  0 & N_\text{I}^-& 0 & 0  \\
0 & 0& \Re[N_\text{I,II}^-] & -\Im[N_\text{I,II}^-]  &  -N_\text{I}^- & 0& 0 & 0 \\
-\Im[N_\text{I,II}^+] & -\Re[N_\text{I,II}^+]   & 0 & 0  &  0 & 0 & 0 & N_\text{II}^- \\
-\Re[N_\text{I,II}^+] & \Im[N_\text{I,II}^+]  & 0 & 0& 0& 0& -N_\text{II}^- & 0
\end{array}\right),
\end{equation}
\end{widetext}
where the entries of the matrix are given by
\begin{subequations}\label{covmatelefermi}
\begin{align}
N_\text{I}^\pm=& 1-2\int \text{d}\Omega\, \frac{|(\psi_\text{I}^\pm , w^\pm_{\text{I}\Omega})|^2}{1+e^{\frac{2\pi\Omega}{a}}} \, ,\label{diagonalelement1}\\
N_\text{II}^\pm=& 1-2\int \text{d}\Omega\, \frac{|(\psi_\text{II}^\pm , w^\pm_{\text{II}\Omega})|^2}{1+e^{\frac{2\pi\Omega}{a}}}\,,\label{diagonalelement2}\\
N_\text{I,II}^\pm=&-2 \int \text{d}\Omega\, \frac{(\psi_\text{I}^\pm, w^\pm_{\text{I}\Omega}) (\psi_\text{II}^\mp, w^\mp_{ \text{II}\Omega})}{1+e^{\frac{2\pi\Omega}{a}}}\,e^{\frac{\pi\Omega}{a}}\label{diagonalelement3}
\end{align}
\end{subequations}
and we used $(\psi_\text{I}^\pm, w^\pm_{\text{I} \Omega})^*=(\psi_\text{I}^\mp, w^\mp_{\text{I} \Omega})$. We note that (\ref{covmatrixnoninertial}) reduces to  (\ref{covmatrixinertial}) in the limit of vanishing accelerations, $a\to 0$. That is, the correlated noise is absent in this limit.

\subsection{Transformation of general fermionic Gaussian states}\label{trafoofgeneralstates}

Using the results of the previous section, we derive the transformation for general fermionic Gaussian states. We recall that the transformed covariance matrix $\sigma^{(d)}$ can be obtained using (\ref{trafocovmatrix}).  Therefore, we can completely describe the transformation of general Gaussian states by the two matrices $M$ and $N$ that we calculate in the following. To obtain $M$, we consider the transformation of the first moments (\ref{firstmoments}). The Bogolyubov transformations can be obtained from (\ref{fieldexpansion}) and read
\begin{subequations}\label{gernaltrafo1}
\begin{align}
\hat{d}_l=&\sum_k (\psi_l^+, \phi_k^+) \hat{f}_k + (\psi_l^+, \phi_k^-) \hat{g}_k^\dagger ,\\
\hat{e}_l^\dagger=&\sum_k (\psi_l^-, \phi_k^+) \hat{f}_k +(\psi_l^-, \phi_k^-)  \hat{g}_k^\dagger  ,
\end{align}
\end{subequations}
where, in the following, we denote the overlaps by $\alpha_\text{I}^+=(\psi_\text{I}^+, \phi_\text{I}^+)$, $\beta_\text{I}^+=-(\psi_\text{I}^+, \phi_\text{I}^-)$, $\alpha_\text{I}^-=(\psi_\text{I}^-, \phi_\text{I}^-)^*$, $\beta_\text{I}^-=-(\psi_\text{I}^-, \phi_\text{I}^+)$, and accordingly for modes $\text{II}$. We note that both the Minkowski wave packets $\phi_{\text{I}/\text{II}}^\pm$ as well as the Rindler wave packets $\psi_{\text{I}/\text{II}}^\pm$ are composed of either only  particles or only  antiparticles. Therefore, the overlaps $\beta_{\text{I}/\text{II}}^\pm$ are vanishing and we can write (\ref{gernaltrafo1}) as
\begin{equation}\label{gernaltrafo2}
\hat{d}_l=\alpha_\text{I}^+ \hat{f}_l,\hspace{10mm} \hat{e}_l=\alpha_\text{I}^-  \hat{g}_l\, .
\end{equation}
The transformation of the first moments can be written as
\begin{equation}
\langle \hat{c}^{(d)}\rangle= M\langle \hat{c}^{(f)}\rangle
\end{equation}
with the block-diagonal $8\times 8$-matrix
\begin{equation}\label{Mmatrix}
M=\left(\begin{array}{cccc}
M(\alpha_\text{I}^+) & 0 &0&0     \\
0 & M(\alpha_\text{II}^+)  &0&0\\
0&0&M(\alpha_\text{I}^-)&0    \\
0&0&0&M(\alpha_\text{II}^-)
\end{array}\right),
\end{equation}
where we defined $M(\cdot)$ as the matrix given by
\begin{equation}
M(\cdot)=\left(\begin{array}{cc}
\Re[\cdot] & \Im[\cdot]      \\
-\Im[\cdot] & \Re[\cdot]      \\
\end{array}\right).
\end{equation}
We now gathered all quantities that we need to specify the transformation of a general Gaussian state characterized by $\langle \hat{c}^{(f)}\rangle$ and $\sigma^{(f)}$. First, we  explicitly calculate the noise matrix $N$. The matrix describing the noise arising due to acceleration is given by
\begin{equation}\label{noisematrix}
N=\sigma^{(d)}-M \sigma^{(f)}_\text{M} M^T,
\end{equation}
where $\sigma^{(d)}$ is given in (\ref{covmatrixnoninertial}). Therefore, in total we obtain
\begin{equation}\label{explnoise}
N =\left(\begin{array}{cc}
N_{+} & N_{\text{vac}}    \\
\tilde{N}_{\text{vac}}   & N_-
\end{array}\right),
\end{equation}
where $N_{\text{vac}}$ and $\tilde{N}_{\text{vac}}$ are given by the off-diagonal blocks of (\ref{covmatrixnoninertial}) and
\begin{widetext}
\begin{equation}
N_{\pm}=\left(\begin{array}{cccc}
0 & N_\text{I}^\pm-|\alpha_{\text{I}}^\pm|^2& 0 & 0     \\
-N_\text{I}^\pm+|\alpha_{\text{I}}^\pm|^2 & 0& 0 & 0    \\
0 & 0 & 0 & N_{\text{II}}^\pm-|\alpha_{\text{II}}^\pm|^2     \\
0& 0& -N_{\text{II}}^\pm+|\alpha_{\text{II}}^\pm|^2 & 0
\end{array}\right).
\end{equation}
\end{widetext}
With the matrices $M$ and $N$ we now have all data that is necessary to completely characterize the effect of acceleration to fermionic Gaussian states. In the following, we use the developed formalism to quantify vacuum entanglement and, as a example of a nonvacuum state, we study a maximally entangled state, the Bell state $|B\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)$.

\section{Vacuum entanglement}\label{vacuumcase}

\subsection{Vacuum entanglement for fermions}

In this section, we study the entanglement in the vacuum as seen by an accelerated observer. As a measure of entanglement we employ the logarithmic negativity. In \cite{eisler2015partial}, it was shown that the partial transpose of a fermionic  Gaussian state is, in general, not Gaussian and therefore it is difficult to calculate the negativity exactly. However, it was shown that a lower bound can be obtained \cite{eisert2016entanglement, eisler2015partial}. Building on this construction, we derive the respective bound  $\tilde{\mathcal{E}}_\mathcal{N}$ for our case in Appendix \ref{lognegativity}. It is given by
\begin{align}\label{lowerboundfermiond0}
\tilde{\mathcal{E}}_\mathcal{N}=&\ln(\frac{1}{2}(1+N_\text{I}^+N_\text{II}^-+|N_\text{I,II}^+|^2  +\nonumber\\
+& \Re[\sqrt{(N_\text{I}^+-N_\text{II}^-)^2-4|N_\text{I,II}^+|^2}]+\nonumber\\
+& \Im[\sqrt{(N_\text{I}^+-N_\text{II}^-)^2-4|N_\text{I,II}^+|^2}] ))
\end{align}
where $N_\text{I}^+$, $N_\text{II}^-$ and $N_\text{I,II}^+$ are the elements of the covariance matrix (\ref{covmatrixnoninertial}) that are given by (\ref{diagonalelement1}) - (\ref{diagonalelement3}).

\begin{figure}[t]
\includegraphics[scale=1.05]{figure2.pdf}
\caption{Negativity $\tilde{\mathcal{E}}_\mathcal{N}$ of the Minkowski vacuum  as a function of the proper accelerations. The parameters are $m=0.1$, $L=2$ and $\Omega_0=4.71$. The vacuum entanglement increases with increasing acceleration.}
\label{picnegvac}
\end{figure}

To obtain the entanglement of the vacuum, we choose the modes as described in Sec.\ \ref{modes} and solve (\ref{lowerboundfermiond0}) numerically. We show the results for the entanglement between two modes of different accelerations in Fig.\ \ref{picnegvac}.

Firstly, it is important to emphasize that $\tilde{\mathcal{E}}_\mathcal{N}$ is a lower bound for the negativity and therefore it is difficult to make strong quantitative statements. However, it is reasonable to assume that the actual value of the negativity is close to the lower bound $\tilde{\mathcal{E}}_\mathcal{N}$. Therefore, in the following, we refer to $\tilde{\mathcal{E}}_\mathcal{N}$ as negativity, keeping this subtlety in mind.

The entanglement between the state of a mode in Rindler wedge I and its respective counterpart in wedge II increases with increasing acceleration, see Fig.\ \ref{picnegvac}. For the parameter regime we are interested in, the entanglement is of the order $10^{-12}$, as measured by the logarithmic negativity. That is, the particles that are produced due to the Unruh effect are correlated across the acceleration horizon. Interestingly, the vacuum entanglement has an operational meaning, as it can, in principle, be extracted by suitable detectors \cite{Salton:2014jaa, Lorek:2014dwa}.

\subsection{Comparison to bosons}\label{bosons}

The aim of this section is to compare our results for fermionic states to the results in the  case of localized modes of a massive scalar field studied in \cite{Ahmadi:2016fbd}. For bosonic Gaussian states the covariance matrix $\sigma^{(d)}_{\text{bos}}$ of the vacuum takes the form
\begin{equation}\label{covmatrixnoninertialbosons}
\sigma^{(d)}_{\text{bos}}=
\left(\begin{array}{cccc}
N_\text{I}^b & 0& \Re[N_\text{I,II}^b] & \Im[N_\text{I,II}^b]      \\
0 & N_\text{I}^b& \Im[N_\text{I,II}^b] & -\Re[N_\text{I,II}^b]     \\
\Re[N_\text{I,II}^b] & \Im[N_\text{I,II}^b] & N_\text{II}^b &       \\
\Im[N_\text{I,II}^b] & \-Re[N_\text{I,II}^b] & 0 & N_\text{II}^b
\end{array}\right),
\end{equation}
where the matrix elements are given by
\begin{subequations}\label{covmatbosonele}
\begin{align}
N_{\text{I}/\text{II}}^b=& 1+2\int \text{d}\Omega \frac{|(\psi_\text{I/II}, w_{\text{I}/\text{II} \Omega})|^2 }{e^{\frac{2\pi\Omega}{a}}-1},\\
N_\text{I,II}^b=& 2\int \text{d}\Omega \frac{(\psi_\text{I}, w_{\text{I} \Omega}) (\psi_\text{II}, w_{\text{II} \Omega})}{e^{\frac{2\pi\Omega}{a}}-1}\,e^{\frac{\pi\Omega}{a}},
\end{align}
\end{subequations}
where $\psi_\text{I}$, $\psi_\text{II}$ are localized modes and the $w_{\text{I}/\text{II} \Omega}$ are solutions of the Klein-Gordon equation spanning the Hilbert space for regions I and II, respectively. The effect of acceleration on bosons is manifest in the covariance matrix (\ref{covmatrixnoninertialbosons}). While the Minkowski vacuum is described by the covariance matrix $\sigma^{(\text{bos})}_{M}=\mathbb{I}$, for nonvanishing acceleration, correlations build up and the off-diagonal elements increase.  As expected, the occupations are governed by the Bose-Einstein distribution that are explicit in the matrix elements (\ref{covmatbosonele}). This contrasts with the case of fermions where the occupation is characterized by the Fermi-Dirac distribution, (\ref{covmatelefermi}). Furthermore, due to the presence of antiparticles for fermions, the correlations that are build up in fermionic Gaussian states are between particles and antiparticles, while for bosons particles and antiparticles are identical.

Besides the differences outlined above, the results for fermions qualitatively agree with the findings for bosons. Quantitatively,  the entanglement we observed for fermions is less than in the bosonic case. In the parameter regime we studied, bosons develop ten times more entanglement \cite{Ahmadi:2016fbd}. However, one should keep in mind that we were calculating a lower bound for the negativity and that the amount of entanglement also depends on the explicit choice of modes.

In the next section, we go beyond vacuum entanglement and study the degradation of entanglement for Bell states.

\section{Entanglement in Bell states}\label{bellstates}

Interestingly, Bell states of fermions are fermionic Gaussian states, while maximally entangled states of two bosons are nonGaussian. This fact, enables us to apply our formalism to Bell states. Therefore, as an example for a nonvacuum state, we consider a maximally entangled state $|B\rangle$ of two particles in the following. We choose
\begin{equation}\label{bellstate}
|B\rangle=\frac{1}{\sqrt{2}}(|00\rangle+|11\rangle)=\frac{1}{\sqrt{2}}(1+\hat{f}_1^\dagger \hat{f}_2^\dagger)|0\rangle.
\end{equation}
This state is an even Gaussian state and it is described by the following covariance matrix (cf.\ Appendix \ref{appminkbell})
\begin{equation}
\sigma^{(f)}_{\text{Bell}}=\left(\begin{array}{cccccccc}
0 & 0&0 &  1 & 0& 0& 0& 0    \\
0 & 0& 1& 0 & 0& 0& 0& 0   \\
0 & - 1 & 0 & 0  & 0& 0& 0& 0  \\
- 1 & 0 & 0 & 0& 0& 0& 0& 0   \\
0 & 0 & 0 & 0& 0& 1& 0& 0   \\
0 & 0 & 0 & 0& -1& 0& 0& 0   \\
0 & 0 & 0 & 0& 0& 0& 0& 1   \\
0 & 0 & 0 & 0& 0& 0& -1& 0
\end{array}\right).
\end{equation}
The transformed state is obtained according to equation (\ref{trafocovmatrix}) with the matrix $M$ given in (\ref{Mmatrix}) and $N$ given in (\ref{explnoise}). Then the transformed covariance matrix, neglecting some terms much smaller than 1, is given by
\begin{widetext}
\begin{equation}\label{transformedcovmatbell}
\sigma^{(d)}_\text{Bell}=\left(\begin{array}{cccc}
0 & N_\text{I}^+-|\alpha_\text{I}^+|^2 & (\sigma^{(d)}_\text{Bell})_{13} & (\sigma^{(d)}_\text{Bell})_{14}   \\
_-N_\text{I}^++|\alpha_\text{I}^+|^2& 0&   (\sigma^{(d)}_\text{Bell})_{23}   &     (\sigma^{(d)}_\text{Bell})_{24}     \\
- (\sigma^{(d)}_\text{Bell})_{13} & -  (\sigma^{(d)}_\text{Bell})_{23} & 0 & N_\text{II}^+ -|\alpha_\text{II}^+|^2    \\
-(\sigma^{(d)}_\text{Bell})_{14}  & -  (\sigma^{(d)}_\text{Bell})_{24} & -N_\text{II}^++|\alpha_\text{II}^+|^2 & 0
\end{array}\right),
\end{equation}
\end{widetext}
where
\begin{subequations}
\begin{align}
(\sigma^{(d)}_\text{Bell})_{13}=&\Im[\alpha_\text{II}^+] \Re[\alpha_\text{I}^+] + \Im[\alpha_\text{I}^+] \Re[\alpha_\text{II}^+],\\
(\sigma^{(d)}_\text{Bell})_{14}=& \frac{1}{2} \left(\alpha_\text{I}^+ \alpha_\text{II}^+ + (\alpha_\text{I}^+\alpha_\text{II}^+)^*\right), \\
(\sigma^{(d)}_\text{Bell})_{23}=&\frac{1}{2} \left(\alpha_\text{I}^+ \alpha_\text{II}^+ + (\alpha_\text{I}^+\alpha_\text{II}^+)^*\right),\\
(\sigma^{(d)}_\text{Bell})_{24}=&-\Im[\alpha_\text{II}^+] \Re[\alpha_\text{I}^+] - \Im[\alpha_\text{I}^+] \Re[\alpha_\text{II}^+]
\end{align}
\end{subequations}
with $\alpha_{\text{I}/\text{II}}^+$ as defined in (\ref{gernaltrafo2}). In (\ref{transformedcovmatbell}), we neglected contributions from the vacuum noise matrix $N$, as these are insignificant corrections to $\sigma^{(d)}_\text{Bell}=M\sigma^{(f)}_\text{Bell} M^T$.

To quantify the entanglement we employ a lower bound for the negativity given by
\begin{align}\label{lowerboundfermionbell}
\tilde{\mathcal{E}}_\mathcal{N}=\ln(\frac{1}{2}+&\frac{\Re[\nu^+\nu^- +\nu^+-\nu^-]}{2}+\nonumber\\
+&\frac{\Im[\nu^+\nu^-+\nu^+-\nu^-]}{2}),
\end{align}
where the $\nu^\pm$ can be calculated from (\ref{transformedcovmatbell}); cf.\ (\ref{lowerboundneg}) in Appendix \ref{lognegativity}.

\begin{figure}
\centering
\includegraphics[scale=1]{figure3.pdf}
\caption{Negativity $\tilde{\mathcal{E}}_\mathcal{N}$ of a Bell state as a function of the proper acceleration $\mathcal{A}=\mathcal{A}_\text{I}=\mathcal{A}_\text{II}$. The parameters are $m=0.1$, $\Omega_0=4.71$ and $L=2$. We normalized $\tilde{\mathcal{E}}_\mathcal{N}$ such that for vanishing accelerations $\tilde{\mathcal{E}}_\mathcal{N}=1$. The blue continuous line gives the numerical results, while the gray dotted line is an extrapolation to $\mathcal{A}=0$. From an inertial perspective Bell states are maximally entangled. However, for an accelerated observer, entanglement gets degraded with increasing acceleration.}
\label{picbell}
\end{figure}

We show the results for the entanglement of state $|B\rangle$ in Fig.\ \ref{picbell}. For vanishing acceleration the modes are maximally entangled, while the entanglement decreases with increasing acceleration. The constraint $\mathcal{A}L\lesssim 1$ that ensures that a single acceleration can be associated to a localized mode prevents us from studying arbitrarily large accelerations.

The physical reason for the degradation of entanglement is that the overlap between the inertial wavepackets $\phi^+_k$ and the accelerated ones $ \psi_\Omega^+$ decreases with increasing acceleration. This is an inevitable effect of the Bogolyubov transformation connecting Minkowski and Rindler solutions. Therefore, with increasing acceleration the mismatch between the modes increases and the entanglement in the initially maximally entangled Bell state decreases. Interestingly, the effect of Unruh radiation is negligibly small. That is, the production of thermal particles due to the Unruh effect has only a very small effect to entanglement degradation of states of localized modes.

\section{Conclusions}\label{conclusions}

In this work, we developed a general framework to describe the effect of acceleration on arbitrary fermionic Gaussian states of two localized modes. We formulated the transformation between the quantum state seen by initial observers and the state seen by uniformly accelerated ones as the action of a fermionic Gaussian channel, and completely characterized this channel. This enabled us to study how the entanglement of the vacuum and the entanglement in Bell states is affected by acceleration.

We found that vacuum entanglement is enhanced by acceleration. That is, the correlations in the vacuum that build up due to the Unruh effect lead to entanglement between particles and antiparticles. In particular, as our framework is employing localized modes, entanglement has an operational meaning and can be extracted by suitable local detectors; a process that is sometimes referred to as ``entanglement harvesting''.  We also compare these findings  with the case of bosons \cite{Ahmadi:2016fbd} and find a qualitative agreement.

Furthermore, we quantified the entanglement degradation in a localized maximally entangled state of fermions. The degradation is due to an increasing mismatch between the initial modes and the accelerated modes that causes their overlap to decrease and, therefore, entanglement decreases as well. The effect of particle creation on entanglement degradation is negligible.

We emphasize that due to the localization of the modes, the framework presented in this work can be applied to quantum information protocols, in scenarios where local gravitational effects or effects due to acceleration are not negligible. In particular, the entanglement can be exploited as a resource.

For the future, we are interested in extending our studies to localized modes of fermions with arbitrary mutual separation, as well as studying localized Gaussian states in higher dimensional spacetimes.

\begin{acknowledgements}

B.R. thanks the Institute of Theoretical Physics at the University of Warsaw for hospitality during his visit. B.R. and Y.O. acknowledge support from Funda\c{c}\~{a}o para a Ci\^{e}ncia e a Tecnologia (Portugal), namely through programmes PTDC/POPH/POCH and projects UID/EEA/50008/2013, IT/QuSim, IT/QuNet, ProQuNet, partially funded by EU FEDER, from the EU FP7 project PAPETS (GA 323901), and from the JTF project NQN (ID 60478). Furthermore, B.R. acknowledges the support from the DP-PMI and FCT through scholarship SFRH/BD/52651/2014. K.L. and A.D. would like to acknowledge financial support from the National Science Center, Sonata BIS Grant No. DEC-2012/07/E/ST2/01402.

\end{acknowledgements}


\section{Introduction}

Traditionally, wireless services have been offered by Service Providers (SPs) that own the infrastructure they are operating on.  Nowadays SPs are divided into (i) Mobile Network Operators (MNOs) that own the infrastructure, and (ii)
Mobile Virtual Network Operators (MVNOs)  that do not own the infrastructure they are operating on, and use the resources of one or more MNOs based on a business contract. MVNOs can distinguish their plans from MNOs by  bundling their service with other products,  offering different pricing plans  for End-Users (EUs), or building  a good reputation through a better customer service. In recent years, the number of MVNOs is rapidly growing. According to \cite{GSMA}, between June 2010 and June 2015, the number of MVNOs increased by 70 percent worldwide,  reaching 1,017 as of June 2015.  Even some MNOs developed their own MVNOs. An example of which is Cricket wireless which is owned by At\&t and offers a prepaid wireless service to EUs.

In this work, we consider the economics of the interaction between  MVNOs and MNOs.  We investigate the incentives of an MNO for offering some of her resources to an MVNO instead of using the resources for her own. Thus, we consider competition and cooperation between an MNO  and an MVNO. Note that it is not apriori known  how much an MNO is willing to invest on the infrastructure and how much an MVNO is willing to lease. More importantly, it is not apriori clear that under what conditions the MNO prefers to generate her revenue through the MVNO by leasing the resources to her, and therefore letting the MVNO to attract EUs.

Many works have considered the economics of resource/spectrum sharing and the subsequent profit sharing between SPs. Examples are \cite{cano2016cooperative,berry2013nature,le2009pricing,singh2012cooperative,korcak2012collusion
,banerjee2009voluntary,lotfi2016uncertain,sponsoring_journal}. In these works, authors model the environment using game theory and seek to provide intuitions for the pricing schemes, the split of EUs/benefits between SPs, and the investment level of different SPs. In this work, however, as mentioned before, we focus on the competition between the MNOs and MVNOs.

We consider a market with one MNO (e.g. At\&t) and  one MVNO (e.g. Google's Project Fi). We assume that the MVNO is active, i.e. has already some resources, and is  willing to lease additional resources from the MNO   in exchange of a fee. The MNO decides on the investment level. Subsequently, knowing this decision, the MVNO decides on the number of resources she wants to lease from the MNO. Then, simultaneously, both MNO and MNVO decide on their pricing strategies for EUs.   Finally, EUs choose one of the MNO or MVNO to buy the wireless plan from. We assume the per resource fee that the MVNO pays to the MNO to be fixed and discuss about the framework for determining this fee.

To model EUs, we consider a continuum of undecided EUs that decide which of the SPs they want to buy their wireless plan from. We assume that EUs have different preferences for each SP. These preferences can be because of different services that SPs offer. For example, the MVNO can bundle the wireless service with a free or cheap international call plan (through VoIP) to make her service more favorable for some EUs. Moreover, the preferences can be because of the reluctance of  an EU to buy her wireless plan from a particular SP (e.g. an  SP with an infamous customer service). We assume that the preferences also depend on the investment level of the SPs. In other words, the higher the investment level of an SP, the lower would be the reluctance of an EU for that SP.

Different preferences of EUs for SPs enables the possibility of different outcomes for the market. For example, the MVNO who bundles her service with free international call may attract some of EUs even with a higher price than what MNO is offering. Thus,  instead of competing for EUs (by lowering her price), the MNO can lease some of her resources to this MVNO and receives her share of profit through the MVNO.

We formulate the game as a sequential game, and seek the
Sub-game Perfect Nash Eequilibrium (SPNE) of the sequential
game using backward induction. We show that there exists a unique interior  SPNE, i.e. the SPNE by which both SPs receive a positive mass of EUs, and characterize it.  Moreover, we consider a benchmark case in which the MNO does not invest in her infrastructure, and subsequently the MVNO cannot lease any resources from the MNO. We  characterize the unique SPNE of this case. We  compare the results of our model to the benchmark case to assess the incentive of the MNO to invest in her infrastructure and to offer it to the MVNO.

Analytic and numerical results reveal that although the number of EUs and the price that the MNO quotes for EUs might be lower than those of the MVNO, the MNO would still generate a higher payoff than the MVNO if the per resource fee that the MVNO pays to the MNO is  sufficiently higher than the marginal cost of investment on the infrastructure.  In addition, results reveal that if the MNO knows that EUs are reluctant joining her or have a high preference for the MVNO, then a better strategy for her would be to invest in her infrastructure, lease a portion (or the whole) of the new resources to the MVNO, and collect the benefit through the fees she charges for the resources leased to the MVNO.

The rest of the paper is organized as follows. In Section~\ref{section:model}, we present the model. In Section~\ref{section:analysis}, we find the SPNE strategies. We present the SPNE outcome in Section~\ref{section:outcome}. In Section~\ref{section:benchmark}, we prove the results of the benchmark case. Finally, in Section~\ref{section:numerical}, we present the numerical results and discuss about the results.

\section{Model}\label{section:model}
We consider one MNO and one MVNO that compete in attracting a pool of undecided EUs.  We assume that the MVNO is active, i.e. has already some resources. Thus, the model is set up such that even when the MVNO does not lease any resources from the MNO, it can still attract some of EUs.

\subsection*{SPs:}
We denote the MNO by SP$_{L}$ ($L$ represents Leader, since this SP is the leader of the game), and denote the MVNO by SP$_{F}$ ($F$ represents Follower, since this SP is the follower of this leader/follower game). SP$_L$ owns the infrastructure, invests in her infrastructure to attract EUs, and can lease parts of the new resources to SP$_F$ in exchange of money.

We denote the fee per resources that SP$_{F}$ pays to SP$_{L}$ by $s$. The utility of SP$_F$ is increasing with the revenue from EUs, which depends on the number of EUs ($n_F$) and the access fee ($p_F$). This utility is decreasing with  respect to the fee she pays to SP$_L$ to reserve resources, i.e. $sI^2_F$, where $I_F$ is the number of resources that SP$_F$ reserves from SP$_L$:
\begin{equation} \label{equ:payoffF}
\pi_F=n_F(p_F-c)-sI^2_F
\end{equation}
, where $c$ is the marginal cost associated with each user. Note that naturally, we expect the cost of investing in infrastructure to be convex, i.e. the cost of investment per resources increases with the amount of resources. For simplicity in analysis, we consider the cost of reserving resources  to be  quadratic with respect to $I_F$.\footnote{The overall intuitions of the model are expected to be the same in the case of considering any convex function of $I_F$}

The utility of SP$_L$ is:
\begin{equation}\label{equ:payoffL}
\pi_L=n_L(p_L-c)+sI^2_F-\gamma I^2_L
\end{equation}
where $n_L$, $p_L$, and $I_L$ are the number of EUs with SP$_L$, the access fee quoted by SP$_L$ for EUs, and the number of resources that SP$_L$ add to her infrastructure. Note that $\gamma$ is the marginal cost of investment. Thus, the utility of SP$_L$ is increasing with the revenue from EUs
($n_L(p_L-c)$) and the fee received from SP$_F$ ($sI^2_F$), and  is decreasing with  respect to the investment cost ($\gamma I^2_L$). Note that we consider the cost of investment to be a quadratic function of $I_L$.\footnote{Intuitions of the model are expected to be the same when considering any convex function of $I_L$.}

Trivially, we assume that the fee per resources ($s$) is greater than or equal to the marginal cost of investment on the infrastructure ($\gamma$), i.e. $s\geq \gamma$. Also, note that $I_F\leq I_L$. To have a non-trivial problem, we also assume $I_L>0$ and $I_F\geq 0$.

The strategies of SP$_L$ is to choose the access fee for EUs ($p_L$) and the level of investment ($I_L$). The strategies of SP$_F$ is to choose the access fee for EUs ($p_F$) and the number of resources she leases from SP$_L$ ($I_F$). We assume that the per resource fee $s$ is pre determined, possibly through a bargaining framework between SP$_L$ and SP$_F$.

\subsection*{EUs:}  The strategy of an EU is to choose one of the SPs to buy wireless plan from. We model the EUs using a hotelling model. We assume that the SP$_L$ is located at 0, SP$_F$ is located at 1, and EUs are distributed uniformly along the unit interval $[0,1]$. The closer an EU to an SP, the more this EU prefers this SP to the other. Note that the notion of closeness and distance is used to model the preference of EUs, and may not be the same as  physical distance.

Specifically, the EU located at $x\in[0,1]$ incurs a \emph{transport cost} of $t_L x$  (respectively, $t_F(1-x)$) when joining SP$_L$ (respectively, SP$_F$), where $t_L$ (respectively, $t_F$) is the marginal transport cost for SP$_L$ (respectively, SP$_F$). In sum, we  consider $t_L$ and $t_F$ as the reluctance of EUs for connecting to SP$_L$ and SP$_F$, respectively.

We  assume the transport costs of EUs to depend on the investment level of  SPs. Thus, we assume that the higher the investment level of an SP in comparison to the other SP, the higher would be the reluctance of EUs for joining the other SP. This model is used to captures the factor of quality in the decision of EUs. Thus, the transport cost of  SP$_L$ and SP$_F$ are considered to be $t_L=\frac{I_F}{I_L}$ and $t_F=1-t_L=\frac{I_L-I_F}{I_L}$, respectively.

Therefore, the utility of an EU located at distance $x_j$ of  SP$_{j}$, $j\in\{L,F\}$ is:
\begin{equation}
u_{x_j}=v^*-t_jx_j-p_j
\end{equation}
, where $v^*$ is a common valuation that captures the  value of buying a wireless plan for EUs regardless of the SP chosen.

\subsection*{Formulation:}
We model the interaction between SPs and EUs using a four-stage sequential game. Naturally, we assume that SP$_L$ makes the first move and is the leader of the game. The timing and the stages of the game are as follows:
\begin{enumerate}
\item SP$_{L}$ decides on the investment on the infrastructure ($I_L$).
\item SP$_{F}$ decides on the investment, i.e. number of resources to lease from SP$_L$ ($I_F$ such that $I_F\leq I_L$).
\item SP$_{L}$ and SP$_{F}$ determine the access fees for EUs (respectively, $p_L$ and $p_F$).
\item EUs decide to subscribe to one of the SPs.
\end{enumerate}

We assumed the selection of investments by SPs ($I_L$ and $I_F$) to happen before the selection of prices for EUs ($p_L$ and $p_F$) since investment decisions are long-term decision. These decisions are expected to be kept constant over longer time horizons in comparison to pricing decisions.

In the sequential game framework,  we seek a \emph{Subgame Perfect Nash Equilibrium}  using \emph{backward induction}:

\begin{definition}\emph{Subgame Perfect Nash Equilibrium (SPNE):}
	A strategy is an SPNE if and only if it constitutes a Nash Equilibrium (NE) of every subgame of the game.
\end{definition}

\begin{definition}\emph{Backward Induction:} Characterizing the equilibrium strategies starting from the last stage of the game and proceeding backward.
\end{definition}

We also assume the full market coverage of EUs by SPs. This means that  each EU chooses exactly one SP to subscribe to.  This  assumption is common in hotelling models and is necessary to ensure competition between SPs. An equivalent assumption is to consider the common valuation $v^*$ to be sufficiently large so that the utility of EUs for buying a wireless plan  is positive regardless of the choice of SP.

 One might think that SP$_L$ can eliminate competition by not offering her resources to SP$_F$ (or simply not investing in  new infrastructure), and use the full market coverage to act as a monopoly. However, this is not the case in our model. Recall that we assumed the SP$_F$ to have resources independent of SP$_L$. Thus, regardless of the strategy of SP$_L$, SP$_F$ would be still active in the market, and the full market coverage would not lead to the monopoly of SP$_L$.

\section{The Sub-Game Perfect Nash Equilibrium}\label{section:analysis}

In the backward induction, we start with Stage 4:

\subsection*{Stage 4:}

In this subsection, we characterize the division of EUs between SPs in the  equilibrium, i.e. $n_L$ and $n_F$, using the knowledge of the strategies chosen by the  SPs in Stages 1, 2, and 3. To do so, we characterize the location of the EU that is indifferent between joining either of the SPs, $x_n$. Thus, EUs located at $[0,x_n)$ join SP$_L$, and those located at $(x_n,1]$ joins the non-neutral SP$_F$ (using the full market coverage assumption).

The EU located at $x_n\in[0,1]$ is indifferent between connecting to SP$_L$ ans SP$_F$ if:
\begin{equation}\label{equ:xn_t_parameter}
\small
	\begin{aligned}
	&v^*-t_{F}(1-x_n)-p_{F}=v^*-t_Lx_n-p_L\\
	&\qquad \qquad \qquad  \Rightarrow x_n=\frac{t_F+p_F-p_L}{t_L+t_F}
	\end{aligned}
\end{equation}

Note that $t_F+t_L=1$. Substituting the value of $t_F$ yields:
\begin{equation}\label{equ:xn}
\small
x_n=\frac{I_L-I_F}{I_L}+p_F-p_L
\end{equation}

Thus, the fraction of EUs with each SP ($n_L$ and $n_F$) is:
\begin{equation} \label{equ:EUs_stage4}
\small
\begin{aligned}
n_L &=
\left\{
\begin{array}{ll}
	0  & \mbox{if } x_n < 0 \\
	\frac{I_L-I_F}{I_L}+p_F-p_L & \mbox{if } 0\leq x_n \leq 1\\
	1 & \mbox{if } x_n>1
\end{array}
\right. \\
n_F&=1-n_L
\end{aligned}
\end{equation}

\subsection*{Stage 3:}

In this stage, SP$_L$ and SP$_F$ determine their  prices for EUs, $p_L$ and $p_F$, respectively, to maximize their payoff. We seek for Nash equilibrium (NE) strategies. Note that in general there might exist several NE strategies: some of them corner equilibria (an extreme case in which one of the SPs receives zero EUs) and some interior equilibria (in which both SPs receive a positive mass of EUs). In practice, the latter equilibria are expected to occur more frequently. Thus, henceforth,  we seek to characterize the interior  equilibria, i.e. when $0<n_L<1$ and $0<n_F<1$.

Thus, we look for all NE by which  $0< x_n<1$ ($x_n$ characterized in the previous stage of the game). In this case, using \eqref{equ:payoffF}, \eqref{equ:payoffL}, and \eqref{equ:EUs_stage4}, the payoffs of SPs would be:
\begin{equation}\label{equ:payoffF_2}
\small
\pi_F=(\frac{I_F}{I_L}+p_L-p_F)(p_F-c)-sI^2_F
\end{equation}
\begin{equation}\label{equ:payoffL_2}
\small
\pi_L=(\frac{I_L-I_F}{I_L}+p_F-p_L)(p_L-c)+sI^2_F-\gamma I^2_L
\end{equation}
\normalsize2
In the following theorem, we prove that the NE uniquely exists and  characterize it:
\begin{theorem}\label{theorem:stage 3}
The NE strategies of Stage 3 by which  $0<n_L<1$ and $0<n_F<1$ are unique, and are $p_F=c+\frac{I_L+I_F}{3I_L}$ and $p_L=c+\frac{2I_L-I_F}{3I_L}$.
\end{theorem}

\begin{remark}
Note that $\frac{d p_L}{d I_L}\geq 0$ and $\frac{d p_F}{d I_L}\leq 0$. Thus, as we expect intuitively, $p_L$ (respectively, $p_F$) is increasing (respectively, decreasing) with respect to $I_L$. Also, $p_L$ (respectively, $p_F$) is decreasing (respectively, increasing) with respect to $I_F$.
\end{remark}

\begin{proof}
In this case, every NE should satisfy the first order condition. Thus, $p^*_F$ and $p^*_L$ should be determined such that $\frac{d \pi_F}{dp_F}=0$ and   $\frac{d \pi_L}{dp_L}=0$.  The first order conditions yield:
$$
\small
\begin{aligned}
2p^*_F-p^*_L&=\frac{I_F}{I_L}+c \qquad \& \qquad 2p^*_L-p^*_F&=\frac{I_L-I_F}{I_L}+c
\end{aligned}
$$
Thus,
\begin{eqnarray}\label{equ:optimump}
\small
\begin{aligned}
p^*_F=c+\frac{I_L+I_F}{3I_L} \qquad \& \qquad p^*_L=c+\frac{2I_L-I_F}{3I_L}
\end{aligned}
\end{eqnarray}
\normalsize
Therefore, $p^*_F$ and $p^*_L$ are the unique NE strategies if they yield $0< x_n< 1$ and no unilateral deviation is profitable for SPs. In Case A, we check the former, and in Case B, we check the latter.

\textbf{Case A:}
We first  check the condition that with $p^*_L$ and $p^*_F$, $0< x_n< 1$. Using \eqref{equ:xn} and \eqref{equ:optimump}:
$$
\small
x_n=\frac{I_L-I_F}{I_L}+p^*_F-p^*_L= \frac{2I_L-I_F}{3I_L}
$$
, which is greater than zero since $I_L\geq I_F$ and $I_L>0$. $x_n$ is also clearly less than one (note that $I_F\geq 0$). Thus, the condition $0< x_n< 1$ holds.

\textbf{Case B:}
Note that $\frac{d^2 \pi_F}{dp^2_F}<0$ and   $\frac{d^2 \pi_L}{dp^2_L}<0$. Thus, any solutions to the first order conditions would maximize the payoff of the SPs when $0<x_n<1$, and no unilateral deviation by which $0<x_n<1$ would be profitable for SPs.

Now, we discuss that any deviation by SPs such that $n_L=0$ and $n_L=1$ (which subsequently yields $n_F=1$ and $n_F=0$, respectively) is not profitable. Note that the payoff of SPs, \eqref{equ:payoffF} and \eqref{equ:payoffL}, are continuous  as $n_L\downarrow 0$, and $n_L\uparrow 1$ (which subsequently yields  $n_F\uparrow 1$ and $n_F\downarrow 0$, respectively). Thus, the payoffs of both SPs when selecting $p^*_L$ and $p^*_F$ (solutions of first order conditions) are greater than or equal to the payoffs when $n_L=0$ and $n_L=1$. Thus, any deviation by SPs such that $n_L=0$ or $n_L=1$ is not profitable for SPs.

This complete the proof that $p^*_F$ and $p^*_L$ are the unique NE strategies by which both SPs are active, i.e. $0<x_n<1$.
\end{proof}

\subsection*{Stage 2:}

In this stage of the game, SP$_F$ decides on the investment level, i.e. the number of resources to be leased from SP$_L$ ($I_F$), with the condition that $I_F\leq I_L$ to maximize $\pi_F$:
\begin{equation} \label{equ:maxpiFIF}
\max_{I_F\leq I_L}\pi_F=\max_{I_F\leq I_L}\big{(}\frac{I_L+I_F}{3I_L}\big{)}^2-sI^2_F
\end{equation}
Note that for the last equality, we used  \eqref{equ:payoffF_2} and Theorem~\ref{theorem:stage 3}.

\begin{theorem}
The optimum investment level of SP$_F$ is:
\begin{equation}\label{equ:optimumI_F}
I^*_F=\left\{
		\begin{array}{ll}
			\frac{I_L}{9sI^2_L-1} & \mbox{if } \quad I_L>\sqrt{\frac{2}{9s}}\\
			 I_L & \mbox{if } \quad  I_L\leq \sqrt{\frac{2}{9s}}
		\end{array}
	\right.
	\end{equation}
\end{theorem}

\begin{remark}
If the fee per resources ($s$) or the investment by SP$_L$ ($I_L$) is low, then SP$_F$ reserves all the available resources. If not, then SP$_F$ reserves a fraction of available resources ($I^*_F<I_L$). Note that in this case, $\frac{d I^*_F}{d I_L}<0$. Thus, the higher the number of available resources, the lower would be the number of resources reserved by SP$_F$.
\end{remark}

\begin{proof}
Note that in \eqref{equ:maxpiFIF}, $\pi_F$ is concave over $I_F$ if $s>\frac{1}{9I^2_L}$, $\pi_F$ is convex if $s<\frac{1}{9I^2_L}$, and is linear if $s=\frac{1}{9I^2_L}$.  We characterize the optimum investment level, i.e. $I^*_F$ under each of these conditions:

\subsection*{$s>\frac{1}{9I^2_L}$:}  In this case, the payoff of SP$_F$ is concave. Thus, the first order condition yields the possible optimum investment level, $\hat{I}_F$:
\begin{equation} \label{equ:hatI_f}
\begin{aligned}
\small
\frac{d \pi_F}{d I_F}|_{\hat{I}_F}=0&\Rightarrow \frac{2(\hat{I}_F+I_L)}{9I^2_L}-2s\hat{I}_F=0\\
&\Rightarrow \hat{I}_F=\frac{I_L}{9sI^2_L-1}
\end{aligned}
\end{equation}
Thus, the optimum investment level, $I^*_F$ is:
\begin{equation}
I^*_F=\min\{\hat{I}_F,I_L\}
\end{equation}
Note that from the assumption of this case $9sI^2_L>1$. Thus, $\hat{I}_F>0$ (using \eqref{equ:hatI_f}). In addition, from the expression of $\hat{I}_F$ \eqref{equ:hatI_f}:
$$
\hat{I}_F<I_L\iff 9sI^2_L>2
$$
Thus, the optimum investment level for SP$_F$ would be:
\begin{equation}
\small
I^*_F=\left\{
		\begin{array}{ll}
			\frac{I_L}{9sI^2_L-1} & \mbox{if } \quad I_L>\sqrt{\frac{2}{9s}}\\
			 I_L & \mbox{if } \quad  I_L\leq \sqrt{\frac{2}{9s}}
		\end{array}
	\right.
\end{equation}
\subsection*{$s<\frac{1}{9I^2_L}$:} In this case, $\pi_F$ would be convex. Thus, the optimum level of investment would be on the boundaries of the feasible set. Thus, $I^*_F$ would be either 0 or $I_L$, whichever yields a higher payoff. Note that:
$$
\begin{aligned}
\pi_F|_{I_F=0}&= \frac{1}{9}\qquad \& \qquad \pi_F|_{I_F=I_L}&=\frac{4}{9}-sI^2_L
\end{aligned}
$$
Thus, $I^*_F=I_L$ if and only if $\frac{4}{9}-sI^2_L\geq \frac{1}{9}$ (Note that we assumed that if $I_F=0$ and $I_F=I_L$ yield the same payoff, then SP$_F$ chooses $I_F=I_L$). Thus, $s\leq \frac{1}{3I^2_L}$ yields $I^*_F=I_L$. Note that from the assumption of the case,  $s<\frac{1}{9I^2_L}<\frac{1}{3I^2_L}$. Thus, $I^*_F=I_L$.

\subsection*{$s=\frac{1}{9I^2_L}$:} In this case, $\pi_F$ would be an increasing linear function of $I_F$. Thus, $I^*_F=I_L$.

Putting all the cases together, the result of the theorem follows.
\end{proof}

\subsection*{Stage 1:} In this stage, SP$_L$ decides on the level of investment $I_L$ to maximize her payoff, $\pi_L$:
\begin{equation}\label{equ:optimizationpil}
\max _{I_L}\pi_L=\max_{I_L} (\frac{2I_L-I^*_F}{3I_L})^2+sI^{*2}_F-\gamma I^2_L
\end{equation}
where for the last equality, we used  \eqref{equ:payoffL_2} and Theorem \ref{theorem:stage 3}, and $I^*_F$ is characterized in \eqref{equ:optimumI_F}.  In the next theorem, we characterize the candidate optimum answers:

\begin{theorem}\label{theorem:stage 1}
The optimum solution to \eqref{equ:optimizationpil} is $I^*_L=\min\{\sqrt{\frac{2}{9s}},\hat{I}_L\}$, where $\hat{I}_L$ is the solution of the first order condition on:
\begin{equation} \label{equ:payoff_1storder}
\pi_{L,I^*_F}=\frac{1}{9}(2-\frac{1}{9sI^2_L-1})^2+s\frac{I_L}{9sI^2_L-1}-\gamma I^2_L
\end{equation}
\end{theorem}

\begin{remark}
Theorem yields that the minimum optimum level of investment by SP$_L$ is $\sqrt{\frac{2}{9s}}$.
\end{remark}

\begin{remark}\label{remark:optimumIL}
Note that the first term in \eqref{equ:payoff_1storder} is increasing with $I_L$ (since the number of EUs is increasing with $I_L$). The second term is the payment from SP$_F$ which is decreasing with respect to $I_L$ (since when $I_L>\sqrt{\frac{2}{9s}}$, $I^*_F$ is decreasing with respect to $I_L$). The third term is the cost of investment which is decreasing with $I_L$. Thus, when either $s$ or $\gamma$ is sufficiently large, we expect the utility \eqref{equ:payoff_1storder} to be decreasing with respect to $I_L$, and the optimum answer to be $I^*_L=\sqrt{\frac{2}{9s}}$. On the other hand, $\hat{I}_L$ is expected to be the optimum answer when both $s$ and $\gamma$ are sufficiently small. We discuss about this in Numerical Results (Section~\ref{section:numerical}). Numerical results reveal that for large sets of parameters, we expect $I^*_F=\sqrt{\frac{2}{9s}}$.
\end{remark}

\begin{proof}
To find the optimum level of investment, $I^*_L$, consider two scenarios:
\subsection*{$I_L\leq \sqrt{\frac{2}{9s}}$:} In this case, $I^*_F=I_L$ by \eqref{equ:optimumI_F}. Thus, the maximization \eqref{equ:optimizationpil} would become:
$$
\max_{I_L\leq \sqrt{\frac{2}{9s}}}(s-\gamma) I^2_L
$$
Thus:
\begin{equation}
I^*_L=\left\{
		\begin{array}{ll}
			\sqrt{\frac{2}{9s}} & \mbox{if } \quad s\geq \gamma\\
			 0 & \mbox{if } \quad  s< \gamma
		\end{array}
	\right.
\end{equation}
Note that we assumed that if $I_L>0$ and $I_L=0$ yield the same payoff for  SP$_L$, then this SP chooses $I_L>0$. This is the reason that for $s=\gamma$, $I_L>0$ is chosen. Also, recall that to have a non-trivial problem, we assumed $s\geq \gamma$. Thus, if $I_L\leq \sqrt{\frac{2}{9s}}$, then $I^*_L=\sqrt{\frac{2}{9s}}$.

\subsection*{$I_L>\sqrt{\frac{2}{9s}}$:} In this case, $I^*_F=\frac{I_L}{9sI^2_L-1}$. Thus, the maximization \eqref{equ:optimizationpil} would become:
\begin{equation}
\small
\max_{I_L>\sqrt{\frac{2}{9s}}}{\pi_{L,I^*_F}=\frac{1}{9}(2-\frac{1}{9sI^2_L-1})^2+s\frac{I_L}{9sI^2_L-1}-\gamma I^2_L}
\end{equation}

Note that $I^*_F$ \eqref{equ:optimumI_F} is continuous at $I_L=\sqrt{\frac{2}{9s}}$. Thus, as $I_L\rightarrow \sqrt{\frac{2}{9s}}$, $\pi_{L,I^*_F}\rightarrow \pi_L|_{I_L=\sqrt{\frac{2}{9s}}}$, which is considered in the previous case.  Thus, the only possible optimum answer is the solution to the first order condition on $\pi_{L,I^*_F}$, i.e. $\hat{I}_L$ if $\hat{I}_L>\sqrt{\frac{2}{9s}}$. Results of the theorem follows.

\end{proof}

\section{The Outcome of the Game} \label{section:outcome}

In this section, we characterize the equilibrium outcome of the game using the results of the previous section and discuss about them. In Corollaries~\ref{corollary:outcome_IL>} and \ref{corollary:outcome_IL<}, we characterize
the equilibrium outcome when $I^*_L>\sqrt{\frac{2}{9s}}$  and $I^*_L\leq \sqrt{\frac{2}{9s}}$, respectively. Note that in each case, there exists a unique outcome and the two cases are mutually exclusive. Thus, the outcome of the game exists and is unique.

\subsection{If $I^*_L>\sqrt{\frac{2}{9s}}$:}
\begin{corollary} [Outcome A]\label{corollary:outcome_IL>}
If $I^*_L>\sqrt{\frac{2}{9s}}$, then:
\begin{itemize}
\item\textbf{Stage 1:}  The optimum level of investment of SP$_L$ is $I^*_L=\hat{I}_L$ (characterized in Theorem~\ref{theorem:stage 1}).
\item \textbf{Stage 2:} The optimum level of investment by SP$_F$ is: $I^*_F=\frac{I^*_L}{9sI^{* 2}_L-1}$.
\item \textbf{Stage 3:} Prices for EUs are: $p^*_F=c+\frac{1+\frac{1}{9sI^2_L-1}}{3}$ and $p^*_L=c+\frac{2-\frac{1}{9sI^2_L-1}}{3}$.
\item \textbf{Stage 4:} The fractions of EUs with each SP are: $n^*_F=\frac{1+\frac{1}{9sI^2_L-1}}{3}$ and $n^*_L=\frac{2-\frac{1}{9sI^2_L-1}}{3}$.
\end{itemize}
\end{corollary}
\begin{proof}
The outcome of Stage 2 directly follows from the results of Stage 2. The outcome of Stage 3 follows from \eqref{equ:optimump} and the outcome of Stage 2. The outcome of Stage 4 follows from  \eqref{equ:EUs_stage4} and the outcomes of Stage 2 and 3.
\end{proof}

Note that in this case, as discussed, $I^*_F$ is decreasing with respect to $I^*_L$ and $s$. In addition, given that $I_L$ is fixed, $p^*_F$ and $p^*_L$ are decreasing and increasing with respect to the per resource fee, $s$, respectively. Also, given that $s$ is fixed, $p^*_F$ and $p^*_L$ are decreasing and increasing with $I^*_L$, respectively.  Regardless of $p^*_F$ (respectively, $p^*_L$) being decreasing (respectively, increasing), the number of EUs with SP$_F$ (respectively, SP$_L$), i.e. $n^*_F$ (respectively, $n^*_L$) is still decreasing (respectively, increasing) with respect to $I^*_L$ (if $s$ fixed). This is because of the increase (respectively, decrease) in $t_F$ (respectively, $t_L$), i.e. the transport cost of SP$_F$ (respectively, SP$_L$). To understand these changes in the transport costs, recall that increasing $I^*_L$, decreases $I^*_F$, and subsequently increases $t_F$.

 In  Section~\ref{section:numerical}, we  observe that $I^*_L$ is dependent on $s$. Thus, the relationship between $s$ and the outcome is more complicated.

\subsection{If $I^*_L\leq \sqrt{\frac{2}{9s}}$:}

\begin{corollary}[Outcome B]\label{corollary:outcome_IL<}
If $I^*_L\leq \sqrt{\frac{2}{9s}}$, then:
\begin{itemize}
\item \textbf{Stage 1:} $I^*_L= \sqrt{\frac{2}{9s}}$.
\item \textbf{Stage 2:} $I^*_F=I^*_L$.
\item \textbf{Stage 3:} $p^*_F=c+\frac{2}{3}$ and $p^*_L=c+\frac{1}{3}$.
\item \textbf{Stage 4:}  $n^*_F=\frac{2}{3}$ and $n^*_L=\frac{1}{3}$.
\end{itemize}
\end{corollary}
Proof is similar to the proof of the previous corollary.

In this case,  SP$_F$ reserves all available resources, and  the investment level of SP$_L$ (which is equal to the number of resources reserved by SP$_F$) is a decreasing function of the fee per resources, i.e. $s$. SP$_F$ quotes a higher price for EUs in comparison to SP$_L$. In spite of the higher price, SP$_F$ would be able to attract more EUs given the better investment level in comparison to SP$_L$ which translates into a lower transport cost.

We  calculate $\pi_L$ and $\pi_F$ using Corollary~\ref{corollary:outcome_IL<} and \eqref{equ:payoffL}:
\begin{equation}
\begin{aligned}
\pi^*_L&=\frac{1}{3}-\frac{2 \gamma}{9 s}\qquad \& \qquad \pi^*_F&=\frac{2}{9}
\end{aligned}
\end{equation}
Thus, the payoff of SP$_L$ would be higher than the payoff of SP$_F$, i.e. $\pi^*_L>\pi^*_F$, if and only if $s>2\gamma$. In this case, although in comparison to SP$_F$, SP$_L$ offers her service to EUs with lower price and  attracts a lower number of EUs, she can still obtain a higher payoff through the per resource fee that she collects from SP$_F$. Thus, SP$_L$ leases the resources to SP$_F$ and instead generates revenue through the fee she charges to SP$_F$.

\section{Benchmark Case} \label{section:benchmark}

Now, we consider the case in which SP$_L$ does not invest in her infrastructure. Similar to the original model, both SP$_L$ and SP$_F$ compete to attract a pool of undecided users. The hotelling model for EUs is the same as before. The only difference is the transport cost. Here, we consider the transport costs to be parameters, and denote them by $t_L>0$ and $t_F>0$ for SP$_L$ and SP$_F$, respectively. We discuss about the effects of these parameters in the Numerical Results (Section~\ref{section:numerical}).

The analysis of stage 4 of the game would be the same as before, and by \eqref{equ:xn_t_parameter}:

\begin{equation}\label{equ:xn_banchmark}
\small
x_n=\frac{t_F+p_F-p_L}{t_L+t_F}
\end{equation}
Thus,
\begin{equation} \label{equ:EUs_stage4_benchmark}
\small
\begin{aligned}
n_L &=
\left\{
\begin{array}{ll}
	0  & \mbox{if } x_n < 0 \\
	\frac{t_F+p_F-p_L}{t_L+t_F} & \mbox{if } 0\leq x_n \leq 1\\
	1 & \mbox{if } x_n>1
\end{array}
\right. \\
n_F&=1-n_L
\end{aligned}
\end{equation}

Note that  in the benchmark case, there is no investment decision. Thus, the payoffs of SPs would be modified as:

\begin{equation} \label{equ:payoff_F_bench}
\small
\pi_F=n_F(p_F-c)
\end{equation}
\begin{equation} \label{equ:payoff_L_bench}
\small
\pi_L=n_L(p_L-c)
\end{equation}

Thus, the only decision of SPs is to determine their pricing for EUs. We use the same approach as Stage 3 of the game to characterize the NE pricing strategies for SPs by which $0<x_n<1$:

\begin{theorem}
The NE pricing strategies of SPs are:
\begin{eqnarray}
\small
\begin{aligned}
p^*_F&=c+\frac{2t_L+t_F}{3}\\
p^*_L&=c+\frac{t_L+2t_F}{3}
\nonumber
\end{aligned}
\end{eqnarray}
\end{theorem}
\begin{proof}
Proof is similar to the proof of Theorem~\ref{theorem:stage 1}.  $p^*_F$ and $p^*_L$ should be determined such that $\frac{d \pi_F}{dp_F}=0$ and   $\frac{d \pi_L}{dp_L}=0$:
\begin{eqnarray}\label{equ:optimump_benchmark}
\small
\begin{aligned}
p^*_F&=c+\frac{2t_L+t_F}{3}\qquad \& \qquad p^*_L&=c+\frac{2t_F+t_L}{3}
\end{aligned}
\end{eqnarray}

Therefore, $p^*_F$ and $p^*_L$ are the unique NE strategies if they yield $0< x_n< 1$ and no unilateral deviation is profitable for SPs. In Case A, we check the former, and in Case B, we check the latter.

\textbf{Case A:}
We first  check the condition that with $p^*_L$ and $p^*_F$, $0< x_n< 1$. Using \eqref{equ:xn_banchmark} and \eqref{equ:optimump_benchmark}:
$$
x_n=\frac{2t_F+t_L}{3(t_F+t_L)}
$$
, which is greater than zero since $t_F>0$ and $t_L>0$. $x_n$ is also clearly less than one. Thus, the condition $0< x_n< 1$ holds.

\textbf{Case B:}
Note that $\frac{d^2 \pi_F}{dp^2_F}<0$ and   $\frac{d^2 \pi_L}{dp^2_L}<0$. Thus, any solutions to the first order conditions would maximize the payoff of the SPs when $0<x_n<1$, and no unilateral deviation by which $0<x_n<1$ would be profitable for SPs.

Now, we discuss that any deviation by SPs such that $n_L=0$ and $n_L=1$ (which subsequently yields $n_F=1$ and $n_F=0$, respectively) is not profitable. Note that the payoff of SPs, \eqref{equ:payoff_F_bench} and \eqref{equ:payoff_L_bench}, are continuous  as $n_L\downarrow 0$, and $n_L\uparrow 1$ (which subsequently yields  $n_F\uparrow 1$ and $n_F\downarrow 0$, respectively). Thus, the payoffs of both SPs when selecting $p^*_L$ and $p^*_F$ (solutions of first order conditions) are greater than or equal to the payoffs when $n_L=0$ and $n_L=1$. Thus, any deviation by SPs such that $n_L=0$ or $n_L=1$ is not profitable for SPs.

This complete the proof that $p^*_F$ and $p^*_L$ are the unique NE strategies by which both SPs are active, i.e. $0<x_n<1$.
\end{proof}

Note that in the absence of  investments, Stage 2 and 1 would be of no importance. In the following corollary, we characterize the outcome of the game in the benchmark case (subscript B stands for Benchmark):
\begin{corollary}
The equilibrium outcome of the benchmark case is as follows:
\begin{itemize}
\item NE pricing strategies are: $p^*_{F,B}=c+\frac{2t_L+t_F}{3}$ and $p^*_{L,B}=c+\frac{t_L+2t_F}{3}$.
\item Fractions of EUs with each SP: $n^*_{L,B}=\frac{2t_F+t_L}{3(t_F+t_L)}$ and $n^*_{F,B}=\frac{t_F+2t_L}{3(t_F+t_L)}$.
\end{itemize}
\end{corollary}

Using \eqref{equ:payoff_F_bench} and \eqref{equ:payoff_L_bench}, the payoffs of SPs in these cases are:
\begin{equation}
\begin{aligned}
\pi^*_{F,B}&=(\frac{t_F+2t_L}{3(t_F+t_L)})^2\qquad \& \qquad \pi^*_{L,B}&=(\frac{2t_F+t_L}{3(t_F+t_L)})^2
\end{aligned}
\end{equation}

Note that if $t_L>t_F$, then $p^*_{F,B}>p^*_{L,B}$, $n^*_{F,B}>n^*_{L,B}$, and subsequently $\pi^*_{F,B}>\pi^*_{L,B}$, and vice versa. Thus, the SP that EUs have lower reluctance for, receives the highest payoff.

\section{Numerical Results  and Discussions}  \label{section:numerical}

In this section, we use numerical simulations (i) to determine whether and under what conditions the outcome in Corollary~\ref{corollary:outcome_IL>} (Outcome A) would emerge, and (ii)  to provide insights for results under different parameters of the model. For all results, we consider $c=1$.\footnote{Note that the choice of $c$ barely affects the results. It may only shift some of the results (e.g. the price charged to EUs) by only a constant.}

\begin{figure}[t]
	\begin{subfigure}{.25\textwidth}
		\centering
		\includegraphics[width=\linewidth]{resources_gamma01_2}

		\label{figure:resources_gamma0.1}
	\end{subfigure}%
	\begin{subfigure}{.25\textwidth}
		\centering
		\includegraphics[width=\linewidth]{resources_gamma0}

		\label{fig:resources_gamma0}
	\end{subfigure}
	\caption{Investment decisions of SPs vs. per resource fee, $s$}\label{figure:resources}
\end{figure}

In Figure~\ref{figure:resources}, we plot the optimum level of investment of SP$_L$ ($I^*_L$), the number of resources that SP$_F$ reserves ($I^*_F$), and the minimum optimum level of investment by SP$_L$ ($\sqrt{\frac{2}{9s}}$), when $\gamma=0.1$ (left) and $\gamma=0$ (right) (recall that,  in \eqref{equ:payoffL}, $\gamma$ is the marginal cost of investment.). The discontinuities  are because of the transition of the outcome of the game from the Outcome A to Outcome B (outcome in Corollary~\ref{corollary:outcome_IL<}). Results reveal that for each outcome, $I^*_L$ and $I^*_F$ are decreasing with the per resource fee, $s$. Thus, the higher the fee per resources that SP$_F$ pays to SP$_L$, the lower would be the number of resources that SP$_F$ reserves, and subsequently the lower  would be the investment of SP$_L$.  Also, when the marginal cost of investment ($\gamma$) is zero, SP$_L$ investment is significantly more in comparison to $\gamma=0.1$. Thus, intuitively, the higher the investment cost, the lower would be the investment.

Results also confirm the intuitions presented in Remark~\ref{remark:optimumIL} that for sufficiently small $s$ and $\gamma$, the optimum level of investment by SP$_L$ would be equal to $\hat{I}_L$ (introduced in Theorem~\ref{theorem:stage 1}) which is  higher than $\sqrt{\frac{2}{9s}}$ (Outcome A). Numerical results also reveal that for $\gamma>0.12$, Outcome A would not occur even for small $s$. Thus, the outcome in which $\hat{I}_L$ is the optimum level of investment only occurs for a small set of parameters. Therefore, for a wide range of parameters, we expect Outcome B to be the outcome of the game.

In Figure~\ref{figure:prices}, we plot the pricing decisions of SPs for EUs, i.e. $p^*_L$ and $p^*_F$, for $\gamma=0.1$ (left) and $\gamma=0$ (right). Similar to Figure~\ref{figure:resources},  the discontinuities in the figures are because of the transition of the outcome of the game from Outcome A to B.  Note that in Outcome B (when $s$ is sufficiently large), $p^*_F$ and $p^*_L$ are constant (given that c is a constant) independent of $\gamma$ and $s$, and the price that SP$_F$ charges is twice the price that SP$_L$ charges. However, in Region A, if $\gamma$, i.e. the marginal cost of investment, is extremely small (e.g. zero), then SP$_L$ would be able to charge a higher price than SP$_F$ (Figure~\ref{figure:prices}-right). The reason is that for $\gamma$ small, SP$_L$ will invest more ($I^*_L$ small). We also stated in Section~\ref{section:outcome} after Corollary~\ref{corollary:outcome_IL>} that  $p^*_F$  and $p^*_L$ are decreasing and increasing with $I^*_L$, respectively. Thus, small investment cost yields a higher investment by SP$_L$ and as a result a higher price for EUs of this SP.

\begin{figure}[t]
	\begin{subfigure}{.25\textwidth}
		\centering
		\includegraphics[width=\linewidth]{prices_gamma01_2}

		\label{fig:prices_gamma01}
	\end{subfigure}%
	\begin{subfigure}{.25\textwidth}
		\centering
		\includegraphics[width=\linewidth]{prices_gamma0}

		\label{fig:prices_gamma0}
	\end{subfigure}
	\caption{Pricing decisions of SPs for EUs versus $s$}\label{figure:prices}
\end{figure}

Also note that  in Outcome A,  depending on $\gamma$, prices can be increasing or decreasing with $s$: If the marginal cost of investment is extremely small, then $p^*_F$ is increasing with $s$ and $p^*_L$ is decreasing with $s$ (Figure~\ref{figure:prices}-right). The opposite is true when $\gamma$ is not small  (Figure~\ref{figure:prices}-left). To describe the reason behind this result, note that from Corollary~\ref{corollary:outcome_IL>}, $p^*_L$ and $p^*_F$ are increasing and decreasing, respectively, with respect to both $s$ and $I^*_L$. On the other hand, $I^*_L$ is itself decreasing with $s$ (Figure~\ref{figure:resources}). Thus, as $s$ increases two factors one increasing (s) and one decreasing ($I^*_L$) affect the prices. Numerical results yield that when $\gamma$ is extremely small (thus, $I^*_L$ is large), the rate of change in $I^*_L$ with $s$ dominates the rate of change in $s$. Thus, as $s$ increases, $p^*_F$ increases and $p^*_L$ decreases.

Note that by Corollaries \ref{corollary:outcome_IL>} and \ref{corollary:outcome_IL<}, the dependency of $n^*_L$ and $n^*_F$ to parameters of the model is  similar to the dependency of $p^*_L$ and $p^*_F$. The only difference is the exclusion of $c$ from the expressions.

In Figure~\ref{figure:payoffs},  we plot the payoffs of SPs, i.e. $\pi^*_F$ and $\pi^*_L$, and the payoff of SP$_L$ in the benchmark case for three scenarios (i) $t_F=0.5$ and $t_L=0.5$, (ii) $t_F=1$ and $t_L=0$, and (iii) $t_F=0$ and $t_L=1$. Results reveal that although for some parameters, the price that SP$_F$ quotes for EUs and the number of EUs that she can attract is higher than the price  of SP$_L$  and the number of EUs of this SP, SP$_L$ can still fetch a higher payoff than SP$_F$. The reason is the fee that SP$_F$ pays to SP$_L$.

\begin{figure}[t]
	\centering
	\includegraphics[width=0.4\textwidth]{payoffs_gamma01_B}
	\caption{Payoffs of SPs for EUs versus $s$}
	\label{figure:payoffs}
\end{figure}

In addition, we  can only expect an investment by SP$_L$ if this investment yields a higher payoff than the benchmark case with no investment. Results in Figure~\ref{figure:payoffs} yield that the higher the transport cost of SP$_L$ in comparison to SP$_F$, the lower would be the payoff of SP$_L$ in the benchmark case, and the higher would be the incentive of this SP to invest and lease her resources to SP$_F$. Thus, the preferences of EUs for SPs, i.e. the transport costs, are one of the main factors in gauging the plausibility of a cooperation between an MNO and an MVNO. If an MNO knows that EUs are reluctant joining her, i.e. EUs have a high transport costs for the MNO, then she may invest in her infrastructure, lease a portion of the new resources to an independent MVNO, which EUs have lower reluctance for, and collect the profit through the fees she charges to the MVNO.

In Figure~\ref{figure:payoffs}, note that if $s$ is small, then SP$_F$ fetches a higher payoff than SP$_L$ since she can get a higher revenue from EUs (Figure~\ref{figure:prices}-left and recalling that $n^*_F$ has the same pattern as $p^*_F$). Overall, increasing $s$ increases the payoff of SP$_L$. Counter intuitively, when $s$ is large enough such that Outcome B occurs (Corollary~\ref{corollary:outcome_IL<}) SP$_F$ receives a higher payoff in comparison to her payoffs when $s$ is small. This occurs since in this case, SP$_F$ reserves all the available resources from SP$_L$. This  enables SP$_F$ to charge a high price for EUs while still attracting EUs. Thus, when $s$ is large enough, both SPs receive high revenue at the expense of EUs paying more. Thus, the welfare of EUs would be low. This highlights the importance  of a framework by which the per resource fee is determined such that all entities of the market can benefit. This framework could be a bargaining game between SPs with some restrictions imposed by a regulator. This is a topic of future work.

\bibliographystyle{IEEEtran}