\section{Introduction}

Quantum Inverse Scattering Method (QISM) is a powerful tool for solving quantum integrable models \cite{TakF79,SklTF79,KorBI93L,FadLH96}. This method
allows one to find spectra of quantum Hamiltonians  via the algebraic Bethe ansatz. The main advantage of the
algebraic Bethe ansatz is that it gives a systematic procedure to describe the spectra of the models, which might have
completely different physical interpretation. This is because this method only deals with the algebra of the monodromy
matrix entries, but not with its specific representation.

The QISM and the algebraic Bethe ansatz also can be used for calculation of form factors and correlation functions \cite{IzeK84,KorBI93L,KitMST05,Sla07}.
Similarly to the problem
of the Hamiltonian spectrum, in many cases this method gives quite general results, which can be used for the study of a wide class of models.
In particular, it was shown recently \cite{PakRS15d,PakRS15b} that form factors of local operators (FFLO) in the models with $\mathfrak{gl}(3)$-invariant $R$-matrix are all expressed in terms of {\it universal form factors} \cite{PakRS15a}. The latter are completely determined by
the $R$-matrix and do not depend on the model under consideration.

Knowing FFLO one can solve the problem of correlation functions via their form factor expansion. In the models,
for which an explicit solution of the quantum inverse scattering problem is known \cite{KitMT99,MaiT00}, the FFLO
are directly related to the ones of the monodromy matrix entries. Using this result, correlation functions of the $XXZ$ spin chain
and other integrable models were studied in the
series of works \cite{KitMST05,KitKMST11,KitKMST12,CauM05,CauPS07}. It is worth mentioning, however, that the existence of an  explicit solution of the
quantum inverse scattering problem is based on a specific representation of the underlying $RTT$-algebra. It can be found for various spin chains, but not
in general.

In the present paper we consider models described by the $\mathfrak{gl}(2|1)$ and $\mathfrak{gl}(1|2)$ superalgebras. Determinant representations for form factors of the monodromy matrix entries in these models were found in \cite{HutLPRS16a}. These formulas can be directly applied to the calculation of
the FFLO in the  supersymmetric t-J model \cite{Sch87,For89,EssK92,FoeK93,Gom02}, because  the quantum
inverse scattering problem for this model was solved in \cite{GohK00}.
Our goal is to generalize these results. Namely, we show that similarly to the $\mathfrak{gl}(N)$ case, the FFLO in the
models with $\mathfrak{gl}(2|1)$ and $\mathfrak{gl}(1|2)$ supersymmetries are proportional to the universal form factors.

Our method is based on the composite model \cite{IzeK84}. It was used in \cite{PakRS15b,PakRS15c} for the calculation of the FFLO
in the $\mathfrak{gl}(3)$-based models. In the present paper we generalize this approach to the case of superalgebras. Due to an isomorphism
between Yangians of $\mathfrak{gl}(2|1)$ and $\mathfrak{gl}(1|2)$ \cite{PakRS17} we consider below the case of $Y\bigl(\mathfrak{gl}(2|1)\bigr)$ only.
It is worth mentioning that our method does not use any specific representation of the $RTT$-algebra. Therefore, the obtained results have the same degree of generality as the Bethe equations which are used to determine the spectra of the quantum Hamiltonians.

The paper is organized as follows. In section~\ref{S-BNN} we introduce the model under consideration and describe the basic
notions and the notation used in the paper. In section~\ref{SS-CM} we recall the notion of the composite model and introduce partial
zero modes of the monodromy matrix. There we explain how the use of the partial
zero modes allows one to relate different FFLO with each other. The main results of the paper are collected in
section~\ref{S-MR}. In section~\ref{S-FFOD} we prove the results of section~\ref{S-MR}. A part of the proof is placed in appendix~\ref{A-CALC}.

\section{Basic notions and notation\label{S-BNN}}

In this section we briefly describe basic notions of $\mathfrak{gl}(2|1)$-based integrable models solvable by the algebraic Bethe ansatz.
The reader can find a more detailed description in \cite{KulS80,PakRS17,BeRa08,HutLPRS17a}.

\subsection{Graded models and Bethe vectors\label{SS-GMBV}}

The central object of the algebraic Bethe ansatz method is a quantum monodromy matrix $T(u)$. It
acts in the tensor product  $V\otimes \mathcal{H}$, where $\mathcal{H}$ is a Hilbert space of a quantum model and
$V$ is an auxiliary space. Commutation relations between the entries $T_{ij} (u)$ are gathered in
an $RTT$-relation
\begin{equation}\label{L-op-com}
R(u,v)\cdot (T(u)\otimes \mathbb{I})\cdot (\mathbb{I}\otimes T(v))=
(\mathbb{I}\otimes T(v))\cdot (T(u)\otimes \mathbb{I})\cdot R(u,v).
\end{equation}
Here $R(u,v)$ is an $R$-matrix acting in the tensor product $V\otimes V$ of the auxiliary vector spaces $V$.  Equation \eqref{L-op-com} holds in the tensor
product $V\otimes V\otimes\mathcal{H}$.

For $\mathfrak{gl}(2|1)$-based models the auxiliary vector space $V$ is a ${\mathbb Z}_2$-graded space ${\mathbb C}^{2|1}$  with a basis $\{{\rm e}_1,
{\rm e}_2,{\rm e}_3\}$. We call the vectors  $\{{\rm e}_1,{\rm e}_2\}$  even, while ${\rm e}_{3}$ is odd. Respectively, we  introduce  a parity function on the set of indices as  $[1]=[2]=0$  and $[3]=1$.

The $R$-matrix in \eqref{L-op-com} has the form
\begin{equation}\label{DYglmn}
R(u,v)\ =\  \mathbb{I}+g(u,v)\mathbb{P}, \qquad g(u,v)=\frac{c}{u-v},
\end{equation}
where $c$ is a constant, $\mathbb{I}$ is the identity matrix in $V\otimes V$, and $\mathbb{P}$ is the graded
permutation matrix \cite{KulS80}. The tensor product in  \eqref{L-op-com}
 is also graded leading to the set of commutation relations between the monodromy matrix entries $T_{ij}$:
\begin{equation}\label{TM-1}
[T_{ij}(u),T_{kl}(v)\}
=(-1)^{[i]([k]+[l])+[k][l]}g(u,v)\Big(T_{kj}(v)T_{il}(u)-T_{kj}(u)T_{il}(v)\Big),
\end{equation}
where we have introduced a graded commutator as
\be{Def-SupC}
[T_{ij}(u),T_{kl}(v)\}= T_{ij}(u)T_{kl}(v) -(-1)^{([i]+[j])([k]+[l])}   T_{kl}(v)  T_{ij}(u).
\end{equation}
The supertrace of the monodromy matrix
 %
 \be{Strace}
 \mathcal{T}(u)=\mathop{\rm str} T(u)=\sum_{i=1}^3 (-1)^{[i]}T_{ii}(u)
 \end{equation}
 %
 is called the transfer matrix. It is a generating function of the integrals of motion of the integrable model under consideration.
 The transfer matrix eigenstates are called {\it on-shell Bethe vectors}. They play an important role in the considerations below.

We assume that the space $\mathcal{H}$, in which the operators $T_{ij}$ act, contains a pseudovacuum vector $|0\rangle$. This vector possesses the following
properties:
\be{vac}
\begin{aligned}
T_{ii}(u)|0\rangle&=\lambda_i(u)|0\rangle,\\
T_{ij}(u)|0\rangle&=0,\qquad i>j,
\end{aligned}
\end{equation}
where $\lambda_i(u)$ are some functions of complex variable $u$. A specific choice of these functions means fixing of a specific
integrable model. For us, however, they remain free functional parameters. This treatment of $\lambda_i(u)$ allows us
to consider a wide class of integrable models within a common framework.

We also assume that the operators $T_{ij}$  act in the dual space $\mathcal{H}^*$ with a dual pseudovacuum vector $\langle0|$. This vector
has analogous properties
\be{dvac}
\begin{aligned}
\langle0|T_{ii}(u)&=\lambda_i(u)\langle0|,\\
\langle0|T_{ij}(u)&=0,\qquad i<j,
\end{aligned}
\end{equation}
where $\lambda_i(u)$ are the same as in \eqref{vac}. Below we use the ratios of these functions
 \be{ratios}
 r_1(w)=\frac{\lambda_1(w)}{\lambda_2(w)}, \qquad  r_3(w)=\frac{\lambda_3(w)}{\lambda_2(w)}.
 \end{equation}

Bethe vectors of $\mathfrak{gl}(2|1)$-invariant models are certain polynomials in operators $T_{ij}(u)$ with $i<j$ acting on the pseudovacuum vector. Their
explicit form was found in \cite{PakRS17} (see also \cite{HutLPRS17a} for the general $\mathfrak{gl}(m|n)$ case). They depend on two sets of variables called Bethe parameters.
We denote the Bethe vectors $\mathbb{B}_{a,b}(\bar u;\bar v)$. Here the Bethe parameters are $\bar u=\{u_1,\dots,u_a\}$ and
$\bar v=\{v_1,\dots,v_b\}$. The subscripts $a$ and $b$ ($a,b=0,1,\dots$) denote the cardinalities of the sets $\bar u$ and $\bar v$ respectively.

Similarly one can construct dual Bethe vectors in the dual space $\mathcal{H}^*$ as polynomials in operators $T_{ij}(u)$ with $i>j$ acting on the
dual pseudovacuum vector $\langle0|$ \cite{PakRS17}. We denote them $\mathbb{C}_{a,b}(\bar u;\bar v)$ with the same meaning of the arguments and
the subscripts.

For generic (dual) Bethe vectors the Bethe parameters $\bar u$ and $\bar v$ are generic complex numbers. If these parameters satisfy a system of Bethe equations
\be{BE}
\begin{aligned}
r_1(u_j)&=\prod_{\substack{k=1\\k\ne j}}^a\frac{f(u_j,u_k)}{f(u_k,u_j)}\prod_{l=1}^bf(v_l,u_j),\qquad j=1,\dots,a,\\
r_3(v_j)&=\prod_{l=1}^af(v_j,u_l),\qquad j=1,\dots,b,
\end{aligned}
\end{equation}
where
\be{f}
f(u,v)=1+g(u,v)=\frac{u-v+c}{u-v},
\end{equation}
then the corresponding (dual) Bethe vector becomes on-shell. It means that this is an eigenvector of the transfer matrix
$ \mathcal{T}(u)$ \eqref{Strace}

\be{Left-act}
\mathcal{T}(w)\mathbb{B}_{a,b}(\bar u;\bar v)= \tau(w|\bar u,\bar v)\,\mathbb{B}_{a,b}(\bar u;\bar v),\qquad
 \mathbb{C}_{a,b}(\bar u;\bar v)\mathcal{T}(w) = \tau(w|\bar u,\bar v)\,\mathbb{C}_{a,b}(\bar u;\bar v),
\end{equation}
where the eigenvalue $\tau(w|\bar u,\bar v)$ is
\be{tau-def}
\tau(w|\bar u,\bar v)=\lambda_1(w)\prod_{j=1}^af(u_j,w)+ \lambda_2(w)\prod_{j=1}^af(w,u_j)\prod_{k=1}^b f(v_k,w)
 -\lambda_3(w)\prod_{k=1}^bf(v_k,w).
\end{equation}

Besides the monodromy matrix $T(u)$ we also consider a twisted monodromy matrix $T_\kappa(u)=\kappa T(u)$, where $\kappa$ is a
diagonal matrix $\kappa=\mathop{\rm diag}\{\kappa_1,\kappa_2,\kappa_3\}$, and $\kappa_i$ are complex numbers. The supertrace
$\mathcal{T}_\kappa(u)=\mathop{\rm str} T_\kappa(u)$ is called the twisted transfer matrix. A generic Bethe vector becomes an eigenstate of the
twisted transfer matrix, if the Bethe parameters satisfy a system of twisted Bethe equations
\be{TBE}
\begin{aligned}
r_1(u_j)&=\frac{\kappa_2}{\kappa_1}\prod_{\substack{k=1\\k\ne j}}^a\frac{f(u_j,u_k)}{f(u_k,u_j)}\prod_{l=1}^bf(v_l,u_j),\qquad j=1,\dots,a,\\
r_3(v_j)&=\frac{\kappa_2}{\kappa_3}\prod_{l=1}^af(v_j,u_l),\qquad j=1,\dots,b.
\end{aligned}
\end{equation}
The corresponding (dual) Bethe vector is then called the twisted (dual)  on-shell Bethe vector. The twisted transfer matrix eigenvalue
on the vector $\mathbb{B}_{a,b}(\bar u;\bar v)$ is given by
\be{ktau-def}
\tau_\kappa(w|\bar u,\bar v)=\kappa_1\lambda_1(w)\prod_{j=1}^af(u_j,w)+ \kappa_2\lambda_2(w)\prod_{j=1}^af(w,u_j)\prod_{k=1}^b f(v_k,w)
 -\kappa_3\lambda_3(w)\prod_{k=1}^bf(v_k,w).
\end{equation}
The use of the twisted monodromy matrix and the twisted on-shell Bethe vectors allows us to construct a special generating functional for
FFLO (see section~\ref{S-FFOD}).

\subsection{Shorthand notation\label{SS-SN}}

We denote sets of variables by bar:  $\bar u$, $\bar v$ etc. If necessary, the cardinalities of the sets are
given in special comments.
Individual elements of the sets are denoted by latin subscripts: $u_j$, $v_k$ etc.
We say that $\bar x=\bar x'$,
if $\#\bar x=\#\bar x'$ and $x_i=x'_i$ (up to a permutation) for $i=1,\dots,\#\bar x$. We say that $\bar x\ne \bar x'$ otherwise.

Below we consider partitions of the sets into subsets.  The notation $\bar u\Rightarrow\{\bar u_{\scriptscriptstyle \rm I},\bar u_{\scriptscriptstyle \rm I\hspace{-1pt}I}\}$ means that the set $\bar u$
is divided into two disjoint subsets. As a rule, we use roman numbers for subscripts of subsets: $\bar u_{\scriptscriptstyle \rm I}$, $\bar v_{\rm ii}$ etc.
However, if we deal with a big quantity of subsets, then we use standard arabic numbers for their notation. In such cases we give special comments  to avoid ambiguities.

To lighten long formulas we use a shorthand notation for products of  some functions.
Namely, if the functions $r_k$ \eqref{ratios} or the functions $g$ and $f$  depend
on sets of variables, this means that one should take the product over the corresponding set.
For example,
 \be{SH-prod}
 r_1(\bar u)=\prod_{u_k\in\bar u} r_1(u_k);\quad
 g(z, \bar w)= \prod_{w_j\in\bar w} g(z, w_j);\quad
 f(\bar u,\bar v)=\prod_{u_j\in\bar u}\prod_{v_k\in\bar v} f(u_j,v_k).
 \end{equation}
By definition any product with respect to the empty set is equal to $1$. If we have a double product, then it is also equal to $1$ if at least
one of the sets is empty.

In section~\ref{SS-CM} we shall introduce several new scalar functions and will extend
the convention \eqref{SH-prod} to their products.

\subsection{Universal form factors\label{SS-UF}}

In this paper we reduce FFLO to the universal form factors of the monodromy matrix entries. The latter  are defined as
follows
\be{Univ-FF}
\mathfrak{F}^{(i,j)}\left(\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a',a}_{b',b}=
\frac{ \mathbb{C}_{a',b'}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})T_{ij}(z)\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})}
{\tau(z|\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C})-\tau(z|\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B})}.
\end{equation}
Here both $\mathbb{C}_{a',b'}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})$ and $\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$ are on-shell
Bethe vectors, and we assume that the Bethe parameters of the two Bethe vectors are different: $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\}\ne \{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$.
The parameter $z$ is an arbitrary complex  number. It was proved in \cite{HutLPRS16a} that the ratio
in the r.h.s. of \eqref{Univ-FF} does not depend on $z$.

The form factors \eqref{Univ-FF} are called universal, because they are completely determined by the $R$-matrix of the model. They do not depend
on a specific representation of the $RTT$-algebra, in particular, they do not depend on the vacuum eigenvalues $\lambda_i(u)$ \eqref{vac}. In other words,
if two different integrable models are described by the $R$-matrix \eqref{DYglmn}, then they have the same universal form factors. Explicit determinant
representations for the universal form factors in the $\mathfrak{gl}(2|1)$-invariant models were obtained in \cite{HutLPRS16a}.

\section{Composite model\label{SS-CM}}

In order to access FFLO we introduce a composite model \cite{IzeK84,PakRS15b,Fuk17}.
Most naturally the composite model arises in the lattice models, where the  monodromy matrix $T(u)$  is equal to the product of local $L$-operators
\be{mat-T}
T(u)=L_M(u)\cdots L_1(u).
\end{equation}
Here $M$ is the number of the lattice sites, and every $L$-operator satisfies $RTT$-relation with $R$-matrix \eqref{DYglmn}. Let us fix a
site $m$ ($1\le m<M$) and define two partial monodromy matrices $T^{(1)}(u)$  and $T^{(2)}(u)$ as
\be{mat-T12}
T^{(1)}(u)=L_m(u)\cdots L_1(u),\qquad T^{(2)}(u)=L_M(u)\cdots L_{m+1}(u).
\end{equation}
Then
\be{T-TT}
T(u)=T^{(2)}(u)T^{(1)}(u).
\end{equation}
Every $T^{(l)}(u)$ obviously satisfies $RTT$-relation \eqref{L-op-com} and has its own pseudovacuum vector $|0\rangle^{(l)}$, such that $|0\rangle= |0\rangle^{(1)}\otimes|0\rangle^{(2)}$. The operators $T_{ij}^{(2)}(u)$ and $ T_{kl}^{(1)}(v)$ supercommute with each other, as they act in different spaces.

Continuous quantum models can be obtained from the lattice ones in the limit $M\to\infty$. Obviously, the determining relation \eqref{T-TT} remains
unchanged. The partial monodromy matrices still satisfy the $RTT$ relation, and the entries of the different partial  monodromy matrices mutually supercommute.
Thus, continuous quantum models also can be considered in the framework of the composite model.

Let
\be{eigen}
T_{ii}^{(l)}(u)|0\rangle^{(l)}= \lambda_{i}^{(l)}(u)|0\rangle^{(l)}, \qquad l=1,2,
\end{equation}
where $\lambda_{i}^{(l)}(u)$ are new free functional parameters. We also introduce
\be{rk}
r_{k}^{(l)}(u)=\frac{\lambda_{k}^{(l)}(u)}{\lambda_{2}^{(l)}(u)}, \qquad l=1,2, \qquad k=1,3.
\end{equation}
Obviously
\be{lr}
\lambda_{i}(u)=\lambda_{i}^{(1)}(u)\lambda_{i}^{(2)}(u), \qquad
r_{k}(u)=r_{k}^{(1)}(u)r_{k}^{(2)}(u).
\end{equation}
Below we  express form factors in terms of $r_{k}^{(1)}(u)$, therefore we introduce a special notation for
these functions
\be{ell}
r_{k}^{(1)}(u)=\ell_k(u),\qquad \text{and hence,}\qquad r_{k}^{(2)}(u)=\frac{r_k(u)}{\ell_k(u)}, \qquad k=1,3.
\end{equation}
We extend the convention on the shorthand notation \eqref{SH-prod} to the products of the functions $\ell_k(u)$.

Any monodromy matrix in \eqref{T-TT} possesses its own Bethe vectors.
The total Bethe vector is a bilinear combination of the partial Bethe vectors \cite{Fuk17}:
\begin{equation} \label{BVcomp}
\mathbb{B}_{a,b}(\bar u;\bar v) = \sum r_1^{(2)}(\bar u_{\scriptscriptstyle \rm I}) r_3^{(1)}(\bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I}) \frac{f(\bar u_{\scriptscriptstyle \rm I\hspace{-1pt}I}, \bar u_{\scriptscriptstyle \rm I})g(\bar v_{\scriptscriptstyle \rm I},\bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I})}{f(\bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I},\bar u_{\scriptscriptstyle \rm I})} \mathbb{B}_{a_2,b_2}^{(2)}(\bar u_{\scriptscriptstyle \rm I\hspace{-1pt}I}; \bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I}) \mathbb{B}^{(1)}_{a_1,b_1}(\bar u_{\scriptscriptstyle \rm I}; \bar v_{\scriptscriptstyle \rm I}).
\end{equation}
Here the sum is taken  over all partitions $\bar u \Rightarrow \{\bar u_{\scriptscriptstyle \rm I}, \bar u_{\scriptscriptstyle \rm I\hspace{-1pt}I}\}$ and $\bar v \Rightarrow \{\bar v_{\scriptscriptstyle \rm I}, \bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I}\}$. The
cardinalities of the subsets satisfy $a_1+a_2=a$ and $b_1+b_2=b$. Recall that here we have used the convention \eqref{SH-prod} for the products
of the functions $r_k^{(\ell)}$, $f$, and $g$.

Similarly, the dual total Bethe vector is a bilinear combination of the partial dual Bethe vectors:
\begin{align} \label{dBVcomp}
&\mathbb{C}_{a,b}(\bar u;\bar v )= \sum r_1^{(1)}(\bar u_{\scriptscriptstyle \rm I\hspace{-1pt}I}) r_3^{(2)}(\bar v_{\scriptscriptstyle \rm I}) \frac{f(\bar u_{\scriptscriptstyle \rm I}, \bar u_{\scriptscriptstyle \rm I\hspace{-1pt}I}) g(\bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I},\bar v_{\scriptscriptstyle \rm I})}{f(\bar v_{\scriptscriptstyle \rm I},\bar u_{\scriptscriptstyle \rm I\hspace{-1pt}I})} \mathbb{C}_{a_1,b_1}^{(1)}(\bar u_{\scriptscriptstyle \rm I}; \bar v_{\scriptscriptstyle \rm I}) \mathbb{C}^{(2)}_{a_2,b_2}(\bar u_{\scriptscriptstyle \rm I\hspace{-1pt}I}; \bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I}).
\end{align}
Here the sum is the same as in \eqref{BVcomp}.

Observe that if the total (dual) Bethe vector is on-shell (i.e. the set $\{\bar u,\bar v\}$ satisfies Bethe equations), then the partial (dual) Bethe vectors
generically are not on-shell, because the subsets $\{\bar u_{\scriptscriptstyle \rm I},\bar v_{\scriptscriptstyle \rm I}\}$ and $\{\bar u_{\scriptscriptstyle \rm I\hspace{-1pt}I},\bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I}\}$ do not satisfy Bethe equations.

Comparing these formulas with the formulas for the total (dual) Bethe vectors in $\mathfrak{gl}(3)$-based models \cite{PakRS15d} one can see that
the difference is very small. Namely, replacing the product of functions $g(\bar v_{\scriptscriptstyle \rm I},\bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I})$ in \eqref{BVcomp} with  the product $f(\bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I},\bar v_{\scriptscriptstyle \rm I})$
we obtain the expression for the total Bethe vector in the models with $\mathfrak{gl}(3)$ symmetry. Similarly, for the dual Bethe vectors
one should make the replacement $g(\bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I},\bar v_{\scriptscriptstyle \rm I})\;\to\;f(\bar v_{\scriptscriptstyle \rm I},\bar v_{\scriptscriptstyle \rm I\hspace{-1pt}I})$.
This similarity makes it possible to calculate FFLO by the
same methods as in the case of the $\mathfrak{gl}(3)$-based models \cite{PakRS15b}.

\subsection{Zero modes\label{SS-ZM}}

We assume a standard representation of the local $L$-operators in \eqref{mat-T}:
\be{L-dep-u}
L_n(u)=\mathbf{1}+\frac cu L_n[0] +o(u^{-1}),\qquad u\to\infty.
\end{equation}
Here $\mathbf{1}$ is the identity operator in $\mathbb{C}^{2|1}\otimes\mathcal{H}$. The matrix elements
$\bigl(L_n[0]\bigr)_{ij}$ depend on the local operators of the model. Due to \eqref{mat-T} we conclude that
the total monodromy matrix $T(u)$ and both partial monodromy matrices $T^{(l)}(u)$ have the standard expansion over $c/u$:
\be{zero-modes}
\begin{aligned}
&T(u)=\mathbf{1}+ \sum_{n=0}^\infty T[n]\,\left(\frac cu\right)^{n+1},\\
&T^{(l)}(u)=\mathbf{1}+ \sum_{n=0}^\infty T^{(l)}[n]\,\left(\frac cu\right)^{n+1},\qquad l=1,2.
\end{aligned}
\end{equation}
The operators $T[0]$ and $T^{(l)}[0]$ respectively are called the total and the partial zero modes. Obviously, the partial zero mode $T^{(1)}[0]$ is equal to
\be{part-zm}
T^{(1)}[0]=\sum_{n=1}^m L_n[0].
\end{equation}

Consider a form factor of the partial zero mode $T^{(1)}_{ij}[0]$
\be{FF-partZM}
\mathcal{M}^{(i,j)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a',a}_{b',b}=\mathbb{C}_{a',b'}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})T^{(1)}_{ij}[0]\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B}),
\end{equation}
where  $\mathbb{C}_{a',b'}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})$ and $\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$ are on-shell
Bethe vectors. We have stressed that this form factor depends on the number  $m$ of the bulk site in \eqref{mat-T12}.
Then due to \eqref{part-zm} we obtain
\be{FF-partL}
\mathbb{C}_{a',b'}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})\bigl(L_m[0]\bigr)_{ij}\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})=\mathcal{M}^{(i,j)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a',a}_{b',b}- \mathcal{M}^{(i,j)}\left(m-1\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a',a}_{b',b}.
\end{equation}
Thus, knowing the form factors of the partial zero mode $T^{(1)}[0]$ we can find the form factors of the local operators
$\bigl(L_m[0]\bigr)_{ij}$. It is clear that
in the case of continuous models the finite difference in the r.h.s. of \eqref{FF-partL} turns into the derivative over a space variable.

The use of the zero modes also allows one to obtain simple relations between different form factors. It follows from the commutation relations
\eqref{TM-1} that
\be{Comm-ZM}
\begin{aligned}
 {}[T_{ij}[0],T_{kl}[0]\} & = (-1)^{[i][j]+[i][l]+[j][l]}(\delta_{il} T_{kj}[0] -\delta_{kj} T_{il}[0]),\\
[T_{ij}^{(s)}[0],T_{kl}^{(s)}[0]\} & = (-1)^{[i][j]+[i][l]+[j][l]}(\delta_{il} T_{kj}^{(s)}[0] -\delta_{kj} T_{il}^{(s)}[0]),
\qquad s=1,2.\\
\end{aligned}
\end{equation}
Using $T_{ij}[0]=T_{ij}^{(1)}[0]+T_{ij}^{(2)}[0]$ and $[T_{ij}^{(1)}[0],T_{kl}^{(2)}[0]\} =0$ we arrive at
\be{Comm-ZM12}
[T_{ij}^{(1)}[0],T_{kl}[0]\}  = (-1)^{[i][j]+[i][l]+[j][l]}(\delta_{il} T_{kj}^{(1)}[0] -\delta_{kj} T_{il}^{(1)}[0]).
\end{equation}
Equation \eqref{Comm-ZM12} yields
\begin{multline}\label{FFlk-FFij}
\delta_{il} \mathcal{M}^{(k,j)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a',a}_{b',b}- \delta_{kj} \mathcal{M}^{(i,l)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a',a}_{b',b}\\
=(-1)^{[i][j]+[i][l]+[j][l]}\mathbb{C}_{a',b'}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})[T_{ij}^{(1)}[0],T_{kl}[0]\}\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B}).
\end{multline}
The actions of the total zero modes $T_{kl}[0]$ onto (dual) on-shell Bethe vectors were studied in \cite{HutLPRS16a}. Under
this action a (dual) on-shell vector either vanishes or remains on-shell. Therefore, the expectation value  in the r.h.s.
of \eqref{FFlk-FFij} is related to the form factor of the partial zero mode $T_{ij}^{(1)}[0]$. Equation \eqref{FFlk-FFij}, thus,
allows us to express this form factor in terms of $\mathcal{M}^{(k,j)}$ and $\mathcal{M}^{(i,l)}$. We consider specific examples
of these relationships in section~\ref{S-FFOD}.

\section{Main results\label{S-MR}}

We have shown in the previous section that FFLO can be reduced to the form factors of the partial zero
modes $T^{(1)}_{ij}[0]$. Studying these form factors one should distinguish between two cases. In the first case an on-shell Bethe vector
$\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$ and a dual on shell Bethe vector $\mathbb{C}_{a',b'}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})$ correspond to different eigenvalues of the transfer
matrix. Then we say that $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\} \neq \{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$. Overwise, if $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\} =\{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$, then both vectors
correspond to the same eigenvalue. The latter case occurs for form factors of the diagonal elements $T^{(1)}_{ii}[0]$ only.
\begin{thm}
Let $\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$ and $\mathbb{C}_{a',b'}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})$ be on-shell (dual) Bethe vectors such that $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\} \neq \{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$. Then
\begin{equation}\label{FF-gener}
  \mathcal{M}^{(i,j)}
\left( m\Bigr|\begin{smallmatrix}
  \bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a',a}_{b',b}

  = \left( \frac{\ell_1(\bar u^{\scriptscriptstyle C})\ell_3(\bar v^{\scriptscriptstyle B})}{\ell_1(\bar u^{\scriptscriptstyle B})\ell_3(\bar v^{\scriptscriptstyle C})} -1 \right) \mathfrak{F}^{(i,j)}
\left( \begin{smallmatrix}
  \bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a',a}_{b',b},
\end{equation}
where $\mathfrak{F}^{(i,j)}$ is the universal form factor of the total monodromy matrix element $T_{ij}(z)$
\end{thm}

\begin{thm} \label{FF-diag}
Let $\mathbb{B}_{a,b}(\bar u;\bar v)$ be an on-shell Bethe vector and $\mathbb{C}_{a,b}(\bar u;\bar v)$ be its dual on-shell Bethe vector. Let $\mathbb{C}_{a,b}(\bar u(\bar\kappa);\bar v(\bar\kappa))$ be a twisted on-shell deformation of $\mathbb{C}_{a,b}(\bar u;\bar v)$ such that $\bar u(\bar\kappa)=\bar u$, $\bar v(\bar\kappa)=\bar v$ at\footnote{Here and below
$\bar\kappa=\{\kappa_1,\kappa_2,\kappa_3\}$, and $\bar\kappa=1$ means $\kappa_1=\kappa_2=\kappa_3=1$.} $\bar\kappa=1$. Then
\begin{equation}\label{FFdiag-eq}
\mathcal{M}^{(i,i)}
\left( m\Bigr|\begin{smallmatrix}
  \bar u & \bar u \\ \bar v & \bar v
\end{smallmatrix}\right)^{a,a}_{b,b}
=  \left(  \lambda^{(1)}_i[0]
 + (-1)^{[i]}  \frac{\mathrm{d}}{\mathrm{d} \kappa_i} \log \frac{\ell_1(\bar u(\bar\kappa))}{\ell_3(\bar v(\bar\kappa))} \Big|_{\bar \kappa=1}   \right) ||\mathbb{B}_{a,b}(\bar u;\bar v)||^2,
\end{equation}
where $\lambda^{(1)}_k[0]$ can be found from the expansion
\be{zero-modesl}
\lambda^{(1)}_k(u)=1+ \lambda^{(1)}_k[0]\,\frac cu +\dots\;,\qquad u\to\infty.
\end{equation}
\end{thm}

We prove theorems~\ref{FF-gener} and~\ref{FF-diag} in section~\ref{S-FFOD}. The most technical part of the proof is
given in appendix~\ref{A-CALC}.  Here we would like to mention only that
the general strategy of the proof is the same as in the $\mathfrak{gl}(3)$ case \cite{PakRS15b}.

Comparing these formulas with the corresponding expressions in the $\mathfrak{gl}(3)$-based models \cite{PakRS15b} we see that the only
difference is the sign factor $(-1)^{[i]}$ in \eqref{FFdiag-eq}. Certainly, the determinant formulas for the universal form factors in the
models with $\mathfrak{gl}(2|1)$ and $\mathfrak{gl}(3)$ symmetries are different, however, the relation between $\mathcal{M}^{(i,j)}$
and $\mathfrak{F}^{(i,j)}$ is the same (modulus the sign factor mentioned above). Most probably, the same relation takes place in the
general $\mathfrak{gl}(m|n)$ case as well. One should  remember, however, that in models with $\mathfrak{gl}(2|1)$-invariant $R$-matrix there exist compact determinant representations for the universal form factors \cite{HutLPRS16a}. These representations can be directly used for analysis of correlation functions.
At the same time, analogous determinant formulas for the universal form factors
in the general $\mathfrak{gl}(m|n)$ case are unknown up to date.

Observe that the dependence on the local site $m$ in \eqref{FF-gener}, \eqref{FFdiag-eq} enters only the functions $\ell_k(u)$ and $\lambda^{(1)}_k[0]$. This dependence follows
from representation \eqref{mat-T12}
\begin{equation}
\lambda^{(1)}_k(u) = \prod_{n=1}^m \lambda_k(u|n), \qquad k=1,2,3,
\end{equation}
\be{lkm}
\ell_k(u)=\prod_{n=1}^m \frac{\lambda_k(u|n)}{\lambda_2(u|n)},\qquad k=1,3,
\end{equation}
where  $\lambda_k(u|n)$ are vacuum eigenvalues of the local $L$-operators entries $\bigl(L_n(u)\bigr)_{kk}$.
Due to  \eqref{L-dep-u} the expansion of the vacuum eigenvalues $\lambda_k(u|n)$ takes the form
\be{exp-ve}
\lambda_k(u|n) = 1+  \lambda_k[0|n] \frac{c}{u} + \dots ,\qquad u\to\infty,
\end{equation}
and hence, the coefficient  $\lambda^{(1)}_k[0]$ in \eqref{zero-modesl} is
\begin{equation}
\lambda^{(1)}_k[0] = \sum_{n=1}^m \lambda_k[0|n] .
\end{equation}
Then, due to \eqref{FF-partL} we find for $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\} \neq \{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$
\be{FF-partL-answ}
\mathbb{C}_{a',b'}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})\bigl(L_m[0]\bigr)_{ij}\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})=\bigl(\mathcal{L}_m-1\bigr)
\Bigl(\prod_{n=1}^{m-1}\mathcal{L}_n\Bigr)\;
\mathfrak{F}^{(i,j)}
\left( \begin{smallmatrix}
  \bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a',a}_{b',b},
\end{equation}
where
\be{cLn}
\mathcal{L}_n=\frac{\ell_1(\bar u^{\scriptscriptstyle C}|n)\ell_3(\bar v^{\scriptscriptstyle B}|n)}{\ell_1(\bar u^{\scriptscriptstyle B}|n)\ell_3(\bar v^{\scriptscriptstyle C}|n)}.
\end{equation}
If $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\} =\{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}=\{\bar u,\bar v\}$, then
\begin{equation}\label{FFdiag-answ}
\mathbb{C}_{a,b}(\bar u;\bar v) \bigl(L_m[0]\bigr)_{ii} \mathbb{B}_{a,b}(\bar u;\bar v)
= \Bigl[\lambda_i[0|m]
+  (-1)^{[i]}\frac{\mathrm{d}}{\mathrm{d} \kappa_i} \log \frac{\ell_1\bigl(\bar u(\bar\kappa)|m\bigr)}
{\ell_3\bigl(\bar v(\bar\kappa)|m\bigr)} \Big|_{\bar \kappa=1} \Bigr]  ||\mathbb{B}_{a,b}(\bar u;\bar v)||^2,
\end{equation}
where $\bar u(\bar\kappa)$ and $\bar v(\bar\kappa)$ are the deformations of $\bar u$ and $\bar v$ described in theorem~\ref{FF-diag}.

\section{Generating functional for form factors of  partial zero modes\label{S-FFOD}}

All the form factors of the partial zero modes $T^{(1)}_{ij}[0]$ can be found from a special generating functional.
Consider an operator
\be{op-Q}
Q_{\bar\beta}=\sum_{i=1}^3 (-1)^{[i]}  \beta_i T^{(1)}_{ii}[0],
\end{equation}
where $\beta_i$ are some complex numbers.

Let  $\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$ be an on-shell Bethe vector.
Let also $\mathbb{C}^{(\kappa)}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})$ be a twisted dual on-shell Bethe vector with the twist  $\kappa=\mathop{\rm diag}\{\kappa_1,\kappa_2,\kappa_3\}$.
We stressed this fact by adding the superscript $(\kappa)$ to the vector $\mathbb{C}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})$. Suppose that
$\kappa_i=e^{\beta_i}$ and consider the following expectation value
\be{genfun-def}
\mathcal{M}^{(\kappa)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a}_{b}
= \mathbb{C}^{(\kappa)}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})\; e^{Q_{\bar\beta}}\;\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B}).
\end{equation}
Taking the derivative of this generating functional over $\beta_i$ at $\bar\kappa=1$ (that is, all $\beta_j=0$) we obtain
\be{genfun-der}
(-1)^{[i]}\mathcal{M}^{(i,i)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a,a}_{b,b}=\frac{\mathrm{d}}{\mathrm{d} \beta_i}\Bigl[
\mathcal{M}^{(\kappa)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a}_{b}
- \mathbb{C}^{(\kappa)}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})\Bigr]_{\bar\kappa=1}.
\end{equation}
It was shown in \cite{HutLPRS16b}  that for $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\}\bigr|_{\bar\kappa=1} \neq \{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$
\be{der-SP}
\frac{\mathrm{d}}{\mathrm{d} \beta_i}
 \mathbb{C}^{(\kappa)}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})\Bigr|_{\bar\kappa=1}=
(-1)^{[i]}\mathfrak{F}^{(i,i)}
\left( \begin{smallmatrix}
  \bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a,a}_{b,b}.
\end{equation}
Hence, we obtain in this case
\be{genfun-FF1}
\mathcal{M}^{(i,i)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a,a}_{b,b}=(-1)^{[i]}\frac{\mathrm{d}}{\mathrm{d} \beta_i}
\mathcal{M}^{(\kappa)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a}_{b}\Bigr|_{\bar\kappa=1}
-\mathfrak{F}^{(i,i)}
\left( \begin{smallmatrix}
  \bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a,a}_{b,b}.
\end{equation}
Thus, calculating the generating functional
\eqref{genfun-def}, we can find the form factors of the diagonal partial zero modes $T^{(1)}_{ii}[0]$ at least
for $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\}\bigr|_{\bar\kappa=1} \neq \{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$. The case $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\}\bigr|_{\bar\kappa=1} =\{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$ will be considered later.

The form factors of the partial zero modes $T^{(1)}_{ij}[0]$ with $i\ne j$ can be obtained via relation \eqref{FFlk-FFij}. Let us give an example.
Let $i=j=l=2$ and $k=1$ in \eqref{FFlk-FFij}. Then this formula takes the form
\begin{equation}\label{FF12-FF22}
\mathcal{M}^{(1,2)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a+1,a}_{b,b}
=\mathbb{C}_{a+1,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})\bigl(T_{22}^{(1)}[0]T_{12}[0]-T_{12}[0]T_{22}^{(1)}[0]\bigr)\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B}).
\end{equation}
Due to the results of \cite{HutLPRS16a} we have
\be{act-ZM}
\mathbb{C}_{a+1,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})T_{12}[0]=0, \qquad T_{12}[0]\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})=\lim_{w\to\infty}
\frac wc\mathbb{B}_{a+1,b}(\{w,\bar u^{\scriptscriptstyle B}\};\bar v^{\scriptscriptstyle B}),
\end{equation}
where we used the fact that both $\mathbb{C}_{a+1,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})$ and $\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$ are on-shell. Due to the Bethe equations
\eqref{BE} we conclude that if $\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$ is on-shell, then $\mathbb{B}_{a+1,b}(\{w,\bar u^{\scriptscriptstyle B}\};\bar v^{\scriptscriptstyle B})$ is also on-shell
at $w\to\infty$. This is because $r_1(w)\to 1$ at $w\to\infty$ according to  expansion \eqref{zero-modes}. Thus, we arrive at
\begin{equation}\label{FF12-FF22-1}
\mathcal{M}^{(1,2)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a+1,a}_{b,b}
=\lim_{w\to\infty}
\frac wc\mathbb{C}_{a+1,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})T_{22}^{(1)}[0]\mathbb{B}_{a+1,b}(\{w,\bar u^{\scriptscriptstyle B}\};\bar v^{\scriptscriptstyle B}),
\end{equation}
and since both vectors in the r.h.s. of \eqref{FF12-FF22-1} are on-shell, we obtain
\begin{equation}\label{FF12-FF22-2}
\mathcal{M}^{(1,2)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a+1,a}_{b,b}
=\lim_{w\to\infty}
\frac wc\mathcal{M}^{(2,2)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \{w,\bar u^{\scriptscriptstyle B}\} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a+1,a+1}_{b,b}.
\end{equation}
Thus, knowing an explicit representation for the form factor $\mathcal{M}^{(2,2)}$ we can find the form factor
$\mathcal{M}^{(1,2)}$ sending one of the Bethe parameters to infinity. Similarly all the other form factors of the off-diagonal
partial zero modes $T_{ij}^{(1)}[0]$ can be found.

\subsection{Calculation of the generating functional\label{SS-CGF}}

The calculation of the generating functional \eqref{genfun-def} is straightforward. First of all, we use explicit expressions
\eqref{BVcomp} and \eqref{dBVcomp} for the total (dual) Bethe vectors in terms of the partial ones. This allows us to find the action
of the operator $\exp\{Q_{\bar\beta}\}$ onto the Bethe vector $\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$. After this we obtain a sum over partitions
of the Bethe parameters involving two scalar products of the partial Bethe vectors. Using an explicit representation for the scalar
product of generic Bethe vectors in the $\mathfrak{gl}(2|1)$-based models \cite{HutLPRS16c} we find an explicit expression for the
generating functional in terms of a sum over partitions. This sum can be further simplified leading eventually to the final result.
We describe here the first steps of this derivation. The most technical part is shifted to appendix~\ref{A-CALC}.

We start with equation \eqref{BVcomp} for the total Bethe vector.
If this vector is on-shell, then we can present the product of functions $r_{1}^{(2)}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})$ as
$r_{1}^{(2)}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})=r_1(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})\ell^{-1}_1(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})$ (see \eqref{ell}) and express $r_1(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})$ in terms of Bethe equations. Then we obtain
\be{BV-BV-1}
\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})=\sum \frac{\ell_{3}(\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I})}{\ell_{1}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})} f(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I})
g(\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I},\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I})f(\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})\;
\mathbb{B}_{a_{\scriptscriptstyle \rm I\hspace{-1pt}I},b_{\scriptscriptstyle \rm I\hspace{-1pt}I}}^{(2)}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I})\mathbb{B}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}^{(1)}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I}) .
\end{equation}
Similarly, if $\mathbb{C}^{(\kappa)}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})$ is a twisted on-shell Bethe vector, then we can present the product of functions $r_{3}^{(2)}(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})$ as
$r_{3}^{(2)}(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})=r_3(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})\ell^{-1}_3(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})$  and express $r_3(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})$ in terms of the twisted Bethe equations. Then we have
\be{CV-CV-eig}
\mathbb{C}^{(\kappa)}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})=\sum \left(\frac{\kappa_2}{\kappa_3}\right)^{b_{\scriptscriptstyle \rm I}}\frac{\ell_{1}(\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I})} {\ell_{3}(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})}f(\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I}) g(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I},\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})f(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})\;
\mathbb{C}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}^{(1)}(\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I}) \mathbb{C}_{a_{\scriptscriptstyle \rm I\hspace{-1pt}I},b_{\scriptscriptstyle \rm I\hspace{-1pt}I}}^{(2)}(\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I};\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I}).
\end{equation}

Now we should compute the action of $\exp\{Q_{\bar\beta}\}$ onto the Bethe vector $\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$. Obviously,  the partial zero modes
$T^{(1)}_{ii}[0]$ act only on the partial Bethe vectors $\mathbb{B}^{(1)}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}$ and do not act on $\mathbb{B}^{(2)}_{a_{\scriptscriptstyle \rm I\hspace{-1pt}I},b_{\scriptscriptstyle \rm I\hspace{-1pt}I}}$.
The explicit action formulas are \cite{HutLPRS16a}:
\be{actBV}
\begin{aligned}
T^{(1)}_{11}[0]\mathbb{B}^{(1)}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})&=(\lambda^{(1)}_1[0]-a_{\scriptscriptstyle \rm I})\mathbb{B}^{(1)}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I}), \\
T^{(1)}_{22}[0]\mathbb{B}^{(1)}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})&=(\lambda_2^{(1)}[0]+a_{\scriptscriptstyle \rm I}-b_{\scriptscriptstyle \rm I})\mathbb{B}^{(1)}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I}), \\
T^{(1)}_{33}[0]\mathbb{B}^{(1)}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})&=(\lambda^{(1)}_3[0] -b_{\scriptscriptstyle \rm I})\;\mathbb{B}^{(1)}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I}).
\end{aligned}
\end{equation}
Here $\mathbb{B}^{(1)}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})$ is a generic partial Bethe vector. Thus, we find
\begin{multline}\label{act-BV-BV}
e^{Q_{\bar\beta}}\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})=\sum  e^{\mathcal{Q}_{\bar\beta}+ a_{\scriptscriptstyle \rm I}(\beta_2-\beta_1) +b_{\scriptscriptstyle \rm I}(\beta_3-\beta_2)}
\frac{\ell_{3}(\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I})}{\ell_{1}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})} \\
\times f(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I}) g(\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I},\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I}) f(\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})
\mathbb{B}_{a_{\scriptscriptstyle \rm I\hspace{-1pt}I},b_{\scriptscriptstyle \rm I\hspace{-1pt}I}}^{(2)}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I})\mathbb{B}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}^{(1)}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I}) ,
\end{multline}
where
\be{qb}
\mathcal{Q}_{\bar\beta}=\sum_{i=1}^3(-1)^{[i]}\beta_i\lambda_i^{(1)}[0]
\end{equation}
is the eigenvalue of the operator $Q_{\bar\beta}$ on the vector $|0\rangle$.
Substituting \eqref{CV-CV-eig} and \eqref{act-BV-BV} into \eqref{genfun-def} we arrive at
\begin{multline}\label{Mat-elS}
\mathcal{M}^{(\kappa)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a}_{b}=\sum
e^{\mathcal{Q}_{\bar\beta}+ a_{\scriptscriptstyle \rm I}(\beta_2-\beta_1)}
\frac{\ell_1(\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I})\ell_3(\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I})}{\ell_1(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})\ell_3(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})}\;
f(\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I})f(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I})g(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I},\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I}) g(\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I},\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I})\\
\times
f(\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I})f(\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I},\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})
\; \mathbb{C}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}^{(1)}(\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I})\mathbb{B}_{a_{\scriptscriptstyle \rm I},b_{\scriptscriptstyle \rm I}}^{(1)}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I}) \cdot \mathbb{C}_{a_{\scriptscriptstyle \rm I\hspace{-1pt}I},b_{\scriptscriptstyle \rm I\hspace{-1pt}I}}^{(2)}(\bar u^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I};\bar v^{\scriptscriptstyle C}_{\scriptscriptstyle \rm I\hspace{-1pt}I})\mathbb{B}_{a_{\scriptscriptstyle \rm I\hspace{-1pt}I},b_{\scriptscriptstyle \rm I\hspace{-1pt}I}}^{(2)}(\bar u^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I};\bar v^{\scriptscriptstyle B}_{\scriptscriptstyle \rm I\hspace{-1pt}I}).
\end{multline}
Thus, the problem of calculating the generating functional is reduced to the calculation of the scalar products and further summation over partitions.
Further derivation is quite technical, therefore, we give the details in appendix~\ref{A-CALC}. We would like to mention only that this derivation
goes along the same lines as in the $\mathfrak{gl}(3)$-case \cite{PakRS15b} with minor modifications. Here we formulate the final result only.

\begin{prop} \label{thmM}
Let $\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})$ be an on-shell Bethe vector and $\mathbb{C}^{(\kappa)}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})$ be a twisted dual on-shell Bethe vector.
Then
\be{result}
\mathcal{M}^{(\kappa)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a}_{b}= e^{\mathcal{Q}_{\bar\beta}}
\frac{\ell_1(\bar u^{\scriptscriptstyle C})\ell_3(\bar v^{\scriptscriptstyle B})}{\ell_1(\bar u^{\scriptscriptstyle B})\ell_3(\bar v^{\scriptscriptstyle C})}\;\mathbb{C}^{(\kappa)}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B}).
\end{equation}
\end{prop}

Suppose that $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\}\bigr|_{\bar\kappa=1} \neq \{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$. Then, it follows immediately from \eqref{result} that
\be{der-result}
\frac{\mathrm{d}}{\mathrm{d} \beta_i}\mathcal{M}^{(\kappa)}\left(m\Bigr|\begin{smallmatrix}
\bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\
\bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a}_{b}\Bigr|_{\bar\kappa=1}=
\frac{\ell_1(\bar u^{\scriptscriptstyle C})\ell_3(\bar v^{\scriptscriptstyle B})}{\ell_1(\bar u^{\scriptscriptstyle B})\ell_3(\bar v^{\scriptscriptstyle C})}\;
\frac{\mathrm{d}}{\mathrm{d} \beta_i}\mathbb{C}^{(\kappa)}_{a,b}(\bar u^{\scriptscriptstyle C};\bar v^{\scriptscriptstyle C})\mathbb{B}_{a,b}(\bar u^{\scriptscriptstyle B};\bar v^{\scriptscriptstyle B})\Bigr|_{\bar\kappa=1}.
\end{equation}
Indeed, for $\bar\kappa=1$ the scalar product in \eqref{result} turns into the scalar product of two different on-shell Bethe vectors. Hence,
it vanishes. Therefore, the $\beta_i$-derivative must act on this scalar product only, otherwise we obtain zero contribution. Using then
\eqref{der-SP} and \eqref{genfun-FF1} we immediately obtain
\begin{equation}\label{FFii-gener}
  \mathcal{M}^{(i,i)}
\left( m\Bigr|\begin{smallmatrix}
  \bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a,a}_{b,b}

  = \left( \frac{\ell_1(\bar u^{\scriptscriptstyle C})\ell_3(\bar v^{\scriptscriptstyle B})}{\ell_1(\bar u^{\scriptscriptstyle B})\ell_3(\bar v^{\scriptscriptstyle C})} -1 \right) \mathfrak{F}^{(i,i)}
\left( \begin{smallmatrix}
  \bar u^{\scriptscriptstyle C} & \bar u^{\scriptscriptstyle B} \\ \bar v^{\scriptscriptstyle C} & \bar v^{\scriptscriptstyle B}
\end{smallmatrix}\right)^{a,a}_{b,b},
\end{equation}
for $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\} \neq \{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}$. Thus, we reproduce \eqref{FF-gener} for $i=j$. The form factors of the off-diagonal partial zero modes
then can be derived via \eqref{FFlk-FFij}.

Consider finally the case $\{\bar u^{\scriptscriptstyle C},\bar v^{\scriptscriptstyle C}\}\bigr|_{\bar\kappa=1} = \{\bar u^{\scriptscriptstyle B},\bar v^{\scriptscriptstyle B}\}=\{\bar u,\bar v\}$. Then \eqref{result} yields
\begin{multline}\label{genfun-FF10}
\frac{\mathrm{d}}{\mathrm{d} \beta_i}
\mathcal{M}^{(\kappa)}\left(m\Bigr|\begin{smallmatrix}
\bar u & \bar u \\
\bar v & \bar v
\end{smallmatrix}\right)^{a}_{b}\Bigr|_{\bar\kappa=1}= \frac{\mathrm{d}}{\mathrm{d} \beta_i}\mathbb{C}^{(\kappa)}_{a,b}(\bar u(\bar\kappa);\bar v(\bar\kappa))\mathbb{B}_{a,b}(\bar u;\bar v)\Bigr|_{\bar\kappa=1}
\\
+
\left( (-1)^{[i]}\lambda_i^{(1)}[0]  +  \frac{\mathrm{d}}{\mathrm{d} \beta_i} \log \frac{\ell_1(\bar u(\bar\kappa))}{\ell_3(\bar v(\bar\kappa))} \Big|_{\bar \kappa=1}   \right) ||\mathbb{B}_{a,b}(\bar u;\bar v)||^2.
\end{multline}
Comparing this equation with \eqref{genfun-der} we arrive at
\begin{equation}\label{FFdiag-eq0}
\mathcal{M}^{(i,i)}
\left( m\Bigr|\begin{smallmatrix}
  \bar u & \bar u \\ \bar v & \bar v
\end{smallmatrix}\right)^{a,a}_{b,b}
=  \left( \lambda_i^{(1)}[0]
  + (-1)^{[i]}  \frac{\mathrm{d}}{\mathrm{d} \beta_i} \log \frac{\ell_1(\bar u(\bar\kappa))}{\ell_3(\bar v(\bar\kappa))} \Big|_{\bar \kappa=1}   \right) ||\mathbb{B}_{a,b}(\bar u;\bar v)||^2.
\end{equation}
Finally, using
\be{der-der}
\frac{\mathrm{d}}{\mathrm{d} \beta_i}\Big|_{\bar \kappa=1}=\frac{\mathrm{d}}{\mathrm{d} \log\kappa_i}\Big|_{\bar \kappa=1}
=\frac{\mathrm{d}}{\mathrm{d} \kappa_i}\Big|_{\bar \kappa=1},
\end{equation}
we reproduce \eqref{FFdiag-eq}.

\section*{Conclusion}

In this paper we have calculated FFLO in the integrable models with $\mathfrak{gl}(2|1)$-invariant
$R$-matrix. We have shown that these form factors are proportional to the universal form factors. Our derivation is not based
on a specific representation of the $RTT$-algebra, and thus, it is valid for a wide class of integrable models solvable by the
algebraic Bethe ansatz (besides the t-J model see e.g. \cite{BraGLZ95,BedF95,RamM96,PfaF96}). Due to an isomorphism between
$Y\bigl(\mathfrak{gl}(2|1)\bigr)$ and $Y\bigl(\mathfrak{gl}(1|2)\bigr)$ these results also can be applied to the
$\mathfrak{gl}(1|2)$-invariant models.

It is worth mentioning that for calculating the FFLO it was enough to find a special generating functional
\eqref{genfun-def}. All the form factors then can be found by taking $\kappa$-derivatives of this generating functional and sending
some of the Bethe parameters to infinity. These close relations between different form factors occur due to the commutation relations between
zero modes $T_{ij}[0]$ in the expansion \eqref{zero-modes}.
Generically, local $L$-operators not necessarily have the form \eqref{L-dep-u}, therefore, the monodromy matrix
might have an asymptotic expansion different from \eqref{zero-modes}.
This may lead to a redefinition of the zero modes and changing of the commutation relations between them, which in turn
yields a modification of the relations between the FFLO (see e.g. \cite{PakRS15c}). Nevertheless,
the latter still are proportional to the universal form factors. Therefore, the results of this paper for the most part apply
to the models in which the monodromy matrix does not have the asymptotic expansion \eqref{zero-modes}.

Our results are in complete analogue with the ones obtained in \cite{PakRS15b} for the models with the
$\mathfrak{gl}(3)$-invariant $R$-matrix. One can conjecture that the relationship between the FFLO
and the universal form factors remains true in the general $\mathfrak{gl}(m|n)$ case. However, compact representations for the universal
form factors in the $\mathfrak{gl}(m|n)$-invariant models are not known for today. At the same time, in the models with
$\mathfrak{gl}(2|1)$ or $\mathfrak{gl}(1|2)$ symmetries the universal form factors were calculated in \cite{HutLPRS16a} in terms
of determinants. These representations allow one to use our results for studying correlation functions via the form factor expansion.

\section*{Acknowledgements}

The work of J.F. has been supported by the Grant Agency of the Czech Technical
University in Prague, grant No. $SGS15/215/OHK4/3T/14$, and by the Grant of
the Plenipotentiary of the Czech Republic at JINR, Dubna. He is also grateful to A. P. Isaev and \v{C}. Burd\'{i}k for discussions and
support.

A part of this work, section~\ref{S-FFOD}, was performed in Steklov Mathematical Institute
of Russian Academy of Sciences by N.A.S. and he was supported by the Russian Science Foundation
under grant 14-50-00005.