Yes. There is a connection, one that is articulated in the following paper:
Your Classifier is Secretly an Energy Based Model and You Should Treat it Like One. Will Grathwohl, Kuan-Chieh Wang, Jörn-Henrik Jacobsen, David Duvenaud, Mohammad Norouzi, Kevin Swersky. ICLR 2020.
In particular, we can think of neural networks with softmax as inspired by an energy-based formulation, and in particular, by a Boltzmann distribution. Basically, we try to fit a Boltzmann distribution to the data, and then use that construct a classifier.
The paper suggests that we should think about neural networks as follows: we should think of them as a way of defining an energy-based model for estimating the joint probability $p(x,y)$, via the formula
$$p(x,y) = c \cdot e^{f(x)[y]},$$
where here $f(x)[y]$ represents the logit output for class $y$, when neural network $f$ is given input $x$, and $c$ is an appropriate normalizing constant. This is basically a way to parametrize the probability distribution $p(x,y)$ in a learnable way: the problem of learning a classifier is to learn a function $f$ so that $c \cdot e^{f(x)[y]}$ is a good approximation to the true joint distribution $p(x,y)$.
Notice that this formula means that we are approximating $p(x,y)$ by a Boltzmann distribution. In particular, the Boltzmann distribution would be $p(x,y) = c \cdot e^{-\epsilon_{x,y}/(kT)}$. Here we take the special case of $kT=1$. Moreover, we use the neural network to predict the appropriate $\epsilon_{x,y}$ values: in particular, we define $\epsilon_{x,y}$ to be the logit value produced by the neural network.
Given this definition, we can now use this to build a discriminative classifier. In particular, we can calculate
$$\begin{align*}
p(y|x) &= {p(x,y) \over p(x)}\\
&={p(x,y) \over \sum_{y'} p(x,y')}\\
&={c \cdot e^{f(x)[y]} \over \sum_{y'} c \cdot e^{f(x)[y']}}\\
&={e^{f(x)[y]} \over \sum_{y'} e^{f(x)[y']}}
\end{align*}$$
Looking closely, we see that the latter expression is exactly the softmax applied to the logits $f(x)$. Therefore, using this energy-based approach to fitting the distribution $p(x,y)$ yields exactly the same result as constructing a neural network that ends with a softmax layer. It's just two different ways to look at the same thing.
Hopefully this now makes clear the connection. See the paper above for more details and discussion.