OK, here goes.

We start with changing the notation ($z\to 20z^2$, $-z-3\to r$, $20rz\to y$ means that what was denoted by $z$ will be denoted by $20z^2$ from now on, $r$ is $-z-3$ with *new* $z$, so it is $-\sqrt{z/20}-3$ in terms of old $z$, and $y$ denotes $20rz$ with just (re)defined $r,z$).

So, execute the following sequence (the order matters!)
$$
x\to e^{-x},\ a\to e^{-a},\ mkja\to T,\ ja\to a
$$
Denote $\varphi(t)=-\log(1-e^{-t})$. For a function $f$ denote by $LRS(f,u,v; h)$ and $RRS(f,u,v;h)$ the left and the right Riemann sums of $f$ on the interval $[u,v]$ with step $h$ respectively (we will always assume that $(v-u)/h\in \mathbb N\cup\{\infty\}$).

Then the above substitutions induce the substitution
$$
-\log f_{mk}(a^jx)\to \frac{LRS(\varphi,x,\infty;a)+xT-RRS(\varphi,0,T;a)}a
$$

The claim to prove reduces to the statement that this function decreases in $a$ as long as we run $a$ along admissible steps for $[0,T]$ with fixed $T$.

We have already seen that replacing Riemann sums by integrals results in the area $A\ge 0$ of the "excess triangle" formed by the vertical line through $x$, the horizontal line through $T$, and the graph of $\varphi$ in the coordinate system where $x$ runs over the horizontal axis and $T$ over the vertical one in the numerator. Thus, the function to investigate is
$$
\frac Aa+\frac {\Phi(a)}a+\frac{\Psi(a)}a
$$
where $\Phi(a)=LRS(\varphi,x,\infty;a)-\int_x^\infty \varphi(t)\,dt$ and $\Psi(a)=\int_0^T\varphi(t)\,dt-RRS(\varphi,0,T,a)$

Note that we have a nice series expansion
$$
\varphi(t)=\sum_{k\ge 1}\frac 1ke^{-kt}
$$
The integration and the Riemann sum computation are linear operations with respect to the function, so we have
$$
\Phi(a)=\sum_{k\ge 1}\frac 1k\Phi_k(a),\qquad \Psi(a)=\sum_{k\ge 1}\frac 1k\Psi_k(a)
$$
where $\Phi_k$ and $\Psi_k$ are defined in the same way but using $e^{-kt}$ instead of $\psi(t)$.

Now, we (or, if you prefer, at least I) cannot integrate $\varphi$ or to sum it along an arithmetic progression. However everybody can do it for an exponential function. The direct computation yields
$$
\Phi_k(a)=e^{-kx}\left[\frac a{1-e^{-ka}}-\frac 1k\right],\qquad \Psi_k(a)=(1-e^{-kT})\left[\frac 1k-\frac a{e^{ka}-1}\right]
$$
Note that in the computation of $\Psi_k(a)$ we have used the fact that $T/a$ is an integer. However the resulting answer is formally defined for all $a>0$. So we will not use the "fitting condition" anywhere below.

Notice also that
$$
\left[\frac a{1-e^{-ka}}-\frac 1k\right]+\left[\frac 1k-\frac a{e^{ka}-1}\right]=a
$$
so
$$
\frac{\Phi_k(a)+\Psi_k(a)}a=\frac{1-e^{-kx}-e^{-kT}}k\left[\frac 1a-\frac k{e^{ka}-1}\right]+ \operatorname{const}
$$
Thus we need to show that
$$
\frac Aa+\sum_{k\ge 1}\frac{1-e^{-kx}-e^{-kT}}{k^2}\left[\frac 1a-\frac k{e^{ka}-1}\right]
$$
is decreasing or, equivalently, that
$$
\frac A{a^2}+\sum_{k\ge 1}\frac{1-e^{-kx}-e^{-kT}}{k^2}\left[\frac 1{a^2}-\frac {k^2}{e^{ka}+e^{-ka}-2}\right]\ge 0\,.
$$
When the excess triangle is absent or is above the graph, we have $e^{-x}+e^{-T}\le 1$ and the factor $k\ge 1$ can only improve this inequality, so each term is non-negative and the estimate is trivial, thus justifying my first remark.

The other case is slightly more interesting.

Let $t>0$ satisfy $e^{-t}+e^{-x}=1$. We have now $t>T$. Write
$$
A=\int_T^{t}\varphi(s)\,ds-(t-T)\varphi(t)=\sum_{k\ge 1}
\frac 1k\left[\int_T^{t}e^{-ks}\,ds-(t-T)e^{-kt}
\right]
\\
=\sum_{k\ge 1}\frac{e^{-kT}-e^{-kt}-k(t-T)e^{-kt}}{k^2}
$$
We want to combine these terms multiplied by $\frac 1{a^2}$ with the corresponding terms in the main sum. It will be more convenient to diminish them first by multiplying them not by the full $\frac 1{a^2}$ but just by the expressions in the brackets in the main sum. Then the result can be rewritten as
$$
\sum_{k\ge 1}\frac{1-e^{-kx}-e^{-kt}-kte^{-kt}}{k^2}\left[\frac 1{a^2}-\frac {k^2}{e^{ka}+e^{-ka}-2}\right]+
\sum_{k\ge 1}\frac{Te^{-kt}}{k}\left[\frac 1{a^2}-\frac {k^2}{e^{ka}+e^{-ka}-2}\right]
$$
The second sum is clearly non-negative. Let's look at the first one.
We have
$$
\int_u^\infty\varphi(s)\,ds=\sum_{k\ge 1}\frac 1k\int_u^\infty e^{-ks}\,ds=\sum_{k\ge 1}\frac{e^{-ku}}{k^2}\,.
$$
Using this with $u=0,x,t$ and flipping the geometric picture for the integral of $\varphi$ over $[t,\infty)$ to the vertical axis, as usual, we see that
$$
\sum_{k\ge 1}\frac{1-e^{-kx}-e^{-kt}}{k^2}=tx
$$
(the area of the remaining rectangle).
On the other hand,
$$
\sum_{k\ge 1}\frac{kte^{-kt}}{k^2}=t\varphi(t)=tx
$$
as well. So the first factors form a sequence whose numerators are (obviously, because $u\mapsto(1+u)e^{-u}$ is decreasing on $[0,+\infty)$) increasing and whose sum equals $0$. Thus, they break from $-$ to $+$ just once when $k$ goes up. However the second factors are increasing in $k$. Thus, the sum we are interested in is non-negative as well.

The End