This is essentially (a special case of) the statement of Poincaré's lemma. The relevant statement here goes as follows:
Let $U\subseteq \mathbb R^3$ be star-shaped, and let $\vec B$ be a smooth vector field on $U$ such that $\nabla \cdot \vec B = 0$. Then there exists a smooth vector field $\vec A$ on $U$ such that $\nabla \times \vec A = \vec B$.
Proof: The proof works by explicit construction. For the moment, let $U$ be star-shaped with respect to the origin. We then define
$$\vec A(\vec r) := \int_0^1 t \left( \vec B(t\vec r) \times \vec r\right) \mathrm dt$$
Crucially, note that we are assuming that for each $\vec r\in U$, we have that $t\vec r\in U$ for each $t\in[0,1]$. This is where the requirement that $U$ is star-shaped with respect to the origin comes in. Computing the curl yields
$$\nabla \times \vec A = \int_0^1 t \left(3 \vec B(t\vec r) - \underbrace{t\vec r \big(\nabla \cdot \vec B\big)(t\vec r)}_{=0} + t\big([\vec r \cdot \nabla]\vec B\big)(t\vec r) - \underbrace{\big(\vec B(t\vec r)\cdot\nabla \big)\vec r}_{=\vec B(t\vec r)}\right) \mathrm dt$$
$$=\int_0^1\left( 2t \vec B(t\vec r) + t^2 \big([\vec r\cdot \nabla]\vec B\big)(tr) \right) \mathrm dt = \int_0^1 \frac{d}{dt} \left(t^2 \vec B(t\vec r)\right) \mathrm dt = \vec B(\vec r)$$
as promised. If $\vec B$ is star-shaped with respect to some point $\vec R$ instead of the origin, the construction is straightforwardly modified:
$$\vec A(\vec r) := \int_0^1 t \left(\vec B\big(t\vec r+ (1-t)\vec R\big) \times \vec r\right)\mathrm dt$$
This produces the same result, with slightly messier algebra.
It's worth noting that Poincaré's lemma can be extended to all contractible domains. This proof is somewhat more annoying so I will omit it, but it follows essentially the same spirit; the salient difference is that where we consider the straight line segments $t\vec r$ and $t\vec r + (1-t)\vec R$ for $t\in[0,1]$, the more general proof allows for more general smooth curves.
Lastly, the implication $\nabla \times \vec F = 0 \implies \vec F = \nabla \chi$ is another specific case of Poincaré's lemma, and it is proved by construction in the same way. It's an excellent exercise in elementary vector calculus to figure out how to modify our prior construction for this case; the answer is given below.
$$\chi(\vec r) := \int_0^t \big[\vec F(t\vec r) \cdot \vec r \big] \mathrm dt$$
It should also be noted that a stronger statement is possible in this case - the implication holds even if $U$ is merely simply-connected. This can be proven from Stokes's theorem fairly straightforwardly, but this answer is already long enough.
For some additional intuition, consider the following counterexamples to Poincaré's lemma.
Let $\vec B = \frac{\alpha}{r^2}\hat r$ in spherical coordinates, on the domain $U=\mathbb R^3-\{\vec 0\}$, which would be the magnetic field of a magnetic monopole. Note that $\nabla \cdot \vec B = 0$ on $U$, but there is no $\vec A$ such that $\vec B = \nabla \times \vec A$. Poincaré's lemma does not apply here because $U$ is not contractible. However, one could cover $U$ with two overlapping but individually contractible domains $U_1$ and $U_2$, and on each we could find a vector potential $\vec A_1$ and $\vec A_2$ which are related by a gauge transformation.
Next, let $\vec B = \frac{\alpha}{r} \hat \varphi$ in cylindrical coordinates on the domain $\mathbb R^3$ minus the $\hat z$-axis. You may recognize as the magnetic field outside of a straight, infinitely long wire which lies along the $\hat z$-axis. Despite the fact that $\nabla \times \vec B = 0$ on $U$, this field cannot be written as the gradient of a scalar potential (feel free to try) on $U$, which is made obvious by the fact that its line integral along a closed loop containing the $\hat z$-axis is nonzero.
We can get close by letting $\psi = \alpha \varphi$, but $\varphi$ has a branch cut somewhere and therefore this does not constitute a valid vector potential on all of $U$.
Once again, the reason is that the expression written above applies only to the region outside of the wire, which is not contractible. If the wire has a finite diameter, then $\vec B$ will fail to be differentiable at the wire's surface which invalidates our assumptions. If the wire is infinitely thin, then $\vec B$ will fail to even exist at the origin, once again invalidating our assumptions.
Then $\vec{\nabla} \cdot \vec{B} = 0 \Rightarrow -\nabla^2 \phi = 0.$ However, this does not imply that $\phi \equiv 0$. So why can't magnetic fields be described in terms of some vector potential $\vec{A}$ but also some scalar potential $\phi$ that has $0$ Laplacian?
It doesn't, but that's not what you asked. $\nabla \cdot \vec B$ implies that there exists some $\vec A$ such that $\nabla \times \vec A= \vec B$ when the domain in question is star-shaped (or more generally, contractible). It does not imply that it is impossible to express $\vec B$ in any other way.
Indeed, as described above if additionally $\nabla \times \vec B = 0$ then we may locally write $\vec B = \nabla \psi$ for some scalar function $\psi$. This magnetic scalar potential is a frequently-used tool in magnetostatics.