\section{Introduction}
In the past decade, physics of  polaritons, a quantum superposition of excitons and cavity photons,  attracts significant
attention due to fundamental interest to macroscopic collective quantum phenomena in this system \cite{Carusotto:13} as well as
due to potential applications in photonics and optoelectronics \cite{Gibbs:11}.
Thanks to their small effective mass, polaritons
demonstrate the  Bose-Einstein condensation (BEC) and superfluidity at temperatures about six orders higher than that for
ultra-cold atomic gases \cite{Deng:02,Kasprzak:06,Balili:07,Utsunomiya:08,Amo:09a}.
 Numerous applications of
polariton condensates  in optical computing \cite{Bajoni:08,Menon:10,Gao:12,Ballarini:13,Berman:14}, nonlinear
interferometry \cite{Sturm:14}, novel light sources \cite{Deng:02a,Malpuech:02a,Szymanska:02,Bajoni:07,Christopoulos:07},
and  atomtronics \cite{Restrepo:14}
have recently been proposed.

Transition-metal dichalcogenide monolayer crystals  provide unprecedented opportunity to transfer
 these achievements to the room temperature scale.
In contrast to another near-perfect two-dimensional material,  graphene,
single-layer molybdenum and tundsen-based transition metal dichalcogenides
are direct-zone semiconductors with the sizable bandgap
$E_g\sim 1-2$ eV \cite{Das:14}.
Extremely  large light-matter interaction strength (the Rabi splitting $\sim 50$ meV)  \cite{Liu:15} and  large
exciton binding energy $\sim 0.3-0.9$ eV  in monolayer
transition metal dichalcogenides \cite{Qiu:13,Chernikov:14} make these materials exceptional candidates
 for observation of room temperature quantum collective phenomena.
Observation of polaritons in a microcavity with embedded molybdenum disulfide (MoS$_2$)
monolayer at the room temperature has recently been reported in Ref.~\cite{Liu:15}.

However, it is  {\it a priory} unclear  if the conditions for a superfluid transition can be satisfied in actual experiments
for a polariton gas in cavities with transition-metal dichalcogenide monolayers at room temperatures.
One of the important factors is significant temperature
fluctuations in the system
that results in spreading of a polariton cloud in a trap and thus, in a drop of the polariton density.
We note that temperature spreading  is not essential for conventional low-temperature polariton experiments
with semiconductor microcavities where the leading factors are quantum fluctuations
and polariton-polariton interactions ~\cite{Berman:08}.

 The goal of this work is to determine a range of experimentally relevant parameters,
for which superfluidity of a polariton gas can be achieved at room temperatures
in an optical microcavity with an embedded MoS$_2$ monolayer.
Specifically, by considering the quasiclassical diffusive dynamics of a trapped
polariton gas
we determined the values of the trapping potential strength
and laser pump power, for which the room temperature
polariton superfluidity can be observed.

 \begin{figure}[t]
 \includegraphics[width=8cm]{tc.eps}
 \caption{\label{fig:tc}
Temperature of the Kosterlitz-Thouless (KT) phase transition, $T_c$, as a function of the polariton gas density $n$
estimated from Eq. (\ref{eq:tc}) where the cutoff length $r_0$ is set equal to the exciton Bohr radius
in molybdenum disulfide $a_B \approx 1$ nm and the polariton mass is $m = 3.14 \times 10^{-5} m_0$, see text.
It is seen that the KT transition at the room temperature $T=300$ K occurs at the density
$n \approx 4.2 \times 10^{12}$ m$^{-2}$.
}
\end{figure}

\section{Superfluid phase transition in a  polariton gas}
A homogenous two-dimensional Bose gas can only transits to a Bose condensed state, in which the occupation
of the lowest energetic level is macroscopically large, at zero temperature since the temperature fluctuations
destroy the long-range order in the system.
In contrast, an ideal trapped quantum Bose gas undergoes the BEC transition at a finite temperature
$T_0 = \hbar (\pi k_B)^{-1}\sqrt{6k N/m}$ where $k$ is the effective trapping potential strength, $N$ is the total number of particles, and
$m$ is the particle mass \cite{Bagnato:91}.
At low temperatures $T \ll T_0$ virtually all particles occupy the lowest energetic level that is,
a condensate is formed \cite{Dalfovo:99}.
However, for a Bose gas of interacting particles, the condensate can be significantly depleted due to mutual particle scattering:
the condensate fraction for a system of $N \gg 1$ particles falls as  $N_c/N \propto \log(N)/N^{1/2}$ \cite{Holzmann:07}.
In effect, despite $T_0$ for an ideal gas of polaritons can be estimated as a few hundred K because of a small mass $m$,
the actual BEC transition temperature of a trapped polariton gas lies in the helium temperature range~\cite{Berman:08}.

Nevertheless, a Bose system with interactions can transit to a superfluid state even if the condensate fraction is small or  absent et all.
The famous example is superfluid $^4$He (He-II), in which the occupation of the lowest energetic
level is only $\sim 10$\%, and $\sim 90$\% of particles are ``pushed'' out of the condensate to higher energetic levels
due to strong interatomic interactions \cite{Atkins:59}. The relation between the Bose condensation and superfluidity is
delicate and is discussed in a number of
publications (e.g., \cite{Dalfovo:99,Yukalov:11,Berman:11a}).
The superfluid transition in a two-dimensional system with interactions is a Kosterlitz-Thouless (KT)  phase transition, which is related to pairing of
quantized vortices \cite{Kosterlitz:73}. The KT transition temperature for a  gas of interacting polaritons
is given by the following equation~\cite{Lozovik:06}
\begin{equation}
T_c = \frac{2 \pi \hbar^2 n}{m \log \log (1/n r_0^2)} \label{eq:tc}
\end{equation}
where $n$ is the gas density and $r_0$ is the small-scale cutoff length. The  dependence (\ref{eq:tc}) for polaritons in a cavity with embedded
molybdenum disulfide is shown in Fig.\ \ref{fig:tc}. It is seen in Fig.\ \ref{fig:tc} that a polariton
gas  can undergo a superfluid transition at the room temperature $T=300$ K if the gas density is
high enough, $n>n_c \approx 4.2 \times 10^{12}$ m$^{-2}$. At $n\sim n_c$, the
thermal wavelength of polaritons  $\lambda_T = 2^{1/2}  \pi \hbar / (m k_B T)^{1/2} \approx 1.3$ $\mu$m
is about five times larger than the inter-particle distance $r_s = (\pi n)^{-1/2} \approx 0.28$ $\mu$m and thus,
quantum collective effects are important at densities above $n_c$. However, the gas parameter is small,
$n_c a_B^2 \sim 4 \times 10^{-6}$,
that justifies the application of the weakly interacting gas  model (\ref{eq:diffusion}), see the next sections.

In what follows we identify the range of experimental parameters at which the KT superfluid transition can be observed
in a trapped polariton gas in microcavities with an embedded transition-metal dichalcogenide monolayer.
In our consideration, we focus on microcavities with an embedded molybdenum disulfide
monolayer that is, on a system where room-temperature polaritons have  recently been
 observed \cite{Liu:15}.

\section{Simulation Method}
The dynamics of a normal  polariton gas at the scales much larger than the polariton mean free path can be
described by the driven diffusion equation
\begin{equation}
{\partial n\over \partial t} = - \mu \, {\rm div} (\bm{F} n) +  D \Delta n - \Gamma n + P_n
\label{eq:diffusion}
\end{equation}
where $n\equiv n(\bm{r},t)$ is the two-dimensional density of the polariton gas, which depends on the
spatial coordinate $\bm{r}=(x,y)$ in the microcavity plane and time $t$, $\mu$ is the polariton mobility,
$\bm{F}\equiv \bm{F}(\bm{r},t)$ is the external force acting on polaritons, $D$ is the diffusivity,
$\Gamma$ is the width of the polariton resonance
and $P_n$ is the source (see Ref. \cite{Carusotto:13} for extensive review). In what follows, we
consider the dynamics of the polariton gas in an in-plane trapping potential $U_{trap}(\bm{r}) = {1 \over 2} k r^2$, where
$k$ is the effective strength of the trap and $r = |\bm{r}|$. The main mechanisms of the diffusive collective motion of
polaritons include their mutual interactions, scattering on phonons and scettring on excitonic and photonic defects
in the microcavity \cite{Bley:98}.  (We consider an undoped cavity.)
In our simulations, we take into account the density-dependent polariton interactions in the mean-field approximation
where the effective repulsive interaction potential is $U_{pp} (\bm{r},t)= g n(\bm{r},t)$, which describes exactly the same
interactions that lead to the blue shift of the polariton excitation spectrum in a superfluid state  \cite{Carusotto:13}.
Here, $g=6 E_B a_B^2 |X_{k}|^4$ is the effective polariton-polariton interaction strength \cite{Carusotto:13},
with $E_B$ and $a_B$ to be the respective binding energy and Bohr radius of the excitons, while $X_{k}$ is the
Hopfield coefficient that describes the excitonic fraction in the polariton wave function.

In our approach, the driven-diffusive dynamics of the polariton gas is captured via a quasiclassic
stochastic differential equation for the  center-of-mass coordinate $\bm{r}(t)$ of the polariton wave packets,
\begin{equation}
 d\bm{r}(t) = \mu \bm{F}(\bm{r}(t), t) dt + \sqrt{2D}d\bm{W}(t). \label{eq:stochastic}
\end{equation}
The first term in Eq.\ (\ref{eq:stochastic}) describes the effect of the  force $\bm{F}(\bm{r},t)$ on the
polariton diffusion and the second term is the Brownian contribution with $d\bm{W}(t)$ to be
the differential of a Wiener process with  unit variance.  The ensemble-averaged
distribution of ``tracer'' particles  (or quasiparticles)  that obey Eq.\ (\ref{eq:stochastic}) reproduces the  driven-diffusive dynamics (\ref{eq:diffusion}) \cite{Szymczak03}.
The advantage of the model  (\ref{eq:stochastic})  is that, in contrast to  Eq.\ (\ref{eq:diffusion}),
the stochastic equation describes the fluctuations of the particle density.
The model  (\ref{eq:stochastic}) is off lattice and is computationally effective for systems with a
large number of particles.
The stochastic model  (\ref{eq:stochastic})
was applied for the description of diffusive properties of quantum gases  (see Ref. \cite{Carusotto:13})
as well as for nanoparticle diffusion in fluids \cite{Szymczak03}.

We numerically integrated Eq.\ (\ref{eq:stochastic})  by using  the  Euler method.
The force acting on polaritons was set equal to $\bm{F} = - \nabla U$ where $U= U_{trap} + U_{pp}$.
The finite-size polariton source at the pumping laser light spot was modeled by adding  particles at each numerical
time step $\Delta t$
with the Gaussian spatial probability distribution centered at the spot center.
 The finite lifetime of polaritons was taken into account by randomly eliminating
particles from the ensemble at each   time step with the probability $w = \Gamma \Delta t$.

\section{Model parameters}
In our studies, we use the parameters relevant for the A-exciton polaritons
in a cavity with an embedded
molybdenum disulfide monolayer  \cite{Liu:15}. In our simulations, we
focus on the polariton dynamics at the room temperature $T=300$ K.

 The effective polariton mass was calculated as
\begin{equation}
m=2/(m_{ex}^{-1} + m_{ph}^{-1}) .
\end{equation}
The exciton mass in MoS$_2$ was taken equal to $m_{ex}=m_e + m_h = 0.78 m_0$
where $m_e=0.35 m_0$ and $m_h = 0.43 m_0$ are the electron and hole mass respectively
(with $m_0$ to be the free electron mass) \cite{Ramasubramaniam:12,Cheiwchanchamnangij:12}.
The effective photon mass is
$m_{ph}={\pi \sqrt{\epsilon} \hbar / c L_C}$,
where the cavity length $L_C$
was found from the condition of the resonance of the photonic and excitonic states,
$E_{ph} = {\pi \hbar c / \sqrt{\epsilon} L_C}$, where $E_{ph}$
is the energy of the cavity photon resonance. The latter was taken equal to $E_{ph}=1.87$ eV \cite{Liu:15}.
The dielectric constant of MoS$_2$ is $\epsilon = 4.3$ \cite{Ramasubramaniam:12}.

\begin{table}[tb]
\centering
\caption{\bf Simulations parameters for a cavity with embedded MoS$_2$ layer}
\begin{tabular}{|lll|}
\hline
Exciton mass & $m_{ex}$ & $0.78 m_0$ \\
Photon mass & $m_{ph}$ & $1.57 \times 10^{-5} m_0$ \\
Polariton mass & $m$ & $3.14 \times 10^{-5} m_0$ \\
Polariton momentum & $\tau_{pol} $ & $9.60 \times 10^{-13}$ s \vspace{-0.1cm}\\
\hspace{0.15cm} relaxation time & & \\
Polariton diffusion  & $D$ &  $139$ m$^2$/s \vspace{-0.1cm}\\
\hspace{0.15cm} coefficient & & \\
Polariton interaction & $g$ & $1.44 \times 10^{-3}$   \vspace{-0.1cm}\\
\hspace{0.15cm}  strength & &\hspace{0.15cm}  meV$\mu$m$^2$  \\
Polariton mobility & $\mu$ &  $3.36 \times 10^{22}$  \vspace{-0.1cm}\\
 & & \hspace{0.15cm} s/kg \\
Dielectric constant & $\epsilon$ & 4.3 \\
Polariton life time & $\tau_{life}^{(pol)}$ & 84 ps \\
Trapping potential  & $k$ & $15-300$ \vspace{-0.1cm}\\
\hspace{0.15cm} strength & &\hspace{0.15cm}  eV/cm$^2$ \\
Efficiency of polariton & $q $& $10^{-3}$ \vspace{-0.1cm} \\
\hspace{0.15cm} generation & & \\
Numerical unit of  & $\Delta x$ & 1 $\mu$m \vspace{-0.1cm}\\
\hspace{0.15cm}length & & \\
Numerical time step & $\Delta t =$ & 272 fs \vspace{-0.1cm}\\
& \hspace{0.15cm}${m  \Delta x^2}\hbar ^{-1}$ & \\
\hline
\end{tabular}
\label{tab:params}
\end{table}

The exciton diffusion coefficient was taken equal to  $D_{ex} = 14$ cm$^2$/s \cite{Ceballos:15}.
This value is comparable with the exciton diffusion coefficient in semiconductor quantum wells  \cite{Carusotto:13}.
The polariton  diffusivity was calculated as \cite{Bley:98,Berman:10b}
\begin{equation}
D=|X|^{-4} {m_{ex} \over m} D_{ex}.
\end{equation}
In what follows, we consider zero detuning between the photonic and excitonic resonances thus, we set  $X=1/\sqrt{2}$.
However, the detuning in the experiments with polaritons in a cavity with MoS$_2$ monolayer can be varied in a wide range up to
$\sim 50$ meV \cite{Liu:15}. Thus, the diffusive properties of polaritons can be controlled by
adjusting the detuning value.

The polariton mobility was calculated  as follows \cite{Bley:98,Berman:10b},
\begin{equation}
\mu={\tau_{pol} \over m}.
\end{equation}
The respective momentum relaxation times for excitons  and polaritons  was estimated as
$\tau_{ex}={D_{ex} m_{ex} / k_B T}$ and $\tau_{pol} = |X|^{-4} \tau_{ex}$.

The polariton resonance width was estimated as
\begin{equation}
\Gamma= |X|^2 \Gamma_{ex} + |C|^2  \Gamma_{ph} \label{eq:gamma}
\end{equation}
where $|C|^2 = 1-|X|^2$ is the photon Hopfield coefficient.
According to Ref.\ \cite{Palummo:15}, the radiative exciton lifetime $\tau_{life}^{(ex)} \equiv \Gamma_{ex}^{-1}$  in MoS$_2$ monolayer
increases from $\tau_{life}^{(ex)} = 0.23$ ps at $T\rightarrow 0$  to $\tau_{life}^{(ex)} \approx 4$ ps at  $T=10$ K and then, to
270 ps at the room temperature. Even larger radiative exciton life time of 850 ps was reported for MoS$_2$ monolayer crystals at the
room temperature \cite{Shi:13}.
On the other hand, large broadening of the  excitonic resonance $\hbar \Gamma_{ex} = 30$ meV
(with the respective time scale of 22 fs) was reported in Ref.\ \cite{Liu:15}.  The respective
short-time processes are probably related to non-radiative  trapping of excitons at impurities or  hot
carrier thermalization, in agreement with Ref.\ \cite{Palummo:15}.
In our simulation, we set $\tau_{life}^{(ex)} = 270$ ps.
The photon resonance width was taken equal to $\Gamma_{ph} = (50\, {\rm ps})^{-1}$ \cite{Nelsen:13}.
 The polariton life time was estimated  from Eq.\ (\ref{eq:gamma}) as $\Gamma^{-1} = \tau_{life}^{(pol)} =  84$~ps,
 which is in agreement with the observations for semiconductor-based microcavities \cite{Nelsen:13}.

In the simulations, we varied the trapping potential strength from $k=15$ eV/cm$^2$ to 300  eV/cm$^2$ that lies in the experimentally relevant parameter range \cite{Negoita:99,Berman:08}.
 The  trapping potential can be formed  by applying inhomogeneous
stress \cite{Negoita:99},  static electric or magnetic field, or laser radiation
(see Ref.\ \cite{High:12} and references therein).
 The polariton  interaction constant was calculated as $g=6  E_B a_B^2 |X|^4 = 1.44 \times 10^{-3}$ meV $\mu$m$^2$
 given the exciton binding energy $E_B = 0.96$ eV and the exciton Bohr radius $a_B \approx 1$ nm
 \cite{Liu:15,Shi:13}. It is worth noting that, despite the exciton binding energy $E_B$ in transition metal dichalcogenides
 is significantly larger than
 that in quantum wells, the polariton interaction strength is of the same order in the both cases.

 The polariton density at a point $\bm{r}$ at the moment $t$
was calculated from the positions of polariton
quasiparticles $\bm{r}_i (t)$ (where the index $i$ labels the quasiparticles) as follows
\begin{equation}
n(\bm{r},t) = \sum_i f(|\bm{r}_i (t) - \bm{r}|)
\end{equation}
with the weight function $f(r) =  (2\pi r_{av}^2)^{-1} \exp(-r^2 / 2 r_{av}^2)$ and the
averaging length $r_{av} = 10$ $\mu$m.

Qualitative analysis of experimental data in
Refs.~\cite{Christopoulos:07,Negoita:99,Feldmann:87,Butov:02a,Malpuech:02a} for quantum well structures
shows that the efficiency of the exciton
generation by a continuous wave (cw) light source is $q\sim 0.1$\%; in other words, only 1 of $10^3$ incident photons creates
a bound exciton, which is, in its turn, forms a polariton.
A single-layer transition metal dichalcogenide monolayer can capture up to $10$\% of incident photons \cite{Mak:10,Splendiani:10}.
However, the efficiency of the  polariton generation for cw pumping in a cavity with MoS$_2$ is not known; we suggest that
the value is close to that obtained for the quantum well structures.
Thus, in our model,
the polariton injection rate, $P_n$, was calculated as
\begin{equation}
P_n={q  P  \over E_{ph}}, \label{eq:pn}
\end{equation}
where $P$ is the pumping laser power.
In the simulations, we consider a Gaussian laser spot profile with full width at half maximum (FWHM) of 30
$\mu$m.

The simulations parameters for $T=300$ K are summarized in Table \ref{tab:params}.

\section{Results and Discussion}

First,  we explore the case where the excitation laser spot is positioned at the center of the trap.
The relaxation of the total number of polaritons in the system to its steady-state value
is shown in Fig.\ \ref{fig:np}. It is seen that the system reaches a steady state  $\sim 400$ ps
after the excitation is switched on.

 \begin{figure}[h]
 \includegraphics[width=8cm]{np.eps}\vspace{-0.5cm}
 \caption{\label{fig:np}
Relaxation to a steady state of a trapped polariton cloud after the laser pump is turned on at the moment $t=0$.
The data points show the total number of polaritons  as a function of time $t$.
The excitation power of the laser pump is $P=550$ mW, the trapping potential
strength is  $k=50$ eV/cm$^2$.
The fluctuations of the total number of polaritons
in the steady state $t \geq 400$ ps are seen in the figure.}
\end{figure}

The steady-state spatial distribution of polaritons in the trap for the trapping potential strength $k=50$ eV/cm$^2$ is shown in Fig.\ \ref{fig:profile}.
It is seen that the characteristic polariton cloud size reaches $\sim 200$ $\mu$m that is, about an order larger than the excitation
spot size 30 $\mu$m.
It is also much larger than the Thomas-Fermi size of a cloud at zero temperature
$a_{TF}=\sqrt{2 n_{max} g / k}\approx 1-10$ $\mu$m for the maximum density at the center
$n_{max} \approx n_c$ \cite{Berman:08}.
 The main reason for cloud spreading is thermal diffusion of polaritons towards the edges of
the cloud from the excitation area. Thus, at room temperatures,
 thermal  spreading   is more pronounced than the  effect of the interparticle interactions.

Figure \ref{fig:nmax} shows the maximum density of polaritons, $n_{max}$, reached
at the center of the cloud as a function of the excitation
power $P$. The data are shown for the trapping potential strength  $k=50$ eV/cm$^2$ and temperature $T=300$ K.
It is seen that the maximum density $n_{max}$ nearly linearly grows with the increase of the power $P$.
For the excitation power $P=P_c \approx 600$ mW, the maximum density of polaritons reaches
the critical density $n_c\approx 4.2\times 10^{12}$ m$^{-2}$. At the pumping power $P>P_c$ the central part
of a polariton system at room temperatures should demonstrate superfluid properties previously known for
a low-temperature polariton superfluid in a semiconductor microcavity (see review \cite{Carusotto:13}).
The wings of the spatial distribution always have a density lower than $n_c$ thus,
the superfluid ``lake'' at the cloud center is surrounded by a low-density
gas of polaritons in a normal state.

\begin{figure}[tb]
 \hspace{1 cm } \includegraphics[width=6.5cm]{profile2d.eps}\\
 \includegraphics[width=7cm]{profile.eps}\vspace{-0.5cm}
 \caption{\label{fig:profile} (Color online)
 Spatial distribution of trapped polaritons in a cavity with MoS$_2$ monolayer at $T=300$ K.
 The power of the laser pump and the trapping potential strength are marked in the figure.
 (a) The two-dimensional polariton density profile
 $n(\bm{r})$ averaged over the time domain $t=400-2700$ ps, in which  the system has reached the steady state.
 The color bar shows the polariton density in the units of $10^{12}$ m$^{-2}$.
 (b) A one-dimensional cut of the graph (a) made along the line  $y=0$.
 Points show the results of the simulations,
 vertical bars show the standard deviation about the mean for the density fluctuations.
 In the simulations, the full width at half maximum (FWHM) of the laser spot is 30 $\mu$m; the center of the laser spot
is positioned at the center of the parabolic trap  $\bm{r}=(0,0)$. }
\end{figure}

We emphasize that, after a part of the polariton system transits to a superfluid state,
quantum correlations in the polariton dynamics should be taken into account.
However, if the superfluid density is small compared to the total density of the system,
$n_{max} - n_c < n_c$,
the superfluid effects on the over-all large-scale density distribution in a steady state  can be disregarded
to the first approximation. In particular, the presence of a small superfluid fraction
will not  modify the value of the  maximum density in the cloud and thus,
the criterion of the superfluidity formation will not be changed.
(We note that a few-percent density jump can occur during the superfluid transition, in a similarity with the
normal-to-superfluid transition in He-II \cite{Khalatnikov:65}; in our case however, this jump
is less than or comparable  with a typical value of the thermal  fluctuations of the density $\sim 4$\% estimated from
Fig.\ \ref{fig:profile} and therefore, it can be disregarded.)
We also disregard the dissipative dynamics at a very close vicinity of the superfluid transition
where the superfluid component can appear as a fluctuation mode \cite{Khalatnikov:72}.
Our goal is to pinpoint a domain in the parameter space where the room-temperature superfluidity
of polaritons can be experimentally achieved.

\begin{figure}[h]
\centerline{ \includegraphics[width=8cm]{nmax.eps}}\vspace{-0.5cm}
 \caption{\label{fig:nmax}
Increase in the maximum density of polaritons, $n_{max}$, at the center of the cloud $\bm{r}=(0,0)$
with the rise of the power of the laser pump $P$.
Points show the results of the simulations for $k=50$ eV/cm$^2$, the line segments connecting the points
are shown to guide the eye.
The dashed horizontal line corresponds to the critical density $n_c=4.2 \times 10^{12}$ m$^{-2}$,
at which a gas of interacting polaritons in a cavity with MoS$_2$ undergoes the Kosterlitz-Thouless superfluid transition
at the room temperature.
A polaritons gas is superfluid at the central part of the cloud
above the line ($P>P_c \approx 600$ mW) and  is in the normal state below the line ($P<P_c$).}
\end{figure}

\begin{figure}[tb]
 \centerline{ \includegraphics[width=8cm]{nk.eps}}\vspace{-0.5cm}
 \caption{\label{fig:nk}
Dependence of the maximum density of the polariton cloud on the trapping potential strength $k$.
Points show the results of the simulations,
lines show the theoretical dependence Eq.\ (\ref{eq:mnax}) for the respective parameters.
The  values for the laser pump power are listed in the figure.}
\end{figure}

To better understand the optimal conditions for observations of room-temperature polariton superfluidity,
we studied the dependence of the maximum density  $n_{max}$ on the trapping potential strength.
The results of the numerical simulations are summarized in Fig.\ \ref{fig:nk}.
It is seen that the density  grows with the rise of the trapping potential strength.
From the obtained results it follows that at experimentally realistic pumping power
 $P \sim 600$  mW, the cloud can transit to a superfluid state at $T=300$ K for the trapping strength $k \geq 50$ eV/cm$^2$.
 For a larger power $P = 750$ mW, the polariton density required for the room-temperature superfluidity
 is reached at even  smaller trapping strengths $k\geq 15$ eV/cm$^2$.
 On the other hand, for $P = 300$ mW, the conditions for the room temperature polariton superfluidity
 are not satisfied even for large trapping strengths $k \sim 300$ eV/cm$^2$.

To further investigate the trapped polariton gas dynamics, we also studied
the effect of the spatial offset in the excitation spot position compared to the
center of the parabolic trap.  For this purpose, we numerically studied the polariton spatial distribution
when the center of the laser spot is shifted along the $x$ axis on the distance
$x_0=200$ $\mu$m that is, larger than the excitation spot size of 30 $\mu$m.
The obtained numerical results are shown in Fig.\ \ref{fig:offset}.
Fig.\ \ref{fig:offset}a shows the two-dimensional polariton density distribution for $k=150$ eV/cm$^2$, and
Fig.\ \ref{fig:offset}b shows its cross-section made along the $x$ axis at $y=0$ for $k=150$ eV/cm$^2$
(filled squares) and for $k=50$ eV/cm$^2$ (unfilled triangles).
It is seen that the polaritons diffuse towards the center of the trap $\bm{r}=0$.
It is also seen that the polariton density at the center of the trap
increases with the rise of the trapping potential strength $k$. However, for moderate
strengths $k \leq 150$ eV/cm$^2$, the polariton density at $\bm{r}=0$
remains small compared to the maximum density at the excitation spot area.
Specifically, the density at $\bm{r}=0$ is $\approx 19$\% of the maximum density in the cloud for $k=150$ eV/cm$^2$
and it is only 10\% for $k=50$ eV/cm$^2$.

The geometry with a shifted excitation spot might provide an additional opportunity for the experimental detection of
room-temperature polariton superfluidity. One can expect a larger occupation of the lowest energetic state
near the trap center  $\bm{r}=0$ for a superfluid gas compared to the case of a normal gas (these simulations).
In other words, the superfluid polariton gas would ``easier'' flow towards the bottom of the trap
that, in its turn, should lead to higher densities $n(0)$.
Thus, observation of  ``anomalously'' high polariton population $n(0)$
with the increase of the trapping strength $k$ or the laser power $P$
in experiments with a shifted laser spot can be a simple experimental test of the
normal-to-superfluid transition in the polariton gas.

Finally, we provide an estimate of the  cloud densities that facilitates
the search of  parameters suitable for the polariton superfluidity observation.
As it follows from the above consideration, a key parameter that determines the feasibility of polariton superfluidity
at room temperatures is the maximum polariton density at the center of the cloud $n_{max}$.
To estimate the maximum density $n_{max}$  the following  model can be utilized.
In the case of small trapping potential strengths, the limiting factor for the  cloud size is
a finite polariton life time $\tau_{life}^{(pol)}$.  In this limit, the cloud size can be estimated as
$a=(D \tau_{life}^{(pol)})^{1/2}$. In the opposite case of strong trapping potentials, the
cloud size is limited by thermal fluctuations, $k_B T \sim {1 \over 2} k a^2$, and thus, it
can be estimated as $a = (2 k_B T / k)^{1/2}$. In the intermediate case, the
cloud size can be obtained from the interpolation equation
\begin{equation}
a = {(D\tau_{life}^{(pol)})^{1/2} \over \sqrt{ 1 + \gamma}} \label{eq:a}
\end{equation}
where $\gamma = k D\tau_{life}^{(pol)} / 2 k_B T$ is the dimensionless trapping potential strength.
We note that for the parameters corresponding to an embedded MoS$_2$ monolayer at $T=300$ K (Table \ref{tab:params}), one obtains
 $\gamma = 1$ for $k \approx 440$ eV/cm$^2$.
By approximating the steady-state cloud shape as an axially-symmetric  profile $n(r) = n_{max} /(1+(r/ a)^6)$
and the total number of polaritons in the system as $N \approx P_n \tau_{life}^{(pol)}$, one has the following
expression for the maximum polariton density,
\begin{equation}
 n_{max} = \frac{3 \sqrt{3}}{2 \pi ^2} {P_n \tau_{life}^{(pol)} \over  a^2}
 \approx 0.26 {P_n \tau_{life}^{(pol)} \over  a^2}, \label{eq:mnax}
\end{equation}
in which the parameter $a$ is given in Eq.\ (\ref{eq:a}) and the source $P_n$ is estimated in Eq.\ (\ref{eq:pn}).
The model Eq.\ (\ref{eq:mnax}) is in good agreement with the results of the simulations, as is
seen in Fig.~\ref{fig:nk}.

\begin{figure}[h]
 \hspace{1.5 cm } \includegraphics[width=6cm]{profile2doffset.eps}\\
 \includegraphics[width=7cm]{offset.eps}\vspace{-0.5cm}
 \caption{\label{fig:offset} (Color online)
Polariton diffusion towards the center of the trap  in the case where the excitation
spot is shifted from the trap center. The pumping power is $P=600$ mW.
The offset for the excitation spot center is $x_0=200$ $\mu$m in the positive direction of the
$x$ axis. (a) Two-dimensional polariton density distribution for $k=150$ eV/cm$^2$.
(b) The density distribution calculated along the line $y=0$ for
$k=50$ (unfilled triangles) and 150 eV/cm$^2$ (filled squares). For the reference, the density distribution for
$k=50$ eV/cm$^2$ with zero excitation spot offset, $x_0=0$, is shown by the dashed curve.
Fluctuations of the polariton density ($\sim 4$\%) are not shown.}
\end{figure}

\section{conclusions}
In this paper, we studied the polariton gas dynamics in a parabolic trapping potential by using the quasiclassical
stochastic equation for the polarion wave packets diffusion.
 We showed that in a cavity with an embedded molybdenum disulfide monolayer,
 room-temperature superfluidity in a trapped polariton
cloud can be observed for experimentally
realistic trapping strengths $k \sim 50-300$ eV/cm$^2$ and laser excitation power $P\geq 600$ mW.
The experiments with a shifted laser spot would provide a simple qualitative  test
for superfluidity.
The theoretical estimate  shows that
a trapped polariton gas  is superfluid at the central part of the polariton cloud if
the laser pumping power exceeds  the critical value
\begin{equation}
P_c = {3.8 D n_c \over 1+\gamma} {E_{ph} \over q}.
\end{equation}

In our studies, we are inspired by the recent observation of room temperature polaritons \cite{Liu:15}
as well as by the advances of low-temperature polariton physics \cite{Carusotto:13,Gibbs:11}.
Superfluidity of excitons in transition metal dichalcogenide-based systems
has recently been predicted for temperatures beyond the helium temperature range \cite{Berman:16}.
In addition to the fundamental importance of high-temperature quantum effects,
polariton superfluidity might be a key to  new optoelectronic  technologies including
 those  mentioned in Introduction.
Room-temperature polariton superfluidity would open a new route towards the design of
photonic integrated circuits suitable for the use in optical information processing systems with low energy consumption.
Our hope is that the experimental search for room-temperature quantum  fundamental effects
in transition metal dichalcogenide-based heterostructures would facilitate the progress in this direction.

\section{Funding Information}
National Science Foundation (NSF) (Supplement to 1345219);
Professional Staff Congress -- City University of New York (PSC-CUNY)  (68090-00 46).

The authors thank the Center for Theoretical Physics of the New York City College of
Technology for providing computational resources.


\section{Introduction}\label{sec:intro}

Distributed storage systems form the backbone for modern cloud computing, large--scale data servers, and peer--to--peer systems. The data in these systems is stored in a redundant fashion --- typically via replication (for instance, Hadoop \cite{Hadoop} and Google file systems \cite{GGL03} adopt a triple replication policy) --- to safeguard data against not--so--infrequently occurring disk failures. An alternative approach to storing data on these systems, which highly reduces the redundancy involved in replication, is to use maximum distance separable (MDS) codes such as Reed--Solomon codes. Though MDS codes are the most space--efficient for a targeted worst--case number of simultaneous node failures, they, unlike repetition codes, incur a high repair bandwidth\footnote{\textcolor{black}{A recent work \cite{GW15} revisits this for the case of Reed--Solomon codes, but we do not go into that here.}} when the system undergoes the repair of a single node failure. A new class of erasure codes, called regenerating codes, was recently defined by Dimakis et al. \cite{DGWWR10} over a set of $n$ nodes, which simultaneously optimizes storage efficiency, worst--case resilience and repair bandwidth for single node failures. These codes follow a trade--off curve which is intuitively evidenced by the contrast between repetition codes and MDS codes: the repair bandwidth decreases as the storage redundancy per node increases.

Formally, a file ${\cal F}$ of size $M$, is said to be stored on a DSS consisting of $n$ nodes, each with a storage capacity of $\alpha$, using an $[n,k,d]$-$(\alpha)$ (or, in short, $[n,k,d]$) regenerating code, if it satisfies two properties:
\begin{enumerate}
\item[(a)] {\em data recovery}: the file ${\cal F}$ can be recovered using the contents of any $k$ of the $n$ nodes (this property will also be referred to as the {\em MDS property}); and
\item[(b)] {\em repair property}: the contents of any node can be recovered using the contents of a {\em helper set} of any $d$ other {\em helper} nodes, where each node transmits $\beta$ number of symbols to the replacement node.
\end{enumerate}
An {\em optimal} $[n,k,d]$ regenerating code achieves the optimal value of total repair bandwidth $\gamma = d\beta$ (minimum repair bandwidth) for a given storage capacity $\alpha$ and $M$. This is given implicitly by the following trade--off:
\begin{eqnarray}
M &=& \sum_{i=0}^{k-1} \min\left\{\alpha, (d-i)\beta\right\}.
\end{eqnarray}

Most of the regenerating codes research (e.g.\cite{SRKR12,RSK11,WD09,CDH09,W09,CJMRS13,CHLM11,PDC11,WTB11,WTB12}) is focussed on the extremal points of this trade--off: MBR and MSR codes. {\em Minimum bandwidth regenerating} (MBR) codes achieve the optimal $\alpha$ when the repair bandwidth equals that of a repetition code. This paper concerns {\em minimum storage regenerating} (MSR) codes, often dubbed as {\em optimal bandwidth MDS codes}, because they are optimal regenerating codes that are {\em also} MDS codes\footnote{To be precise, these are vector MDS codes, i.e., MDS codes over $\mathbb{F}_q^{\alpha}$.}. For these codes, $\alpha = M/k$, and the optimal repair bandwidth is given by:
\begin{eqnarray}
\beta &=& \frac{\alpha}{d-k+1}.
\end{eqnarray}

It is easy to see that the total repair bandwidth $d\beta$ is optimized when the number of helper nodes $d = n-1$. However, it is not always practical to connect to {\em all} the remaining nodes to aid the repair of a failed node. We therefore consider the following question:
{\em Are there constructions of $[n,k,d]$ MSR codes, for $d < n-1$?}

\subsection{Previous Work}
This question has not been wholly unanswered. The first
MSR code constructions appeared in \cite{SR10,RSK11}, which roughly correspond to the family of parameters $\{n,k,d\}$ with rate $k/n \le 1/2$. The asymptotic existence of MSR codes for all triples $\{n,k,d\}$ was eventually shown in \cite{CJMRS13} using interference alignment techniques developed for a wireless interference channel; these codes achieve optimality as a regenerating code (as well as approach the MSR point) {\em only} when $\alpha \to \infty$, i.e., $\beta/\alpha \to 1/(d-k+1)$, as $M \to \infty$.

MSR codes, being MDS vector codes, can be expressed as a set of $k$ systematic vectors and $n-k$ parity vectors (the corresponding nodes are referred to as systematic and parity nodes, respectively). For the high--rate ($k/n \ge 1/2$) regime, code constructions
were discovered independently
in \cite{CHLM11,PDC11,TWB11,WTB11} for the specific case of $d = n-1$. Of these, the constructions in \cite{CHLM11,PDC11,TWB11} focus on the relaxation of restricting optimal repair to {\em systematic nodes} in the system; we call the corresponding codes {\em systematic--repair MSR codes}. Practical systems usually store information in a systematic format. Parity nodes may fail, but as in the above works, we do {\em not} require optimal bandwidth repair for such nodes (maybe they are less urgent or critical). Clearly, any node can be repaired by reconstructing the whole file, so this covers the node repairability (even if suboptimally).

\subsection{Contribution \& Outline}
We present the first\footnote{\textcolor{black}{This work was first presented (invited) at the 53rd Annual Allerton
Conference on Communication, Control, and Computing. A simultaneous result was presented at the same venue by Tamo and En Gad \cite{TE15}. Recently and independently, Rawat et al. \cite{RKV16} have constructed MSR codes which optimally repair all nodes. However, the flavor of their construction, which is not systematic in nature, differs from ours.
}} high--rate finite--$\alpha$ constructions for systematic--repair MSR codes for $d < n-1$. We start by describing in Section \ref{sec:preliminaries} the representative code construction that contains the ideas behind those in \cite{CHLM11,PDC11,TWB11}. Leveraging on this, we present our construction in Section \ref{sec:d<n-1}, but restrict to the case when the helper nodes contain the remaining $k-1$ systematic nodes. This restriction is removed in Section \ref{sec:any_helper_nodes}, thus rounding out the code construction. We conclude with some remarks in Section \ref{sec:conclusion}.

\section{Primer: Code Construction for $d = n-1$}\label{sec:preliminaries}

Let $n=k+r$ denote the number of nodes in the distributed storage system, where each node has the capacity to store a vector of size $\alpha$ over $\mathbb{F}_q$.
Throughout this paper, we discuss systematic constructions
and assume that the first $k$ nodes are information nodes and store raw information, while the remaining $r$ nodes correspond to the parities. We use the notation ${\bf x}_i,\, i\in [k],$ for the raw information vectors stored in the systematic nodes. The parity nodes are defined by
\begin{align}\label{eq:parity definition}
{\bf x}_{k+i} = \sum_{j=1}^{k} A_{ij} {\bf x}_{j}, \;\;\;\;\;\;\;\;\; i \in [r],
\end{align}
where $A_{ij}$'s are $\alpha\times\alpha$ encoding matrices. The generator matrix of the code is then given by
\begin{align}\label{eq:generator_matrix}
G = \left[ \begin{array}{ccc}
I 		&		&0		\\
		&\ddots 	&		\\
0		&		&I 		\\
A_{1,1}	&\cdots 	&A_{1,k}	\\
\vdots	 	&\ddots	&\vdots	\\
A_{r,1}	&\cdots 	&A_{r,k}
\end{array} \right].
\end{align}

In this section, we consider MSR codes where $d=n-1$. In other words, when a single node failure occurs, all the remaining nodes aid in its repair.
We also restrict our attention to codes that consider failures only of the systematic nodes, and discuss in this section, a construction that underlies the ideas in \cite{TWB11, PDC11} and \cite{CHLM11}. This construction will inform our generalization for the general parameter triple $\{n,k,d\}$ in Section \ref{sec:d<n-1}.

{\em Remark:} Wang et al. constructed an MSR code for $d = n-1$ in \cite{WTB11} that achieves the optimal repair bandwidth also for parity nodes, albeit at the cost of some other metrics such as the number of symbols read from a node and the complexity of updating parities when systematic data changes. We leave for future the question of whether such a code exists when $d < n-1$.

A commonly adopted strategy in constructing an MSR code is to first guarantee the optimal repair bandwidth property for a single failure (in this case, for a single systematic node failure), and then transform the construction to ensure the MDS property. This is illustrated in Example \ref{ex:423 code} below.

\begin{example}\label{ex:423 code}
Assume $(n,k,d)=(4,2,3)$ and $\alpha=4$. Let the first two nodes ${\bf x}_1$ and ${\bf x}_2$ be the systematic nodes, and let the parity nodes ${\bf x}_3$ and ${\bf x}_4$ be defined as

\[ {\bf x}_3 =\underbrace{ \left( \begin{array}{cccc}
1 	&0	&0	&0	\\
0 	&1	&0	&0	\\
0 	&0	&1	&0	\\
0 	&0	&0	&1
\end{array} \right)}_{I}{\bf x}_1+
\underbrace{\left( \begin{array}{cccc}
1 	&0	&0	&0	\\
0 	&1	&0	&0	\\
0 	&0	&1	&0	\\
0 	&0	&0	&1
\end{array} \right)}_{I}{\bf x}_2,\]
\[{\bf x}_4 = \underbrace{\left( \begin{array}{cccc}
0 	&0	&1	&0	\\
0 	&0	&0	&1	\\
1 	&0	&0	&0	\\
0 	&1	&0	&0
\end{array} \right)}_{P_1}{\bf x}_1+
\underbrace{\left( \begin{array}{cccc}
0 	&1	&0	&0	\\
1 	&0	&0	&0	\\
0 	&0	&0	&1	\\
0 	&0	&1	&0
\end{array} \right)}_{P_2}{\bf x}_2.
 \]

Figure~\ref{fig:423 code}.a depicts the component-wise storage in each node. It can be observed that a single failure in either ${\bf x}_1$ or ${\bf x}_2$ can be reconstructed by downloading $\alpha/2=2$ elements from each of the remaining $d=3$ nodes. However, the data is not recoverable if both ${\bf x}_1$ and ${\bf x}_2$ fail and hence, the code is not MDS. To overcome this problem, we associate a coefficient $\lambda$ with $P_2$ such that $\small{\begin{pmatrix} I & I\\  P_1 & \lambda P_2  \end{pmatrix}}$ is non-singular. Note that,
\begin{align*}
\begin{vmatrix} I &\hspace{-1.5ex} I\\
  P_1 & \hspace{-1.5ex}\lambda P_2  \end{vmatrix} \hspace{-.5ex}=\hspace{-.3ex} \det(\lambda P_2 - P_1)\hspace{-.5ex}=\hspace{-.5ex}
\begin{vmatrix}
0 		&\hspace{-.7ex}\lambda	&\hspace{-.7ex}-1		&\hspace{-.7ex}0	\\
\lambda 	&\hspace{-.7ex}0		&\hspace{-.7ex}0		&\hspace{-.7ex}-1	\\
-1 		&\hspace{-.7ex}0		&\hspace{-.7ex}0		&\hspace{-.7ex}\lambda	\\
0 		&\hspace{-.7ex}-1		&\hspace{-.7ex}\lambda	&\hspace{-.7ex}0
\end{vmatrix}\hspace{-.5ex}=\hspace{-.5ex} (\lambda^2 -1)^2,
\end{align*}
which is non-zero\footnote{In general, if $A,B,C,$ and $D$ are nonsingular $\alpha\times\alpha$ matrices, then $\det\big({\begin{bmatrix}A &B\\ C &\lambda D\end{bmatrix}}\big)$ is given by $\det(A)\det(\lambda D - CA^{-1}B)$, which is a polynomial of degree at most $\alpha$ in $\lambda$. If the field size is large enough, \emph{i.e.} $q>\alpha$, one can always find a value for $\lambda$ so that the $2\times 2$ block matrix becomes non-singular as well. The same approach can be used to prove Lemma~\ref{lemma:MDS property2}. } if $q=5$ and $\lambda = 2$. Figure~\ref{fig:423 code}.b shows the component-wise storage for the resulting MSR code.

\begin{figure}[t]
	\centering
	\includegraphics[scale=.5]{Example1.pdf}
	\caption{{(a) Component wise storage in a $(4,2,3)$ binary array code with optimal repair bandwidth for a single systematic node failure, described by $({\bf x}_1, {\bf x}_2, {\bf x}_1 + {\bf x}_2, P_1{\bf x}_1 + P_2{\bf x}_2)$; (b) A $(4,2,3)$ MSR code in $\mathbb{F}_5$ described by $({\bf x}_1, {\bf x}_2, {\bf x}_1 + {\bf x}_2, P_1{\bf x}_1 +2 P_2{\bf x}_2)$. In both cases, gray cells are accessed to rebuild $C_1$.
	\label{fig:423 code}}}
\end{figure}
\end{example}

Construction \ref{cons:d_n-1_systematic} generalizes the construction given in Example \ref{ex:423 code} for an $(n,k,n-1)$ MSR code. Note that any MSR code construction must specify both the generator matrix of the code as well as the optimal bandwidth repair strategy that is implemented on the code.

\begin{cons}\label{cons:d_n-1_systematic} Let $\alpha = r^k$ and label the $\alpha$ elements $[0:r^k-1]$ by $r$-ary vectors in $\mathbb{Z}_r^k$.  Define permutation $f_j^{\ell}$ on $[0:r^k-1]$ as follows:
\begin{eqnarray*}
\begin{array}{cccl}
f_j^{\ell}:&\mathbb{Z}_r^k &\to& \mathbb{Z}_r^k\\
&v &\mapsto& v + \ell e_j,
\end{array}
\end{eqnarray*}
for $j\in [k]$ and $\ell \in [0:r-1] := \{0,1,\ldots,r-1\}$, where $\{e_1,e_2,\ldots,e_k\}$ is the standard vector basis for $\mathbb{Z}_r^k$.
The mapping $f_j^{\ell}$ is bijective, and therefore, corresponds to a permutation on $[0:r^k-1]$.
Let $P_{\ell,j}$ be the $\alpha\times\alpha$ matrix corresponding to the permutation $f_j^{\ell}$, that is, $P_{\ell,j}\,{\bf x} = {\bf y}$, where ${\bf x}$, ${\bf y} \in \mathbb{F}_q^{\alpha}$, and ${\bf x}(v) = {\bf y}(f_j^{\ell}(v))$. In other words, $P_{\ell,j}$ scrambles the elements of a vector according to the permutation $f_j^{\ell}$. (Notice that $P_{0,j} = I_{\alpha}$.)
\begin{enumerate}
\item {\em MSR Code}: The generator matrix of the code is given by (\ref{eq:generator_matrix}), where $A_{i,j} = \lambda_{i,j}P_{i-1,j}$, $i \in [r]$ and $j \in [k]$.
The non-zero coefficients $\lambda_{i,j}\in \mathbb{F}_q$ will be defined in Section \ref{sec:MDS} to ensure the MDS property.

\item {\em Repair Strategy}: Let $Y_j=\{v\in[0,r^k-1]:v \cdot e_j=0 \}$ denote a subset of $[0:r^k-1]$. $Y_j$ can be interpreted as those elements in $[0:\alpha-1]$ whose label representation in $\mathbb{Z}_r^k$ have a $0$ in their $j^{\text{th}}$ coordinate. If systematic node $j$ fails, it is repaired by accessing the elements corresponding to $Y_j$ from each of the remaining nodes, i.e., by accessing ${\bf x}_i(v)$, where $v \in Y_j$ and $j \neq i \in [n]$.
\end{enumerate}
\end{cons}

Construction \ref{cons:d_n-1_systematic} is obtained by first constructing an $[n,k]$ array code\footnote{By an $[n,k]$ array code, we mean a set of $k$ systematic vectors, and $n-k$ parity vectors defined according to (\ref{eq:parity definition}), which may or may not satisfy any properties.} (Section \ref{sec:IA}) which guarantees the optimal bandwidth repair for a single systematic node failure. The array code is then transformed (Section \ref{sec:MDS}) to an MDS array code (and thereby, a systematic--repair MSR code) by transforming the encoding matrices of the parity nodes, while retaining the repair property.

\subsection{Repair Property: Interference Alignment}\label{sec:IA}
The optimal repair bandwidth property of an $[n,k,n-1]$ MSR code can be viewed as a {\em signal interference} problem: the objective is to retrieve the desired signal --- the contents of the failed systematic node, say, ${\bf x}_i$ --- which, in the repair data downloaded from the remaining nodes, is interfered by partial contents of the remaining systematic nodes, ${\bf x}_j$, where $i \neq j \in [n]$. The solution, turns out to be an {\em interference alignment} strategy, where the repair data associated with the interfering systematic data is aligned, so as to minimize the interference. This is crystallized in the following lemma\footnote{This result is known and has been used in several papers on MSR codes, but we state and prove it for completeness.}.

\begin{lemma}\label{lemma:interference alignment_n-1}
Let ${\bf x}_i$, $i \in [k]$, be the failed systematic node. For an $[n,k,n-1]$ MSR code, the set of $d = n-1$ helper nodes is given by ${\cal D} = \{{\bf x}_j \,|\, j \in [n]\backslash \{i\}\}$. To recover the contents of the failed systematic node with the optimal repair bandwidth, it is necessary and sufficient to find $n-1$ (repair) matrices denoted by
 $\{S_j^i\in \mathbb{F}_q^{\alpha/r \times \alpha}\,|\,j \in [n]\backslash \{i\}\}$, where $r = n-k$, such that, for $j \in [k], j \neq i$, the following two conditions are satisfied:

\noindent (a) signal recovery:
\begin{eqnarray}\label{eq:signal-recovery}
\mathsf{rank}\left(\left(\begin{array}{c}
S_{k+1}^iA_{1,i}	\\
S_{k+2}^iA_{2,i}	\\
\vdots 	\\
S_{k+r}^iA_{r,i}
\end{array}\right)\right) &=& \alpha,
\end{eqnarray}
(b) interference alignment:
\begin{eqnarray}\label{eq:IAcondition}
\mathsf{rank}\left(\left(\begin{array}{c}
S_j^i	\\
S_{k+1}^iA_{1,j}	\\
\vdots 	\\
S_{k+r}^iA_{r,j}
\end{array}\right)\right) &=& \frac{\alpha}{r}.
\end{eqnarray}

Stated otherwise, to optimally repair ${\bf x}_i$, it is necessary and sufficient to find $n-1$ (repair) subspaces of dimension $\alpha/r$, denoted\footnote{Whenever this lemma is referenced, we use the subspace and matrix notation interchangeably as some proofs or expressions are clearer in one of the formats. We accordingly overload the notation $S_j^i$ to refer to both the matrix and the subspace spanned by the row vectors of the matrix.} by $\{S_j^i\,|\,j \in [n]\backslash \{i\}\}$, where $r = n-k$, such that, for $j \in [k], j \neq i$, the following two conditions are satisfied:

\noindent (a) signal recovery:
\begin{eqnarray}\label{eq:signal-recovery2}
S_{k+1}^iA_{1,i} \oplus \cdots \oplus S_{k+r}^iA_{r,i} &\backsimeq& \mathbb{F}_q^{\alpha},
\end{eqnarray}
(b) interference alignment:
\begin{eqnarray}\label{eq:IAcondition2}
S_j^i &\backsimeq& S_{k+s}^iA_{s,j},\,\,\forall \,s \in [r],
\end{eqnarray}
where $\backsimeq$ denotes equality of subspaces, $SA$ is the subspace obtained by operating the subspace $S$ by the matrix $A$, and $\oplus$ denotes the subspace sum.
\end{lemma}

For completeness, we provide a proof for Lemma \ref{lemma:interference alignment_n-1} in Appendix \ref{app:lemmaIA}. Lemma \ref{lemma:interference alignment_general} generalizes Lemma \ref{lemma:interference alignment_n-1} when the number of helper nodes $d < n-1$. This will be used later in Section \ref{sec:d<n-1}.

\begin{lemma}\label{lemma:interference alignment_general}
{\em (Corollary of Lemma \ref{lemma:interference alignment_n-1}.)}
In general, for an $[n,k,d]$ MSR code, if the set of $d=(k-1)+t<n-1$ helper nodes is given by ${\cal D} = \{{\bf x}_j \,|\, j \in {\cal J} = [k]\backslash \{i\}\cup\{b_1,\cdots,b_t\}\}$ (where $b_i \in \{k+1,\ldots,n\}$ denote the $t$ parity nodes in the helper set), it is necessary and sufficient to find $d$ (repair) subspaces of dimension $\alpha/t$ denoted by $\{S_j^i\,|\,j \in {\cal J}\}$,
such that, for $j \in [k], j \neq i$, the following two conditions are satisfied:
\begin{eqnarray}
S_{b_1}^iA_{b_1-k,i} \oplus \cdots \oplus S_{b_t}^iA_{b_t-k,i} &\backsimeq& \mathbb{F}_q^{\alpha},
\end{eqnarray}
\begin{eqnarray}
S_j^i &\backsimeq& S_{b_s}^iA_{b_s-k,j},\,\,\forall \,s \in [t].
\end{eqnarray}
\end{lemma}

The optimal repair property of Construction \ref{cons:d_n-1_systematic} can now be justified.

\begin{figure}[t]
	\centering
	\includegraphics[scale=.33]{Inteference_Alignment.pdf}
	\caption{Visualization of Lemma \ref{lemma:interference alignment_n-1}(a)[left], and Lemma \ref{lemma:interference alignment_n-1}(b)[right] to justify repair optimality in Construction~\ref{cons:d_n-1_systematic}.}
\end{figure}

\begin{lemma}\label{lemma:optimal-repair-d_n-1}
The repair strategy in Construction \ref{cons:d_n-1_systematic} is optimal with respect to repair bandwidth.
\end{lemma}
\begin{IEEEproof}
Define $S_j^i\triangleq S^i\triangleq Y_i, j \neq i$. Notice that the rank of subspace $S^i$ is $r^{k-1}=\alpha/r$. Per definition, the permutation $P_{\ell,i}$ maps $Y_i$ to $Y_i + \ell e_i = \{v\in[0,r^k-1]:v \cdot e_i=\ell \}$. This implies that for any distinct $\ell, \ell' \in [0:r-1]$, the intersection $S^iP_{\ell,i} \cap S^iP_{\ell',i}$ contains only the all-zero vector. Thus the subspaces: $S^i, S^iP_{1,i}, \ldots, S^iP_{r-1,i}$, span the space $\mathbb{F}_q^{\alpha}$ ($\alpha = r^k$) and the signal recovery condition(s) in Lemma \ref{lemma:interference alignment_n-1} are satisfied. Furthermore, applying a permutation $P_{\ell,j}$ corresponding to a different coordinate $j \neq i$ maps $Y_i$ to itself.
This validates the interference alignment condition(s) in Lemma \ref{lemma:interference alignment_n-1}. Finally, note that the two conditions continue to be satisfied when replacing the permutations $P_{i-1,j}$ with any scaled versions $A_{i,j} = \lambda_{i,j}P_{i-1,j}$, because the scaling of the basis vectors does not change the relevant subspaces and thereby does not affect the conditions in Lemma \ref{lemma:interference alignment_n-1}.
\end{IEEEproof}

\subsection{MDS Property}\label{sec:MDS}
This second step relies on the following two lemmas, the proofs of which are left to the reader.
\begin{lemma}\label{lemma:MDS property1}
Let $B$ denote the parity part of the generator matrix for an $[n,k]$ array code denoted by $\mathbb{C}$, where
\[ B = \left[ \begin{array}{ccc}
B_{1,1}	&\cdots	&B_{1,k}	\\
\vdots		&\ddots 	&\vdots	\\
B_{r,1}	&\cdots	&B_{r,k}
\end{array} \right].\]
Given that $B_{i,j}$ is non-singular for all $i,j$, then $\mathbb{C}$ is an MDS array code if and only if any square sub-block-matrix $B'$ of $B$ is also non-singular, where
\[ B' = \left[ \begin{array}{ccc}
B_{i_1,j_1}	&\cdots	&B_{i_1,j_t}	\\
\vdots		&\ddots 	&\vdots	\\
B_{i_t,j_1}	&\cdots	&B_{i_t,j_t}
\end{array} \right],\;\;\]
for some $\{i_1,\cdots,i_t\}\subset [r], \{j_1,\cdots,j_t\}\subset [k]$.
\end{lemma}

\begin{lemma}\label{lemma:MDS property2}
Let $B$ denote the $r\alpha\times k\alpha$ matrix associated with the parity part of the generator matrix for an $[n,k]$ array code, as defined in Lemma~\ref{lemma:MDS property1}.
Given that $B_{i,j}$ is non-singular for all $i \in [r]$, $j \in [k]$, and the field size $q$ is large enough, there exist coefficients $\lambda_{i,j}\in \mathbb{F}_q$,
such that all square sub-block-matrices of $A$  are non-singular, where
\[ A = \left[ \begin{array}{ccc}
\lambda_{1,1}B_{1,1}	&\cdots	&\lambda_{1,k}B_{1,k}	\\
\vdots		&\ddots 	&\vdots	\\
\lambda_{r,1}B_{r,1}	&\cdots	&\lambda_{r,k}B_{r,k}
\end{array} \right].\]
In other words, any parity generator matrix $B$ for an $[n,k]$ array code with non-singular encoding matrices can be transformed into a parity generator matrix $A$ for an $[n,k]$ MDS array code by multiplying the encoding matrices with appropriate scalar coefficients.
\end{lemma}

\begin{IEEEproof}[\textcolor{black}{Proof Sketch}] To obtain a valid set of $\lambda_{i,j}$'s, one may first sort the pairs $(i,j)$ with respect to $i+j$ increasingly, and then recursively choose a value for each $\lambda_{i,j}$ such that all sub-block-matrices with $\lambda_{i,j}A_{i,j}$ on their bottom right corner become non-singular. It suffices to have the field size $q$ greater than the number of such sub-block-matrices at any step multiplied by $\alpha$;
\begin{eqnarray*}
|\mathbb{F}|>q_{\text{\tiny{MDS}}}=\alpha\max_{t}\left\{\displaystyle\hspace{-.5ex}\binom{\hspace{-.2ex}n-k-1\hspace{-.2ex}}{t}\hspace{-.5ex}\times\hspace{-.5ex}\binom{k-1}{t}\bigg|t\in[k]\right\}\hspace{-.6ex}. \hspace{0.55cm
\end{eqnarray*}
\end{IEEEproof}

\section{Code Construction for Restricted Helper Set}\label{sec:d<n-1}
We now move to the construction of $[n,k,d]$ systematic--repair MSR codes for any $n$, $k$, and $d$, where $k+1\leq d \leq n-1$.
In this section, we start with the restricted case when the helper set $\mathcal{D}$ includes all remaining $k-1$ systematic nodes. Let us begin with an example.
\begin{example}\label{example:k+3_k_k+1}
Let us look at the case when $[n,k,d]= [k+3,k,k+1]$ for $k\in\mathbb{N}$. Given a failure at the systematic node $i$, we are interested in repairing it by downloading $\frac{\alpha}{d-k+1} = \frac{\alpha}{2}$ symbols from each node in the helper set $\mathcal{D}_i$. Let us assume that $\mathcal{D}_i$ includes all of the remaining $k-1$ systematic nodes. Hence, there are $\binom{3}{2}=3$ different ways to choose ${\cal D}_i$ depending on which two parity nodes are included in it. Let us use an indicator $a\in[3]$ to differentiate between these scenarios, and denote the helper set for each scenario by $\mathcal{D}_{i,a}$.

\emph{Construction.} Let $\alpha=2^{3k}$ and label the $\alpha$ elements $[0:2^{3k}-1]$ by binary vectors in $\mathbb{Z}_2^{3k}$. Define permutation $f_{j}^{\ell}$ on $[0:2^{3k}-1]$ as follows:
\begin{align*}
&f_{j}^\ell : 	&\hspace{-7ex}\mathbb{Z}_2^{3k} 	&\hspace{2ex}\to 		&\hspace{-9ex}\mathbb{Z}_2^{3k}\hspace{3.3ex}\\
&		&\hspace{-7ex}v 				&\hspace{2ex}\mapsto	&\hspace{-9ex}v+\ell e_i,
\end{align*}
for $j\in[3k]$ and $\ell\in\{0,1\}$, where $\{e_1,e_2,\cdots,e_{3k}\}$ is the standard vector basis for $\mathbb{Z}_2^{3k}$. The mapping $f_j^{\ell}$ is again bijective and therefore corresponds to a permutation on $[0:2^{3k}-1]$. As before, let $P_{\ell,j}$ be the $\alpha\times\alpha$ matrix corresponding to the permutation $f_j^{\ell}$, that is, $P_{\ell,j}{\bf x} = {\bf y}$, where ${\bf x}, {\bf y} \in \smash{\mathbb{F}_{\!q}}^{\alpha}$, and ${\bf x}(v) = {\bf y}(f_j^{\ell}(v))$. (Notice again that $P_{0,j} = I_{\alpha}$.)

\begin{enumerate}
\item {\em MSR Code}: The generator matrix of the code is given by
\[ G = \left[ \begin{array}{ccc}
I 		&		&0		\\
		&\ddots 	&		\\
0		&		&I 		\\
A_{1,1}	&\cdots 	&A_{1,k}	\\
A_{2,1}	&\cdots 	&A_{2,k}	\\
A_{3,1}	&\cdots 	&A_{3,k}
\end{array} \right],\]
where
\begin{align}\label{eq:example_2_construction}
&A_{1,j}=\lambda_{1,j}\hspace{1ex}\times\hspace{1ex} P_{0,3j-2} \hspace{1ex}\times \hspace{1ex} P_{0,3j-1} \hspace{1ex}\times I_{\alpha},\nonumber\\\nonumber
&A_{2,j}=\lambda_{2,j}\hspace{1ex}\times \hspace{1ex}P_{1,3j-2} \hspace{1ex}\times  \hspace{1ex}I_{\alpha}   \hspace{5.15ex}\times P_{0,3j},\\
&A_{3,j}=\lambda_{3,j}\hspace{1ex}\times \hspace{1ex} I_{\alpha} \hspace{5.15ex}\times \hspace{1ex}P_{1,3j-1} \hspace{1ex}\times P_{1,3j},
\end{align}
for $j\in [k]$. The non-zero coefficients $\lambda_{i,j}\in \smash{\mathbb{F}_{\!q}}$ are again selected according to the discussion in Section \ref{sec:MDS} to establish the MDS property.
\item {\em Repair Strategy via  $\mathcal{D}_{i,1} = \{{\bf x}_j | j\in [k+3], j\neq i, k+3 \}$}:

Let $Y_{i,1}=\{v\in[0,2^{3k}-1]:v \cdot e_{3i-2}=0 \}$ denote a subset of $[0:2^{3k}-1]$. $Y_{i,1}$ can be interpreted as those elements in $[0:2^{3k}-1]$ whose label representation in $\mathbb{Z}_2^{3k}$ have a $0$ in their $(3i-2)^{\text{th}}$ coordinate. If systematic node $i$ fails, it can be repaired by accessing the elements corresponding to $Y_{i,1}$ from each of the helper nodes, i.e., by accessing ${\bf x}_j(v)$, where $v \in Y_{i,1}$ and $j \in \{1,2,\cdots,i-1,i+1,\cdots,k,k+1,k+2\}$.
\item {\em Repair Strategy via $\mathcal{D}_{i,2} = \{{\bf x}_j | j\in [k+3], j\neq i, k+2 \}$}:

Similarly, let $Y_{i,2}=\{v\in[0,2^{3k}-1]:v \cdot e_{3i-1}=0 \}$. If systematic node $i$ fails, it can be repaired by accessing ${\bf x}_j(v)$, where $v \in Y_{i,2}$ and $j \in \{1,2,\cdots,i-1,i+1,\cdots,k,k+1,k+3\}$.
\item {\em Repair Strategy via $\mathcal{D}_{i,3} = \{{\bf x}_j | j\in [k+3], j\neq i, k+1 \}$}:

Finally, let $Y_{i,3}=\{v\in[0,2^{3k}-1]:v \cdot e_{3i}=0 \}$ denote the location of the elements that have to get accessed if the systematic node $i$ fails, i.e., node $i$ can be repaired by accessing ${\bf x}_j(v)$, where $v \in Y_{i,3}$ and $j \in \{1,2,\cdots,i-1,i+1,\cdots,k,k+2,k+3\}$. \\
\end{enumerate}
\begin{IEEEproof}[Justification of the repair strategy]
Let ${\bf x}_i$, $i\in[k]$, be the failed systematic node. Define $Q_{u,v}=A_{u,v}\lambda_{u,v}^{-1}$, $u\in[3], v\in[k]$, which is a product of multiple permutation matrices, and hence can be viewed as a permutation matrix itself. In order to justify the repair strategy, it suffices to define the proper subspaces $S_j^i$ that fulfill the two interference alignment conditions in Lemma \ref{lemma:interference alignment_general}. Let $U_{i,a}$ be the complimentary subset of $Y_{i,a}$ in $\mathbb{Z}_2^{3k}$, i.e.,
\begin{align*}
&U_{i,1}=\{v\in[0,2^{3k}-1]:v \cdot e_{3i-2}=1 \},\\
&U_{i,2}=\{v\in[0,2^{3k}-1]:v \cdot e_{3i-1}=1 \},\\
&U_{i,3}=\{v\in[0,2^{3k}-1]:v \cdot e_{3i}=1 \}.
\end{align*}
Given the code construction in (\ref{eq:example_2_construction}), we can verify that
\begin{align}\label{eq:example2_justification}
Y_{i,1}Q_{1,i} &= Y_{i,1},&	Y_{i,2}Q_{1,i} &= Y_{i,2},	&Y_{i,3}Q_{1,i} &= Y_{i,3},\nonumber\\
U_{i,1}Q_{1,i} &= U_{i,1},&	U_{i,2}Q_{1,i} &= U_{i,2},	&U_{i,3}Q_{1,i} &= U_{i,3},\nonumber\\
Y_{i,1}Q_{2,i} &= U_{i,1},&	Y_{i,2}Q_{2,i} &= Y_{i,2},	&Y_{i,3}Q_{2,i} &= Y_{i,3},\nonumber\\
U_{i,1}Q_{2,i} &= Y_{i,1},&	U_{i,2}Q_{2,i} &= U_{i,2},	&U_{i,3}Q_{2,i} &= U_{i,3},\nonumber\\
Y_{i,1}Q_{3,i} &= Y_{i,1},&	Y_{i,2}Q_{3,i} &= U_{i,2},	&Y_{i,3}Q_{3,i} &= U_{i,3},\nonumber\\
U_{i,1}Q_{3,i} &= U_{i,1},&	U_{i,2}Q_{3,i} &= Y_{i,2},	&U_{i,3}Q_{3,i} &= Y_{i,3},\nonumber\\
 \hspace{-2ex}\text{and,}  & & & & \nonumber\\
Y_{i',u}Q_{i,j} &= Y_{i',u}\hspace{.5ex},\hspace{-.5ex}&		U_{i',u}Q_{i,j} &= U_{i',u}\hspace{.5ex},\hspace{-.5ex}		&\text{for } &i'\neq i.
\end{align}
Now we define subspaces $S_{j,a}^i \triangleq S_{a}^i \triangleq Y_{i,a}, j\neq i, a\in [3]$. Let us for simplicity assume $a=1$. The other scenarios follow the proof similarly. Based on (\ref{eq:example2_justification}), we observe that the permutation $Q_{2,i}$ maps the basis $Y_{i,1}$ to its complementary subset $U_{i,1}$ and vice versa, while $Q_{1,i}$ preserves both of them. Hence,
\begin{align*}
\mathsf{rank}\left( \left( \begin{array}{c}
S_1^iQ_{1,i}	\\
S_1^iQ_{2,i}
\end{array} \right) \right) =
\mathsf{rank}\left( \left( \begin{array}{c}
Y_{i,1}	\\
U_{i,1}
\end{array} \right) \right) =\alpha.
\end{align*}
Furthermore, $Y_{i,1}$ remains unchanged under any other permutation $Q_{t,j}, j\neq i$, and hence
\begin{align*}
\mathsf{rank}\left( \left( \begin{array}{c}
S_1^i	\\
S_1^iQ_{1,i'}	\\
S_1^iQ_{2,i'}
\end{array} \right) \right) =
\mathsf{rank}\left( \left( \begin{array}{c}
Y_{i,1}	\\
Y_{i,1}	\\
Y_{i,1}
\end{array} \right) \right) =\frac{\alpha}{2}.
\end{align*}
\end{IEEEproof}
\end{example}
The key element in the construction is to satisfy the two requirements in Lemma \ref{lemma:interference alignment_general} for any systematic failure and any such helper set $\mathcal{D}$. Let $\rho=d-k+1$ denote the number of parity nodes in the helper set of size $d$. There are $\binom{r}{\rho}$ different ways to choose $\rho$ parity nodes during the repair. Let us label these cases with numbers $a\in [\binom{r}{\rho}]$, and set $\mathcal{R}_a$ to be the subset of parity nodes corresponding to case $a$.

Assume that $\mathcal{R}_a=\{{\bf x}_{k+d_1^{(a)}},{\bf x}_{k+d_2^{(a)}},\cdots,{\bf x}_{k+d_{\rho}^{(a)}}\}$ is the ordered representations, where $\{d_1^{(a)},\cdots,d_{\rho}^{(a)}\} \subset [r]$. Finally, define $r$-ary vectors ${\omega}_{a}$ for $a\in[\binom{r}{\rho}]$ as
$$\omega_a(i) = \left\{
	\begin{array}{lll}
		t-1  && \mbox{if } \exists t:  i=d_t^{(a)},\\
		0 && \mbox{otherwise.}
	\end{array}
\right. $$

\begin{cons}\label{cons:fixed_d}
Let $\alpha =\rho^{k\binom{r}{\rho}}$ and label the $\alpha$ elements $[0:\alpha-1]$ by $\rho$-ary vectors in $\mathbb{Z}_{\rho}^{k\binom{r}{\rho}}$.  Define permutation $f_j^{\ell}$ on $[0:\alpha-1]$ as follows:
\begin{eqnarray*}
\begin{array}{cccl}
f_j^{\ell}:&\mathbb{Z}_{\rho}^{k\binom{r}{\rho}} &\to& \mathbb{Z}_{\rho}^{k\binom{r}{\rho}}\\
&v &\mapsto& v + \ell e_j,
\end{array}
\end{eqnarray*}
for $j \in [k\binom{r}{\rho}]$ and $\ell \in [0:\rho-1]$, where $\{e_1,\cdots,e_{k\binom{r}{\rho}}\}$ is the standard vector basis of $\mathbb{Z}_{\rho}^{k\binom{r}{\rho}}$. Let $P_{\ell,j}$ be the $\alpha\times\alpha$ matrix corresponding to the permutation $f_j^{\ell}$.
\begin{enumerate}
\item {\em MSR Code}: The generator matrix of the $[n,k,d]$ code is given by (\ref{eq:generator_matrix}), where
\begin{align*}
A_{i,j}=\lambda_{i,j}\prod_{a\in [\binom{r}{\rho}]}
P_{w_{a}(i),\hspace{.2ex}a+(j-1)\binom{r}{\rho}}, &\hspace{2ex}\text{for } j\in[k],i\in[r].
\end{align*}
The non-zero coefficients $\lambda_{i,j}\in \mathbb{F}_q$ are defined according to Section \ref{sec:MDS} to ensure the MDS property; and later will be modified again in Section \ref{sec:any_helper_nodes}.
\item {\em Repair Strategy}: Let ${\cal R}_a$ correspond to the parity subset of the helper set $\mathcal{D}$. Define $Y_{j,a}\subset[0:\alpha-1]$  as $\{x\in[0,\alpha-1]:x \cdot e_{a+(j-1)\binom{r}{\rho}}=0 \}$. If systematic node $j$ fails, it is repaired by accessing the elements corresponding to $Y_{j,a}$ from helper nodes, i.e., by accessing ${\bf x}_i(v)$, where $i\in\mathcal{D}$, and $v \in Y_{j,a}$.
\end{enumerate}
\end{cons}

\begin{lemma}\label{thm:fixed_d}
The repair strategy in {\em Construction \ref{cons:fixed_d}} is optimal with respect to repair bandwidth.
\end{lemma}

\begin{figure}[t]
\centering
\includegraphics[scale=0.5]{Example2.pdf}
\caption{Relation between $\omega_a$ and (\ref{eq:example_2_construction}).  }
\label{fig:example_2}
\end{figure}

\begin{IEEEproof}
Let us first explain the role of $\omega_a$ by revisiting Example \ref{example:k+3_k_k+1} via Figure \ref{fig:example_2}. Here we assumed that
\begin{align*}
a=1 \rightarrow \hspace{1ex }\mathcal{R}_1=\{{\bf x}_{k+1},{\bf x}_{k+2}\}\hspace{1ex}\rightarrow \hspace{1ex}\omega_1 = (0,1,0)^t,\\
a=2 \rightarrow \hspace{1ex }\mathcal{R}_2=\{{\bf x}_{k+1},{\bf x}_{k+3}\}\hspace{1ex}\rightarrow \hspace{1ex}\omega_2 = (0,0,1)^t,\\
a=3 \rightarrow \hspace{1ex }\mathcal{R}_3=\{{\bf x}_{k+2},{\bf x}_{k+3}\}\hspace{1ex}\rightarrow \hspace{1ex}\omega_3 = (0,0,1)^t.
\end{align*}
In general, the matrix $\Omega=\{\omega_1 | \omega_2 | \cdots | \omega_{\binom{r}{\rho}}\}$ is designed in a way that for any choice of $a\in [\binom{r}{\rho}]$ we can always find a column in $\Omega$, denoted by $\omega_a$, such that its intersection with $r'$ rows associated with scenario $a$, forms $\{0,1,\cdots, \rho-1\}$.

Now assume that node $i$ is failed and we are to perform an optimal systematic repair given parity repairs in $\mathcal{R}_a=\{{\bf x}_{k+d_1^{(a)}},{\bf x}_{k+d_2^{(a)}},\cdots,{\bf x}_{k+d_{\rho}^{(a)}}\}$. It is now clear that if we selected our subspaces as $S_{j,a}^i \triangleq S_{a}^i \triangleq Y_{i,a}=\{x|x \cdot e_{a+(i-1)\binom{r}{\rho}}=0 \}$, then
\begin{align*}
&Y_{i,a}Q_{d_1^{(a)},i} = &\{x|x \cdot e_{a+(i-1)\binom{r}{\rho}}=0 \},\hspace{4ex}\\
&Y_{i,a}Q_{d_2^{(a)},i} = &\{x|x \cdot e_{a+(i-1)\binom{r}{\rho}}=1 \},\hspace{4ex}\\
&\hspace{4ex}\vdots\\
&Y_{i,a}Q_{d_{\rho}^{(a)},i} = &\{x|x \cdot e_{a+(i-1)\binom{r}{\rho}}=\rho-1 \},
\end{align*}
and hence,
\begin{align*}
\mathsf{rank}\left( \left( \begin{array}{c}
S_a^iQ_{d_1^{(a)},i}	\\
S_a^iQ_{d_2^{(a)},i}	\\
\vdots\\
S_a^iQ_{d_{\rho}^{(a)},i}
\end{array} \right) \right) = \rho\times
\mathsf{rank}\left( Y_{i,a}\right) =\rho \frac{\alpha}{\rho} = \alpha.
\end{align*}
The second condition in Lemma \ref{lemma:interference alignment_general} is also automatically satisfied since
\begin{align*}
Y_{i,a} \simeq Y_{i,a}Q_{1,i'} \simeq Y_{i,a}Q_{2,i'} \simeq \cdots \simeq Y_{i,a}Q_{\binom{r}{\rho},i'} \hspace{2ex} \text{for } i'\neq i.
\end{align*}
\end{IEEEproof}

Lastly, we note that optimizing the sub-packetization parameter, $\alpha$, is not the main concern. Although Construction ~\ref{cons:fixed_d} suggests a fairly large value, i.e. $\alpha=\rho^{k\binom{r}{\rho}}$, but it is clear that we do not need $\binom{r}{\rho}$ many columns in $\Omega$ to cover all the $\binom{r}{\rho}$ helper set selection scenarios. Indeed, $\alpha$ in Example \ref{example:k+3_k_k+1} can be reduced to $2^{2k}$, where $\Omega_{\text{new}} = \{\omega_1|\omega_2\}$. We leave the optimizations of this kind to future work.

\section{Code Construction for any Helper Set}\label{sec:any_helper_nodes}

In this section, we show that Construction \ref{cons:fixed_d} in fact holds, even when an arbitrary set of $d$ helper nodes is allowed to be chosen from the $(n-1)$ surviving nodes. This generality merely imposes some additional constraints on the selection of the scaling coefficients $\lambda_{i,j}$ of the encoding matrices $A_{i,j} = \lambda_{i,j}Q_{i,j}$, where $Q_{i,j}$ is the (product) permutation matrix corresponding to $A_{i,j}$, as defined in Construction \ref{cons:fixed_d}. We now arrive at the main theorem.

\begin{theorem}{\label{thm:any_d}} {\em Construction \ref{cons:fixed_d}} gives an $[n,k,d]$ systematic--repair MSR code for any set of $d$ helper nodes, for a large enough field size for the scaling coefficients $\lambda_{i,j}$ for the encoding matrices $A_{i,j}$.
\end{theorem}

\begin{IEEEproof}\emph{ Part 1:} First, we illustrate the proof by fixing $d = k+1$, and taking an example set of helper nodes for an example failure of node ${\bf x}_1$ (or node $1$). Let us denote the (indices of the) helper set by ${\cal D}$, and let ${\cal D} = \{h,h+1,\ldots,k,k+1,\ldots,k+h\}$,
that is, there are $h$ parity nodes and $d-h = k+1-h$ systematic nodes in the helper set. Let $S_i^j({\cal D}){\bf x}_i$ denote the repair information
that node $i$ sends to help in the repair of node $j$ when ${\cal D}$ is the set of helper nodes. (Wherever clear, we ignore the ${\cal D}$ in the notation and simply write $S_i^j$.) When node $1$ fails, the information we therefore have at its replacement node can be written as:
\begin{eqnarray}\label{eq:repl-node-1}
\begin{blockarray}{cccccccc}
{\color{gray}{\bf x}_1} & {\color{gray}{\bf x}_2} & {\color{gray}\hspace{-4ex}\cdots} &\hspace{-2ex} {\color{gray}{\bf x}_h} & {\color{gray}\hspace{-2ex}{\bf x}_{h+1}} & {\color{gray}\hspace{-2ex}\cdots} & {\color{gray}\hspace{-2ex}{\bf x}_{k-1}} & {\color{gray}\hspace{-4ex}{\bf x}_k}\\
    \begin{block}{(cccccccc)}
& & & \hspace{-2ex}S_h^1 & & & &\\
& & & & \hspace{-2ex}S_{h+1}^1 & & &\\
& & & & & \hspace{-2ex}\ddots & &\\
& & & & & & \hspace{-2ex}S_{k-1}^1&\\
& & & & & & &\hspace{-4ex}S_k^1\\
S_{k+1}^1A_{11}&\hspace{-1.2ex}S_{k+1}^1A_{12}&&\hspace{-2ex}\cdots&&\hspace{-2ex}\cdots&&\hspace{-4ex}S_{k+1}^1A_{1k}\\
\vdots&\vdots&&\hspace{-2ex}\vdots&&\hspace{-2ex}\vdots&&\hspace{-4ex}\vdots\\
S_{k+h}^1A_{h1}&\hspace{-1.2ex}S_{k+h}^1A_{h2}&&\hspace{-2ex}\cdots&&\hspace{-2ex}\cdots&&\hspace{-4ex}S_{k+h}^1A_{hk}\\
    \end{block}
  \end{blockarray}
\hspace{-.7ex}\left( \hspace{-1.3ex} \begin{array}{c}
  {\bf x}_1 \\ {\bf x}_2  \\\vdots\\{\bf x}_h \\ {\bf x}_{h+1}\\\vdots \\ {\bf x}_{k-1}\\{\bf x}_k
  \end{array}\hspace{-1.8ex}\right)\hspace{-1.4ex}
\end{eqnarray}

Suppose all $S_j^1$'s in (\ref{eq:repl-node-1}) be replaced by a repair subspace $S^1$ (corresponding to Lemma \ref{lemma:interference alignment_n-1}(b)) that we would have used if ${\cal D} = \{2, 3,\ldots, k, k+1, k+2\}$. Specifically, suppose $S^1A_{1,1}$ and $S^1$ complete the space $\mathbb{F}_q^{\alpha}$.
Since $S^1$ and $S^1A_{i,j}$ denote the same subspace,
for $j \neq 1$, the components of ${\bf x}_i, i \in \{h, h+1,\ldots,k\}$ can be easily subtracted from the information coming from the parity nodes
, using that coming from the systematic nodes $h$ to $k$. Thus, in order to recover ${\bf x}_1$, we can concentrate on the following information at the replacement node:
\begin{eqnarray}\label{eq:repl-node-1b}
\begin{blockarray}{ccccc}
{\color{gray}{\bf x}_1} & {\color{gray}{\bf x}_2} & {\color{gray}\cdots} & {\color{gray}{\bf x}_{h-2}} & {\color{gray}{\bf x}_{h-1}}\\
 \begin{block}{(ccccc)}
S^1A_{11}&S^1A_{12}&\cdots&S^1A_{1,h-2}&S^1A_{1,h-1}\\
S^1A_{21}&S^1A_{22}&\cdots&S^1A_{2,h-2}&S^1A_{2,h-1}\\
\vdots&\vdots&\vdots&\vdots&\vdots\\
S^1A_{h1}&S^1A_{h2}&\cdots&S^1A_{h,h-2}&S^1A_{h,h-1}\\
    \end{block}
  \end{blockarray}
\hspace{-.5ex}\left( \hspace{-1.2ex} \begin{array}{c}
  {\bf x}_1 \\ {\bf x}_2  \\\vdots\\{\bf x}_{h-2} \\ {\bf x}_{h-1}
  \end{array}\hspace{-1ex}\right)\hspace{-2ex}
\end{eqnarray}

Let $S^1A_{i,j} = \lambda_{i,j}S^1Q_{i,j} = \lambda_{i,j}\widetilde{Q}_{i,j}S^1$, where $\widetilde{Q}_{i,j}$ is an $\alpha/2 \times \alpha/2$ matrix, and $(i,j) \neq (1,1)$. It must be noted that not only is $S^1$ dependent on the choice of ${\cal D}$, but so in turn is $\widetilde{Q}_{i,j}$. Let us also denote $S^1{\bf x}_i$ by $\widetilde{{\bf x}}_i$. Then, (\ref{eq:repl-node-1b}) can be rewritten as:
\begin{eqnarray}\label{eq:repl-node-1c}
\left(\begin{array}{cccc}
\lambda_{1,1}S^1Q_{1,1}&\hspace{-1.2ex}\lambda_{1,2}\widetilde{Q}_{1,2}&\hspace{0ex}\cdots&\hspace{-.8ex}\lambda_{1,h-1}\widetilde{Q}_{1,h-1}\\
\lambda_{2,1}\widetilde{Q}_{2,1}S^1&\hspace{-1.2ex}\lambda_{2,2}\widetilde{Q}_{2,2}&\hspace{0ex}\cdots&\hspace{-.8ex}\lambda_{2,h-1}\widetilde{Q}_{2,h-1}\\
\vdots&\hspace{-1.2ex}\vdots&\hspace{0ex}\ddots&\hspace{-.8ex}\vdots\\
\lambda_{h,1}\widetilde{Q}_{h,1}S^1&\hspace{-1.2ex}\lambda_{h,2}\widetilde{Q}_{h,2}&\hspace{0ex}\cdots&\hspace{-.8ex}\lambda_{h,h-1}\widetilde{Q}_{h,h-1}\\
\end{array}\right)
\hspace{-.5ex}\left( \hspace{-1.1ex} \begin{array}{c}
  {\bf x}_1 \\ \widetilde{{\bf x}}_2  \\\vdots\\ \widetilde{{\bf x}}_{h-1}
  \end{array}\hspace{-1.2ex}\right).\hspace{-2ex}
\end{eqnarray}

The matrix in (\ref{eq:repl-node-1c}) --- call it $M$ --- is a square matrix of dimensions $h\alpha/2 \times h\alpha/2$. A sufficient condition to recover ${\bf x}_1$ is that $M$ is invertible. Notice that the determinant of $M$, ${\sf det}(M)$, is a polynomial in the following variables: $\lambda_{i,j}, i \in [h], j \in [h-1]$.
Hence, ${\sf det}(M)$ is a nonzero polynomial of degree $h\alpha/2$ in the given variables.
From Schwartz--Zippel--DeMillo--Lipton lemma,
 if the finite field $\mathbb{F}_q$ over which the determinant is defined has cardinality $|\mathbb{F}_q| = q > h\alpha/2$, there exist $\lambda_{i,j}$'s for which the determinant ${\sf det}(M)$ above is nonzero.

\smallskip
\emph{ Part 2}: Notice that $M$ above is defined for a particular example scenario. In general, let the number of helper nodes be $d$, the failed systematic node be $f \in [k]$, the set of helper nodes by ${\cal D} \subseteq [n]\backslash \{f\}$, the set of systematic helper nodes be ${\cal D}_s \subseteq [k]\backslash \{f\}$, and the set of parity helper nodes be ${\cal D}_p \subseteq [k+1:k+r]$. Let the number of parity helper nodes be denoted by $h$, where $h$ ranges from $d-k+1$ to $r$. Let us represent by ${\cal H}_p$ the set of parity helper nodes but indexed within $[r]$, where $i$ corresponds to node $k+i$ of the system, that is, ${\cal H}_p = \{i \,|\, k+i \in {\cal D}_p\} \subseteq [r]$.

The matrix $M$ in (\ref{eq:repl-node-1c}), in general, can be seen to be a square matrix of dimensions $h\alpha/(d-k+1) \times h\alpha/(d-k+1)$. In particular, $M$ is a function of $f$, ${\cal D}_s$, and ${\cal H}_p$, and the determinant polynomial has degree which is a function of $|{\cal D}_p| = h$ and $d$. For each $f$, ${\cal D}_s$ and ${\cal H}_p$, we obtain a sufficiency condition that the corresponding $M$ is invertible. Therefore, the product of the corresponding determinant polynomials is a nonzero polynomial of degree
\vspace{0ex}\begin{eqnarray*}
q_{\text{ANY}} &=& k\left(\sum_{h=d-k+1}^r {r \choose h}{k-1 \choose d-h} \frac{h\alpha}{d-k+1}\right)\\
&=& \left(\sum_{h=d-k+1}^r h{r \choose h}{k-1 \choose d-h}\right)\frac{k\alpha}{d-k+1};
\end{eqnarray*}
\vspace{.1ex}consequently, there exist $\lambda_{i,j}$'s in $\mathbb{F}$ such that any systematic node is repairable with optimal repair bandwidth using any arbitrary set of $d$ helper nodes, as long as the field size $|\mathbb{F}| > q_{\text{ANY}}$.

\emph{ Part 3:} Finally, using Lemma \ref{lemma:MDS property1}, Lemma \ref{lemma:MDS property2}, and Lemma \ref{thm:fixed_d}, we obtain an $[n,k,d]$ systematic--repair MSR code for any set of $d$ helper nodes, when the field size $q > q_{\text{ANY}} + q_{\text{MDS}}$.
\end{IEEEproof}

\section{Conclusion}\label{sec:conclusion}
In this paper we presented a new construction for systematic--repair MSR codes for all possible values of parameters $[n,k,d]$.

A more generalized construction, where a single $[n,k]$ code simultaneously satisfies the optimal repair for all $d\in\{k+1,\cdots,n-1\}$ will be introduced in a sequel paper. It is to be noted that both these generalizations come at the cost of increasing $\alpha$. A lower bound on $\alpha$ is proved in \cite{GTC14} when $d=n-1$. Whether similar bounds exist for general $[n,k,d]$ or not is left for  future work. So is the question of constructing MSR codes that also optimally repair parity nodes.

\bibliographystyle{IEEEtran}