I'm leaving my previous approach at the end, but here's my much simpler proof of the inequality.
If $u\le v\le w$, we get $0\le(v-u)(w-v)=uv+vw-uw-v^2$:
ie $uv+vw\ge uw+v^2$. We can apply this to $u,v,w=\sqrt{a_i},\sqrt{a_{i+1}},\sqrt{a_{i+2}}$ to get
$$
\sqrt{a_ia_{i+1}}+\sqrt{a_{i+1}a_{i+2}}\ge
\sqrt{a_ia_{i+2}}+a_{i+1}
\quad\text{for}\quad i=1,\ldots,n-2
$$
which makes
$$
\begin{align}
2\left(\sqrt{a_ia_{i+1}}+\sqrt{a_{i+1}a_{i+2}}\right)\ge &
\sqrt{a_ia_{i+1}}+\sqrt{a_{i+1}a_{i+2}}+\sqrt{a_ia_{i+2}}+a_{i+1}\\
\ge & 3\,\sqrt[3]{a_ia_{i+1}a_{i+2}} + a_{i+1}.
\end{align}
$$
This makes
$$
\sum_{i=1}^{n-2} \sqrt[3]{a_ia_{i+1}a_{i+2}}+\frac{a_{i+1}}3
\le \frac{2}{3}\left(\sqrt{a_1a_2}+\sqrt{a_{n-1}a_n}\right)
+\frac43\left(\sqrt{a_2a_3}+\cdots+\sqrt{a_{n-2}a_{n-1}}\right).
$$
Applying this, we get
$$
\begin{align}
\sum_\text{cyc}\sqrt{a_ia_{i+1}}
\ge&
\sqrt{a_1a_n}+\frac{\sqrt{a_1a_2}+\sqrt{a_{n-1}a_n}}{3}
-\frac13\sum_{i=2}^{n-2}\sqrt{a_ia_{i+1}}\\
&+\sum_{i=1}^{n-2}\sqrt[3]{a_ia_{i+1}a_{i+2}}
+\frac13\sum_{i=2}^{n-1}a_i.
\end{align}
$$
Using $D$ to denote the difference between the two sides,
$$
\begin{align}
D:=&\sum_\text{cyc}\sqrt{a_ia_{i+1}}-\sqrt[3]{a_ia_{i+1}a_{i+2}}\\
\ge&
\sqrt{a_1a_n}+\frac{\sqrt{a_1a_2}+\sqrt{a_{n-1}a_n}}{3}
-\frac13\sum_{i=2}^{n-2}\sqrt{a_ia_{i+1}}\\
&-\sqrt[3]{a_1a_2a_n}-\sqrt[3]{a_1a_{n-1}a_n}
+\frac13\sum_{i=2}^{n-1}a_i
\end{align}
$$
to which we apply $\sqrt{a_ia_{i+1}}\le(a_i+a_{i+1})/2$ to obtain
$$
D\ge
\sqrt{a_1a_n}+\frac{\sqrt{a_1a_2}+\sqrt{a_{n-1}a_n}}{3}
-\sqrt[3]{a_1a_2a_n}-\sqrt[3]{a_1a_{n-1}a_n}
+\frac{a_2+a_{n-1}}6.
$$
Next, using $x^sy^{1-s}\le sx+(1-s)y$ for $x,y\ge0$, $0\le s\le1$, entering $x=\sqrt{a_1a_2}, y=a_n, s=2/3$ makes $\sqrt[3]{a_1a_2a_n}\le\frac23\sqrt{a_1a_n}+\frac13\sqrt{a_2}$, etc. This makes
$$
\begin{align}
D\ge&
\frac{\sqrt{a_1a_2}+\sqrt{a_{n-1}a_n}-\sqrt{a_1a_n}}{3}
-\frac{a_2+a_{n-1}}{6}\\
=&
\frac{\left(\sqrt{a_n}-\sqrt{a_2}\right)
\left(\sqrt{a_{n-1}}-\sqrt{a_1}\right)}{3}
-\frac{\left(\sqrt{a_{n-1}}-\sqrt{a_2}\right)^2}{6}\\
\ge&
\frac{\left(\sqrt{a_{n-1}}-\sqrt{a_2}\right)^2}{6}
\ge 0.
\end{align}
$$
And $D\ge0$ is exactly what we wanted to prove.
We also see from this that equality requires $a_1=a_2=a_{n-1}=a_n$ which means the $a_i$ are all equal.
Here is my previous approach, somewhat too briefly described. I'll try to update it since the middle step was somewhat unusual. I have checked the steps by simulations (random $a_i$) and think the method should be sound, but there may be typos and mistakes even so.
Let $A(a_1,\ldots,a_n)=\sum_\text{cyc}\sqrt{a_ia_{i+1}}$ and
$B(a_1,\ldots,a_n)=\sum_\text{cyc}\sqrt[3]{a_ia_{i+1}a_{i+2}}$. We want to show that $(A-B)(a_1,\ldots,a_n)\ge0$ whenever $0\le a_1\le\cdots\le a_n$.
The proof has three steps.
Show that $(A-B)(a_1,\ldots,a_n)\le(A-B)(a_1,\ldots,a_{i-1},a_{i+1},\ldots,a_n)$ for $3\le i\le n-2$.
Starting with $a_1,\ldots,a_n$, for any $\epsilon>0$ we can fill in the gaps between $a_i$ and $a_{i+1}$ forming a new sequence $a'_1\le\cdots\le a'_N$ containing $a_i$ as a subsequence, so that $a'_1=a_1$, $a'_2=a_2$, $a'_{n-1}=a_{N-1}$, $a'_{n}=a_{N}$, and $a'_{i+1}-a'_i<\epsilon$ for $2\le i<N-1$. Then $A(a)-B(a)\ge A(a')-B(a')$ by repeated use of the above inequality starting with $a'_1,\ldots,a'_N$ and removing elements one at a time to obtain $a_1,\ldots,a_n$.
Taking the limit as $\epsilon\rightarrow0$, we get a lower bound on $A(a)-B(a)\ge A(a')-B(a')$ which depends only on $a_1,a_2,a_{n-1},a_{n}$ which can be verified.
Proofs of the inequalities in 1. and 3. use the same approach as explained in the question: express inequality as polynomial in non-negative variables, and check for negative coefficients. I did this using the sympy Python package since I was already doing numerical simulations in Python to test hypotheses, but Sage is a good alternative.
Denoting $a=(a_1,\ldots,a_n)$ and $a^{(i)}=(a_1,\ldots,a_{i-1},a_{i+1},\ldots,a_n)$,
$$
\begin{align}
A(a)-A(a^{(i)}) =& \sqrt{a_{i-1}a_i}+\sqrt{a_ia_{i+1}}-\sqrt{a_{i-1}a_{i+1}},
\\
B(a)-B(a^{(i)}) =&
\sqrt[3]{a_{i-2}a_{i-1}a_{i}}+\sqrt[3]{a_{i-1}a_{i}a_{i+1}}
+\sqrt[3]{a_{i}a_{i+1}a_{i+2}}
-\sqrt[3]{a_{i-2}a_{i-1}a_{i+1}}-\sqrt[3]{a_{i-1}a_{i+1}a_{i+2}}.
\end{align}
$$
If $3\le i\le n-2$, we can eliminate $a_{i-2}$ and $a_{i+2}$ by
$$
\begin{align}
B(a)-B(a^{(i)})
=& \sqrt[3]{a_{i-1}a_{i}a_{i+1}}
-\sqrt[3]{a_{i-2}a_{i-1}}\cdot\left(\sqrt[3]{a_{i+1}}-\sqrt[3]{a_{i}}\right)
+\sqrt[3]{a_{i+1}a_{i+2}}\cdot\left(\sqrt[3]{a_{i}}-\sqrt[3]{a_{i-1}}\right)
\\
\ge& \sqrt[3]{a_{i-1}a_{i}a_{i+1}}
-\sqrt[3]{a^2_{i-1}}\cdot\left(\sqrt[3]{a_{i+1}}-\sqrt[3]{a_{i}}\right)
+\sqrt[3]{a^2_{i+1}}\cdot\left(\sqrt[3]{a_{i}}-\sqrt[3]{a_{i-1}}\right)
\end{align}
$$
If we let $a_{i-1}=x^6$, $a_i=(x+y)^6$, $a_{i+1}=(x+y+z)^6$ with $x,y,z\ge0$, we can express
$$
\left[A(a^{(i)})-B(a^{(i)})\right]
-\left[A(a)-B(a)\right]
\ge P(x,y,z)
$$
where $P(x,y,z)$ is a degree 6 polynomial with all non-negative coefficients. So, for $x,y,z\ge0$, we get $P(x,y,z)\ge0$, and so
$$
A(a^{(i)})-B(a^{(i)}) \ge A(a)-B(a).
$$
The original idea was to use induction starting with $n=4$, but that would require an inequality in the opposite direction: I confirmed this using simulations with random sequences.
Instead, I can insert extra values into the sequence, which reduces the difference $A-B$, and prove that $A-B\ge0$ even as $n$ increases. We can do so, except from the first and last two elements in the list which are not changed, and in such a way that we get a new sequence $a'_1\le\cdots\le a'_N$ as described in 2. Since $a'_{i+1}-a'_i<\epsilon=O(1/N)$ for $3\le i\le N-2$, we get
$$
\sqrt{a'_ia'_{i+1}}-2\sqrt[3]{a'_ia'_{i+1}a'_{i+2}}+\sqrt{a'_{i+1}a'_{i+2}}
= O(\epsilon^2)
$$
as long as we take care to spread the $a'_i$ out evenly, eg $a'_i=\sqrt{a'_{i-1}a'_{i+1}}$ for all $a'_i$ not in the original sequence, and so the sum
$$
\sum_{i=3}^{n-2}
\sqrt{a'_ia'_{i+1}}-2\sqrt[3]{a'_ia'_{i+1}a'_{i+2}}+\sqrt{a'_{i+1}a'_{i+2}}
\rightarrow 0\quad\text{as }\epsilon\rightarrow 0,\; N\rightarrow\infty.
$$
Using the above limit, what remains of the sum $A(a')-B(a')$ becomes
$$
\begin{align}
A(a')-B(a')\rightarrow&
\sqrt{a_1a_2}+\sqrt{a_1a_n}+\sqrt{a_{n-1}a_n}+\frac{a_2+a_{n-1}}{2}
\\
&-\left( \sqrt[3]{a_1a_2^2}+\sqrt[3]{a_1a_2a_n}
+\sqrt[3]{a_1a_{n-1}a_n}+\sqrt[3]{a_{n-1}^2a_n} \right)
\end{align}
$$
where we can enter $a_1=t^6$, $a_2=(t+u)^6$, $a_{n-1}=(t+u+v)^6$,
$a_n=(t+u+v+w)^6$ with $t,u,v,w\ge0$ to obtain a polynomial $G(t,u,v,w)$ for the limit. Again, we expand and check the signs of the coefficients, and there is only one negative coefficient: $-t^4uw$. However, this is balanced out by coefficients $3t^4u^2$ and $3t^4w^2$, and the remaining terms are all non-negative. So $G(t,u,v,w)\ge0$ which proves that
$$
A(a)-B(a)
\ge \lim_{\epsilon\rightarrow0}A(a')-B(a')
= G(t,u,v,w)
\ge 0.
$$