David Carlisle has given the right hints in his comments but he has
not written an answer. So I try to formulate a complete answer: Change
the value of \fontdimen16 for \textfont2. It should get the value that
\fontdimen17 has.
Here are some executable scenarios in plain TeX; they can be executed
standalone:
% the current situation
\centerline{$
f_{k,j} \quad f_{k,j}^{\vphantom{k}} \qquad a_{k,j} \quad a_{k,j}^{\vphantom{k}}
$}
$$
f_{k,j} \quad f_{k,j}^{\vphantom{k}} \qquad a_{k,j} \quad a_{k,j}^{\vphantom{k}}
$$
% change fontdimen 16; note: this is always a global change
\fontdimen16\textfont2=\fontdimen17\textfont2
\centerline{$
f_{k,j} \quad f_{k,j}^{\vphantom{k}} \qquad a_{k,j} \quad a_{k,j}^{\vphantom{k}}
$}
$$
f_{k,j} \quad f_{k,j}^{\vphantom{k}} \qquad a_{k,j} \quad a_{k,j}^{\vphantom{k}}
$$
% increase the fontdimens 16 and 17
\fontdimen17\textfont2=1.25\fontdimen17\textfont2 % 1.25 is a guess
\fontdimen16\textfont2=\fontdimen17\textfont2
\centerline{$
f_{k,j} \quad f_{k,j}^{\vphantom{k}} \qquad a_{k,j} \quad a_{k,j}^{\vphantom{k}}
$}
$$
f_{k,j} \quad f_{k,j}^{\vphantom{k}} \qquad a_{k,j} \quad a_{k,j}^{\vphantom{k}}
$$
\bye
I don't think it is a 'hack' to change the value of the \fontdimen16.
The TeXbook, p. 179, contains the following text, which states more
or less your problem:
Instead of changing the sizes of subformulas, or using |\raise|, you
can also control vertical spacing by changing the parameters that
\TeX\ uses when it is converting math lists to horizontal lists. These
parameters are described in Appendix G; you need to be careful
when changing them, because such changes are global (i.e., not local
to groups). Here is an example of how such a change might be made:
Suppose that you are designing a format for chemical typesetting, and
that you expect to be setting a lot of formulas like
$\rm Fe_2^{+2}Cr_2O_4$.
You may not like the fact that the subscript
in $\rm Fe_2^{+2}$ is lower than the subscript in $\rm Cr_2$; and you
don't want to force users to type monstrosities like
\begintt
$\rm Fe_2^{+2}Cr_2^{\vphantom{+2}}O_4^{\vphantom{+2}}$
\endtt
just to get
the formula $\rm Fe_2^{+2}Cr_2^{\vphantom{+2}}O_4^{\vphantom{+2}}$
with all subscripts at the same level. Well, all you need to do is set
'\fontdimen16\tensy=2.7pt|' and '|\fontdimen17\tensy=2.7pt|', assuming
that |\tensy| is your main symbol font (|\textfont2|); this lowers all
normal subscripts to a position $2.7\pt$ below the baseline, which is
enough to make room for a possible superscript that contains a plus
sign.
The details about the fontdimens are described in Appendix G of The
TeXbook. A very nice TUGboat article by B. Jackowski explains the
rules of this Appendix with diagrams; see figures 9--11 for the
placement of subscripts.