OK, I'll post my argument then. It will be a somewhat long story so I'll do it in a few stretches (I don't think I'll have a single continuous time interval to make the full post, so I apologize for bumping it to the front page a few times when editing to add the stuff).
First of all, let us compute the second derivative of the logarithm of $\sigma_n(x)$. It is
$$
\frac{\sigma_n''(x)}{\sigma_n(x)}-\left[\frac{\sigma_n'(x)}{\sigma_n(x)}\right]^2=\frac{\sum_{k=1}^n k^x\log^2 k}{\sum_{k=1}^n k^x}-
\left[\frac{\sum_{k=1}^n k^x\log k}{\sum_{k=1}^n k^x}\right]^2
\\
=\frac 12\frac{\sum_{k,m=1}^n k^xm^x(\log k-\log m)^2}{\sum_{k,m=1}^n k^xm^x}
\\
=\frac 12\frac{\sum_{k,m:1\le k<m\le n} k^xm^x \log^2 \frac km}{\sum_{k,m:1\le k<m\le n} k^xm^x+\frac 12\sum_{m:1\le m\le n} m^xm^x}
\\
=\frac 12\frac{\sum_{m=1}^n m^x\left(\sum_{k:1\le k<m} k^x \log^2 \frac km\right)}{\sum_{m=1}^nm^x\left(\sum_{k:1\le k<m} k^x+\frac 12 m^x\right)}
=
\frac 12\frac{\sum_{m=1}^n w_mE_m}{\sum_{m=1}^n w_m}\,,
$$
where
$$
w_m=m^x\left(\sum_{k:1\le k<m} k^x+\frac 12 m^x\right),\
E_m=\frac{\sum_{k:1\le k<m} k^x \log^2 \frac km}{\sum_{k:1\le k<m} k^x+\frac 12 m^x}.
$$
Thus, to show that the second derivative of $\log\sigma_{n+1}(x)$ is at least as large as that of $\log\sigma_{n}(x)$, it is sufficient to show that $E_n$ is a non-decreasing sequence in $n$.
Dividing the numerator and the denominator in the definition of $E_n$ by $n^{x+1}$, we see that for $x>0$, $E_n$ is the ratio of the $n$-th trapezoid sums $T_ng$ and $T_nf$ of the functions $g(z)=z^x\log^2z$ and $f(z)=z^x$ on $[0,1]$. It will be convenient from now on to use $p$ instead of $x$, so that $x$ could be used for something else as needed.
The trapezoid sum of a reasonably smooth on $[0,1]$ function $h$ can be also expressed as
$$
T_n h=\sum_{m:n\mid m}\widehat h(m)\,.
$$
If, in addition, $h$ is real-valued, then we can take the real part of both sides and get
$$
T_n h=\widehat h(0)+2\sum_{m:m>0, n\mid m}\Re\widehat h(m)\,,
$$
so, if the real parts of the Fourier coefficients $\widehat h(n)$ with positive indices have certain monotonicity for $n\ge n_0$, the trapezoid sums $T_nh$ will be monotone in the same way for $n\ge n_0$ (note that the direction of the monotonicity of the real parts of the Fourier coefficients of $h$ determines their signs as well since they tend to $0$ at infinity).
It has been shown in this thread that for $p>1$, the real parts of the Fourier coefficients of $f(z)$ are decreasing all the way from $n=1$ to $n=\infty$, so for $p>1$ the denominators $T_nf$ of $E_n$ form a positive decreasing sequence (for $p\ge 2$ one can also give an elementary real variable proof of this fact). Thus, everything will be fine if we show that the numerators $T_ng$ form an increasing sequence.
That is not always true (and, technically, we do not need that much to show that $E_n$ increase, which seems to be a correct statement for all $p>0$) but the numeric experiments indicate that it is, indeed, so for $p\ge 4$ and we will (try to) prove that it is true for $p\ge p_0$ with some large $p_0>0$.
We shall start with the trivial reason that forces $T_ng$ to increase. Suppose that $g$ is convex on $[0,\frac n{n+1}]$. Then, using the representation
$$
\frac kn=\frac{n-k}{n}\frac k{n+1}+
\frac{k}{n}\frac {k+1}{n+1}
$$
for $k=1,\dots,n-1$, we get
$$
g(\tfrac kn)\le \tfrac{n-k}{n}g(\tfrac k{n+1})+
\tfrac{k}{n}g(\tfrac {k+1}{n+1})\,,
$$
so
$$
T_ng=\frac 1n\sum_{k=1}^{n-1}g(\tfrac kn)\le
\frac 1n\sum_{k=1}^{n-1}[\tfrac{n-k}{n}g(\tfrac k{n+1})+
\tfrac{k}{n}g(\tfrac {k+1}{n+1})]
\\
=\frac{n-1}{n^2}\sum_{k=1}^n g(\tfrac k{n+1})=\frac{n^2-1}{n^2}T_{n+1}g\le T_{n+1}g\,.
$$
Let's now find the range of $z\in[0,1]$ for which $g''(z)>0$. We have $g'(z)=z^{p-1}(p\log^2 z+2\log z)$ and $g''(z)=z^{p-2}[p(p-1)\log^2 z+(4p-2)\log z+2]$. The expression in parentheses is positive at $0$ if $p>1$ and hits $0$ when
$$
\log z=-\frac{(2p-1)+\sqrt{2p^2-2p+1}}{p(p-1)}\approx -\frac{2+\sqrt 2}{p}
$$
for large $p>1$. Thus, for every fixed $\delta>0$, the convexity assumption holds as long as $\frac pn>2+\sqrt 2+\delta$ and $p>p(\delta)$.
We shall show that we have $\Re\widehat g(n+1)\ge \Re\widehat g(n)$ when $\frac pn<4-\delta$ and $p>p(\delta)$.
To this end we shall flip the function and write
$$
\Re\widehat g(n+1)-\Re\widehat g(n)=\int_0^1 e^{2\pi i nz}[e^{2\pi iz}-1](1-z)^p\log^2(1-z)\,dz=\int_0^1 G_{p,n}(z)\,dz\,.
$$
We shall do the contour integration again (like in that thread I mentioned already with an escape to $+i\infty$) but the details will be uglier.
First, let us take care of the return, which will again be done over the right side $1+iy$, $y\ge 0$ of the half-strip. We note that
$$
|G_{p,n}(1+iy)|=\left[\log^2 y+\left(\tfrac \pi 2\right)^2\right](1-e^{-2\pi y})y^p e^{-2\pi ny}\,.
$$
Recall that we are interested in $n>p/4$ and large $p$. Let us consider $|G_{1,1/4}(1+iy)|$. Choose $Y'>0$ so that $y\mapsto ye^{-\frac \pi 2y}$ is decreasing on $[Y',+\infty)$. Now choose $Y>0$ so that
$$
\int_{Y'}^Y|G_{1,1/4}(1+iy)|\,dy\ge \int_{Y}^{+\infty}|G_{1,1/4}(1+iy)|\,dy
$$
Then this inequality will be preserved if we multiply the integrand by any decreasing on $[Y',+\infty)$ function of $y$. In particular, we can multiply it by
$[ye^{-\frac \pi 2y}]^{p-1}e^{-2\pi(n-\frac p4)y}$ to get
$$
\int_{Y'}^Y|G_{p,n}(1+iy)|\,dy\ge \int_{Y}^{+\infty}|G_{p,n}(1+iy)|\,dy
$$
for all $p>1$, $n>\frac p4$. Thus the return integral will be at most
$$
2\int_{0}^Y|G_{p,n}(1+iy)|\,dy
$$
with some fixed $Y$ independent of $p$ and $n$ in our range.
Now we want to escape from $0$ to $+i \infty$ along some curvilinear contour. Note that $\left(\frac{\log(1-z)}z\right)^2$ and $\frac{e^{2\pi iz}-1}z$ are nice analytic functions with bounded argument in some fixed size neighborhood of our half-strip minus the disk of that size centered at $1$. Thus the argument of their product is some bounded (in the same neighborhood) harmonic function $U(z)$ with $U(0)=\frac \pi 2$, which then has bounded derivatives as well. Thus, the integral over the escaping contour is
$$
\int |G_{p,n}(z)|e^{i(-p\alpha+2\pi nx+U(z) +3\arg(z)+\theta)}|dz|=\int |G_{p,n}(z)|e^{i\Phi(z)}|dz|\,
$$
in the same notation as in the thread about monotonicity of Fourier coefficients.
We shall choose below the phase $\Phi$ so that $\cos\Phi\ge c\sin\alpha$ with some absolute $c>0$. Let us show now that this will suffice.
Again, we shall compare $|G_{n,p}(z)|$ with $|G_{n,p}(z')|$ where $z'=1+i\Im z$ when $0\le \Im z\le Y$ and $z$ runs over our contour. $|e^{2\pi i nz}|=|e^{2\pi i nz'}|$, $|e^{2\pi i z-1}|\ge |e^{2\pi i z'}-1|$, so these parts are fine.
Now $|1-z'|^p=|1-z|^{p}\frac{\Im z}{|1-z|}\sin^{p-1}\alpha$, so we have a huge gain here if $\alpha$ stays separated from $\frac\pi 2$, which will be true if $x$ stays away from $1$ (recall that $\Im z$ is bounded by $Y$ in that comparison). The problem is the $log$ factor.
In $|G_{p,n}(z)|$ we can guarantee just $\alpha^2$ (the square of the imaginary part) while in $|G_{p,n}(z')|$ we have $\log^2\Im z
+\left(\frac\pi 2\right)^2$. Fortunately, this expression times $\frac{\Im z}{|1-z|}$ stays bounded if we keep away from $1$ and below $Y$, so we see that to get $c\sin\alpha |G_{p,n}(z)|>2|G_{p,n}(z')|$, we just need to ensure that
$$
c\sin\alpha \alpha^2\ge C\sin^{p-1}\alpha\,,
$$
which, indeed, holds for large enough $p$ that depend on $c$, $Y$ and on how well we stay separated from $1$ in $x=\Re z$.
It remains to construct the contour.
To be continued...
