I will start by motivating of the equations. You assume $B$ is the magnetostatic solution with a temporal modulation:
$$
\nabla\cdot B=0\\
\nabla\times B=\mu_0je^{-i\omega t}
$$
You want to solve for $E$ satisfying:
$$
\nabla\cdot E=0\\
\nabla\times E=-\partial_tB
$$
This approach is the quasistatic approximation of the Maxwell equations. It amounts to neglecting the displacement current, and can be viewed as the leading order term in a perturbation series in $c$ (the speed of light). If you’re interested in radiation or the far field, then this approximation is invalid. However if you’re interested in the near field or low frequency regime, then it is certainly useful. Quantitatively, it works as long as:
$$
\omega r\ll c
$$
Your problem is linear, inhomogeneous and you found an inhomogeneous solution:
$$
E = \frac{i\omega}{2\pi}\ln(r)e^{-i\omega t} e_z
$$
You can always add a homogeneous solution. Examples of such homogeneous solutions is given by your three additional terms:
$$
E = \frac{c_1}{r}e_r+\frac{c_2}{r}e_\phi+c_3 e_z
$$
For a physical interpretation of this indeterminacy, your field equations are valid only for $r>0$ due to the change in cylindrical coordinates. There are two main approaches to see the problem. The first one is to only use your equations, but then the domain is space minus the $z$ axis which is topologically non trivial. Another approach is that the RHS of the equations could contain distributions supported on the $z$ axis.
For Gauss’ law, this would mean some electric charges (or higher multipoles) on the $z$ axis which would not be visible in the equation for $r>0$. The term:
$$
E\propto\frac{1}{2\pi r}e_r
$$
is a contribution of a uniform monopole. Without the symmetry requirement, you would have found additional undetermined terms with $\phi,z$ dependence like the uniform dipole along $x$ contribution:
$$
E\propto\frac{1}{2\pi r^2}(\cos\phi e_r+\sin\phi e_\phi)
$$
In your case, there are no charges, so $c_1=0$.
For the curl, the equations are the same as in magnetostatics assuming zero divergence. You can therefore have terms corresponding to magnetic multipole on the axis. Your term:
$$
E\propto\frac{1}{2\pi r}e_\phi
$$
corresponds to the uniform moment along $z$ contribution.
Physically, those multipole terms would correspond to an additional magnetic field localized on the axis. This magnetic field must be generated by other currents according to magnetostatics (still in the quasistatic approximation). Your $c_2$ term corresponds to an infinitely thin solenoid about the $z$ axis. This is why $c_2=0$ in your example.
Thus your symmetry requirement is an efficient way of getting rid of unwanted multipolar fields that are not relevant for your situation. In the same vein, you can use more symmetries to get rid of $c_1,c_2$. Any plane passing by the $z$ axis are planes of symmetry, so $B$ is orthoradial. Furthermore, any plane perpendicular to the $z$ axis are planes of antisymmetry, so $E$ is along $z$.
Finally, your third term arises from logarithm. This divergence analogous to the log divergence of the electric field about a uniformly charged line. Similarly, it shows that a single line current is not physical. They come in pairs of opposing currents like in waveguides (this is a form of confinement). Therefore the leading order logarithms cancel out leaving only a $1/r$ decay at infinity. Continuing the electrostatic analogue, only the 2D dipole contribution survives, the monopole contribution is cancelled, so you’d expect:
$$
E\sim\frac{\mu_0I_0\omega d}{2\pi r}(\cos\phi e_r+\sin\phi e_\phi)
$$
with with the opposing line of current being at a distance $d$ along the $x$ direction with respect to the $z$ axis.
Thus the additive constant $c_3$ must be chosen in accordance with the other line’s constant so that they cancel and the field decays to zero at infinity. A pedestrian approach would be to realize that you need a characteristic length to make the argument of the logarithm dimensionless. In your original problem, there are none, so this is why it is ill defined. Adding another opposing wire introduces a new length: their relative distance.
The logarithmic divergence for large $r$ is an artefact of the quasistatic approximation. The full radiating field does decay at infinity. Using the displacement current, your equations become (setting $\mu_0=c=I_0=1$):
$$
-E_z’=i\omega B_\phi\\
\frac{1}{r}(rB_\phi)’=i\omega E_z
$$
Ampère’s law gives the “initial condition”:
$$
2\pi rB_r \xrightarrow{r\to 0} 1
$$
You usually choose the causal solution:
$$
\begin{align}
B_\phi &=-\frac{\omega}{4}H^{(1)}_1(r\omega)\\
&\sim \frac{1}{2\pi r} &&r\to 0\\
&\sim \frac{1}{2}\sqrt{\frac{\omega}{2\pi r}}e^{ir\omega}&& r\to\infty \\
E_z &= \frac{i\omega}{4}H^{(1)}_0(r\omega) \\
&\sim \frac{i\omega}{2\pi }\ln(r\omega) && r\to 0\\
&\sim \frac{i}{2}\sqrt{\frac{\omega}{2\pi r}}e^{ir\omega}&& r\to\infty
\end{align}
$$
Thus the fields do decay to zero, but the energy diverges linearly. Adding a line of opposite current will give finite energy.
In general, without the causality prescription, there is still an indeterminacy as you can always add homogeneous solutions (cylindrical waves) of the form:
$$
\begin{align}
B_\phi &= iCJ_1(r\omega)\\
E_z &= CJ_0(r\omega)
\end{align}
$$
with an undetermined $C$. This is the same as your indeterminate $c_3$ (or $r_0$). Causality (or equivalently Sommerfeld radiation conditions) allows you to solve the indeterminacy. It gives you the missing characteristic length $1/\omega$ with additional constants:
$$
E_z = \frac{i\omega}{2\pi }(\ln(r\omega)-\ln2+\gamma-i\pi/2)
$$
