Abstract

In this work, we study approximations of supercritical or suction vortices in tornadic flows and their contribution to tornadogenesis and tornado maintenance using self-avoiding walks on a cubic lattice. We extend the previous work on turbulence by A. Chorin and collaborators to approximate the statistical equilibrium quantities of vortex filaments on a cubic lattice when both an energy and a statistical temperature are involved. Our results confirm that supercritical (smooth, “straight”) vortices have the highest average energy and correspond to negative temperatures in this model. The lowest-energy configurations are folded up and “balled up” to a great extent. The results support A. Chorin’s findings that, in the context of supercritical vortices in a tornadic flow, when such high-energy vortices stretch, they need to fold and transfer energy to the surrounding flow, contributing to tornado maintenance or leading to its genesis. The computations are performed using a Markov Chain Monte Carlo approach with a simple sampling algorithm using local transformations that allow the results to be reliable over a wide range of statistical temperatures, unlike the originally used pivot algorithm that only performs well near infinite temperatures. Efficient ways to compute entropy are discussed and show that a system with supercritical vortices will increase entropy by having these vortices fold and transfer their energy to the surrounding flow.

Keywords

Tornadogenesis, Supercritical Vortices, Vortex Filaments, Negative Temperature, Kinetic Energy, Entropy, Statistical Mechanics, Equilibrium Statistics, Self-Avoiding Walks, Cubic Lattice, Monte-Carlo Techniques, Pivot Algorithm

Share and Cite:

1. Introduction

During the formation stage of a tornado, long narrow vortices often appear spontaneously in the region of tornado formation, then fold up and dissipate. These vortices also appear in other stages of the tornado’s existence. An example of such vortices in a strong tornado is shown schematically in Figure 1 [1] . As indicated in the figure, such vortices are called suction vortices. These vortices are straight and narrow, locally concentrating a lot of kinetic energy. Their size is often below that detectable by a radar and even by mobile radars. They are transient and often fold up and dissipate within one complete revolution around the main tornadic flow. Evidence for them can be seen in tornado videos or in the damage surveys done after the tornado has passed. In these surveys, such vortices leave behind tracks, where grass in lawns or plants in farmers’ fields have been ripped from the ground, indicating their high energy density [2] (see Figure 2), sometimes followed by debris dropped along the end of the tracks, due to dissipation of the vortices. Such strong, narrow vortices have been analyzed in [3] as so-called supercritical vortices. Their formation is related to the breakdown in the cyclostrophic balance, where the pressure-gradient force dominates, and the vortex collapses to a narrow filament; as the vortex collapses, the energy density increases and the entropy density decreases. This mirrors the behavior of negative-temperature vortices studied by A. Chorin in his work on turbulence [4] [5] [6] [7] [8] . In these works, the negative-temperature vortices are straight and as they transfer energy to the surrounding flow, they fold up and dissipate. The similarity between the behavior of supercritical vortices studied in [3] and negative-temperature vortices studied by Chorin is striking and is one of the motivations for this paper. Our results will support and illustrate this phenomenon as well.

Recent progress in radar observations and numerical modeling of tornadogenesis and maintenance has revealed the importance of small-scale vorticity’s contributions to development of the large-scale rotation. For example, Figure 11 in the radar observations study [9] shows a stream of small-scale vortices-flowing into the region of a developing tornado. Similarly, Figure 5 and Figure 6 in the numerical modeling study [10] show rivers of vorticity forming in response to surges in outflow. Specifically, it is stated on p. 3044: “There are distinct rivers of cyclonic (and anticyclonic) shear vorticity that originate from the base of downdrafts. Near-ground tornado-like vortices are fed by such cyclonic $\zeta$ rivers during their genesis.”

Another study references [9] [10] and the rivers of vorticity but argues for further study of the phenomenon [11] : “other examples may be available but are not widely recognized. However, these rivers are yet to be seen definitively in

multiple-Doppler-radar analyses, to the best of the authors’ knowledge, perhaps because they are too small scale to be detected (owing to the inherent smoothing during the analysis procedure) or because they occur only over a very shallow layer near the ground, where data collection can be difficult.”

The recent paper [12] considers the supercell environment, it discusses physical processes producing horizontal vorticity and mechanisms for re-orienting it into the vertical, leading to tornadogenesis. The paper suggests ways of using output from a numerical simulation to test for the origins of vorticity in the simulation, thus giving insight into the process of tornadogenesis.

Many storm chaser videos show small-scale vortices entering tornadoes. Examples of small-scale vortices can be seen, for example, in the video in [13] between the times of 2:20 and 2:50 minutes. Small vortices entering the tornado can be seen at two separate instances in the video in [14] , one at 0:35 to 0:40, and one at 0:52 to 2:14 minutes. One also sees the filament-like vortices folding up and dissipating as would be expected for negative-temperature filament-like vortices; this point is one of the main observations of this study.

Supercritical or suction vortices also play an important role in tornadogenesis and tornado maintenance. A deeper discussion of this idea appears in [15] [16] and numerical evidence is provided in the state-of-the-art simulation [17] . A snapshot of the simulated dynamical behavior of these vortices appears in Figure 3, in which intense, narrow, vertical vortices to the right of the developing or existing tornado move into the region where the tornado is forming or has formed. In the simulation, these vortices eventually fold and dissipate and transfer their energy to the surrounding flow. The physical existence of several

intense multiple vortices within the larger flow is revealed by the radar data of the May 31, 2013, El Reno, OK tornado obtained by a Doppler on Wheels mobile radar [18] . Analysis and discussion of such smaller, violent vortices within a large tornadic flow is provided in [19] . Finally, a whole hierarchy of vortices within vortices is discussed in [20] .

Energy and entropy serve as fundamental descriptors of any thermodynamic system with applications in physics, atmospheric science, biochemistry, biophysics, and more. A recent example, reaching as far as DNA genetic mutations, is [21] . Motivated by Chorin’s results, our goal is to reliably compute the energy and entropy of vortex filaments on a cubic lattice and use this knowledge in the context of supercritical suction vortices and their behavior. To this extent, we employ the model developed in [4] [5] [6] [7] [8] , in which a vortex filament is modeled as a self-avoiding walk (SAW) on a cubic lattice and its energy is readily computed. This model was further extended to vortex structures with Brownian cores and fractal cross sections [22] and to Brownian semimartingales [23] . Another model rigorously studies an ensemble of nearly parallel vortices, however, under the restriction that the vortex filaments cannot fold [24] , a restriction that is detrimental to our consideration of vortices that fold up and dissipate. The cubic lattice SAW model used by Chorin has yielded results in a narrow range of statistical temperatures due to the employment of a Markov chain Monte Carlo (MCMC) algorithm (the pivot algorithm [25] [26] ) that is well suited for polymers (infinite temperature or maximum entropy case), but fails to deliver reliable results when this is not the case. We utilize a more flexible algorithm that alleviates most of the problems experienced by the pivot algorithm. This algorithm allows us to compute equilibrium average energies that can be validated when exact values are known, and at the moment we do not have any indication that the results are significantly off in general. Having computed energies that appear reliable, we also propose an efficient way to compute the entropy of the system whose accuracy is mainly affected by how accurately the average energies have been computed.

2. Background Mathematics

Atmospheric flows can be modeled using Euler’s equations relating velocity, pressure, and mass density of the fluid and external body forces. Typical flows are usually incompressible, so the divergence of the velocity field is zero. Isentropic flows are nearly incompressible when one of two scenarios occurs: either the flow speeds or the local changes in flow speeds along streamlines are small compared to the speed of sound in the medium [27] . A numerical study of intense tornadic compressible and incompressible isentropic flows has shown little difference in results [28] .

With $u$ being the fluid’s velocity, $\rho$ its mass density, p its pressure, and $b$ an external body force, all functions of position $x\in {ℝ}^3$ and time $t\in ℝ$ , the governing equations for an incompressible fluid flow are

$\frac{Du}{Dt}=-\frac{1}{\rho }\nabla p+b,\nabla \cdot u=0,\frac{D\rho }{Dt}=0,$ (1)

where $\frac{Df}{Dt}=\frac{\partial f}{\partial t}+\left(u\cdot \nabla \right)f$ denotes the material derivative of a scalar function $f$ . When applied to a vector function, the operator $\frac{D}{Dt}$ applies component wise.

The vorticity field $\xi$ of the velocity field $u$ is given by the curl of $u$ , i.e., $\xi =\nabla \times u$ . By applying curl to the first equation in (1) one can obtain an equation for $\xi$ ,

$\frac{\partial \xi }{\partial t}=\nabla \times \left(u\times \xi \right)+\frac{1}{{\rho }^2}\nabla \rho \times \nabla p+\nabla \times b\mathrm{.}$

In this equation, the first term in the right-hand side corresponds to the “barotropic” generation of vorticity (capturing the advection, stretching, and tilting of the vertical vorticity; see, e.g., [29] ), while the second term corresponds to the “baroclinic” generation of vorticity, i.e., vorticity generation due to the misalignment of the gradients of mass density and pressure. The body force is often conservative (e.g., due to gravity), in which case $\nabla \times b=0$ .

In what follows, we will focus on flows in which vorticity is supported on long, narrow vortex tubes embedded in a larger irrotational flow. These may be baroclinic or barotropic in origin, and may demonstrate themselves, e.g., as the suction vortices discussed earlier [2] .

Given a sufficiently fast decaying vorticity field $\xi \left(x\right)$ in ${ℝ}^3$ , the system

$\xi =\nabla \times u\mathrm{,}\nabla \cdot u=0$

can be solved for the velocity field [30]

$u\left(x\right)=-\frac{1}{4\pi }{\displaystyle {\int }_{{ℝ}^3}}\frac{\left(x-x^{\prime }\right)\times \xi \left(x^{\prime }\right)}{{\left|x-x^{\prime }\right|}^3}\text{d}{x}^{\prime }\mathrm{.}$

This expression is known as the Biot-Savart law. From this, one has an expression for the kinetic energy of the flow [31]

$\begin{array}{c}E=\frac{1}{2}{\displaystyle {\int }_{{ℝ}^3}}\rho \left(x\right){\left|u\left(x\right)\right|}^2\text{d}x\\ =\frac{1}{8\pi }{\displaystyle {\iint }_{{ℝ}^3\times {ℝ}^3}}\rho \left(x\right)\frac{\xi \left(x\right)\cdot \xi \left(x^{\prime }\right)}{\left|x-x^{\prime }\right|}\text{d}{x}^{\prime }\text{d}x+\frac{1}{8\pi }{\displaystyle {\iint }_{{ℝ}^3\times {ℝ}^3}}\frac{\left(\nabla \rho \left(x\right)\times u\left(x\right)\right)\cdot \xi \left(x^{\prime }\right)}{\left|x-x^{\prime }\right|}\text{d}{x}^{\prime }\text{d}x.\end{array}$

For a homogeneous fluid the second term above vanishes and if density is rescaled so that $\rho \left(x\right)\equiv 1$ , the expression for the kinetic energy becomes [8]

$E=\frac{1}{8\pi }{\displaystyle {\iint }_{{ℝ}^3\times {ℝ}^3}}\frac{\xi \left(x\right)\cdot \xi \left(x^{\prime }\right)}{\left|x-x^{\prime }\right|}\text{d}{x}^{\prime }\text{d}x\mathrm{,}$ (2)

where we note, however, that the integral can be taken over $\text{supp}\xi \times \text{supp}\xi$ , so the integrals are evaluated over the vortex tubes only.

Following Chorin’s work, we now define an approximation to the energy (2) for an infinitely thin vortex filament, a “vortex line”. The assumption is that the vortex tube can be approximately divided into N narrow circular cylinders $I_i$ of equal length, for which we then have

$\begin{array}{c}E=\frac{1}{8\pi }{\displaystyle {\iint }_{{ℝ}^3\times {ℝ}^3}}\frac{\xi \left(x\right)\cdot \xi \left(x^{\prime }\right)}{\left|x-x^{\prime }\right|}\text{d}{x}^{\prime }\text{d}x\\ =\frac{1}{8\pi }{\displaystyle \underset{i}{\sum }}{\displaystyle \underset{j\ne i}{\sum }}{\displaystyle {\int }_{I_i}}{\displaystyle {\int }_{I_j}}\frac{\xi \left(x\right)\cdot \xi \left(x^{\prime }\right)}{\left|x-x^{\prime }\right|}\text{d}{x}^{\prime }\text{d}x+\frac{1}{8\pi }{\displaystyle \underset{i}{\sum }}{\displaystyle {\int }_{I_i}}{\displaystyle {\int }_{I_i}}\frac{\xi \left(x\right)\cdot \xi \left(x^{\prime }\right)}{\left|x-x^{\prime }\right|}\text{d}{x}^{\prime }\text{d}x,\end{array}$

where the first term in the last expression corresponds to interactions of one of the cylinders with the others and is referred to as interaction or exchange energy, while the second term is referred to as self-energy and gives the contribution to the kinetic energy arising from interactions of nearby points along the vortex filament. In Chorin’s work, the self-energy term is neglected. For discussions of this term see, e.g., [5] [6] [7] [8] .

Finally, for an infinitely thin vortex line, the exchange energy is approximated in the following way. Assume that the vortex consists of N linear segments $I_i$ of equal length, whose midpoints are denoted by $M_i$ . Since vorticity along a vortex line is parallel to it, let ${\xi }_i$ denote the “vorticity” vector on the segment $I_i$ , which is in the direction of $I_i$ and has magnitude equal to the length of the segment $I_i$ . Then the exchange energy is approximated by

$E=\frac{1}{8\pi }{\displaystyle \underset{i}{\sum }}{\displaystyle \underset{j\ne i}{\sum }}\frac{{\xi }_i\cdot {\xi }_j}{\left|M_i-M_j\right|}\mathrm{.}$ (3)

This scaling corresponds to the vortex having circulation equal to one [5] [6] . In the next section we will discuss the cubic lattice approximation, in which the individual segments $I_i$ connect two nearest neighbors on a cubic lattice, so only angles of 90˚ and 180˚ between neighboring segments will be allowed. Notice that due to the dot products in (3) the largest contribution to the energy comes from nearby segments oriented in the same direction, smallest (negative) contribution from nearby segments oriented in the opposite direction, and orthogonal segments contribute zero to the energy. That is, the largest energy will correspond to straight vortices, while smallest energy will likely correspond to very folded up vortices.

3. Vortex Filaments on a Cubic Lattice

To simplify the set of all possible configurations of a line vortex in 3D, we only consider those constrained to a cubic lattice ${\mathbb{Z}}^3$ with lattice constant one [4] [5] [6] [7] [8] . The vortex filament will then correspond to a self-avoiding random walk (SAW) on this lattice, a concept of great interest in the polymer and protein community (see, e.g., [32] ). Even with this simplifying assumption, the problem of studying various vortex configurations is intractable exactly, since the number of all possible SAWs of N segments appears to grow exponentially as a function of N [33] [34] [35] . The exact enumeration of all possible configurations has been achieved for only small values of N and the number has been increasing slowly with time (N = 26 in 2000 [33] and N = 36 in 2011 [35] ). In Table 1 we show the number of SAWs of N steps on a cubic lattice for $N=1,\cdots ,9$ , 36 (all exact), and 1000 (an estimate based on a best-fit formula [35] ).

In this paper we demonstrate an improvement of the computation of average energy and entropy compared to previous attempts, and thus fairly small values of N (≤1000) will be considered. These are small values compared to the computational state of the art, however the numbers of distinct SAWs, M, are still large, ranging from $M\approx {10}^{67}$ for $N=100$ to $M\approx {10}^{671}$ for $N=1000$ [34] [35] . As such, a statistical approach to our problem is required and an appropriate Monte Carlo technique with a suitable sampling procedure is usually employed.

In the next sections we describe the relevant statistical mechanics and the algorithm used to obtain accurate results.

4. Statistical Mechanics Background

We now briefly review relevant concepts from statistical mechanics. Given a positive integer N corresponding to the number of segments in a self-avoiding walk on a cubic lattice, we denote by M the total number of possible SAWs starting at the origin. We denote by $x_i$ , $i=1,\cdots ,M$ , the possible SAWs, or vortex configurations of length N. The vortex configurations $x_i$ represent individual states in the system ${\mathcal{S}}_N=\left\{x_i\mathrm{<mo>\; :\; </mo>}i=\mathrm{1,}\cdots \mathrm{,}M\right\}$ .

The energy of each configuration, a state energy, will be denoted by $E_i=E\left(x_i\right)$ with E defined in (3). As introduced in [6] [7] [8] and re-derived in [15] [16] , the probability $p_i$ that a system will be in state $x_i$ is given by the Boltzmann probability distribution

$p_i=\frac{{\text{e}}^{-\beta E_i}}{Z},Z={\displaystyle \underset{i=1}{\overset{M}{\sum }}}{\text{e}}^{-\beta E_i},$ (4)

where Z is called the (canonical) partition function and $\beta =\frac{1}{k_{B}T}$ is sometimes

referred to as “coldness” as it is proportional to the reciprocal of the (statistical) temperature, T. The Boltzmann constant, k_B, will be assumed equal to 1 in subsequent computations but we will use it in the general results in this section. Note that when $\beta =0$ all states are equally likely and $p_i=1/M$ for $i=1,\cdots ,M$ . This is sometimes referred to as the polymeric case as equally likely configurations are considered when modeling the behavior of polymers [32] .

Notice that the relationship $\beta =\frac{1}{k_{B}T}$ can be also used to define the temperature $T=\frac{1}{k_B\beta }$ . Consider β increasing from −∞ to +∞. When β increases from

−∞ to 0, this corresponds to T decreasing from 0 to −∞; when β increases from 0 to +∞, this corresponds to T decreasing from +∞ to 0. Put together, when one considers the temperature T, temperature first increases from $T=0$ to $T=+\infty$ ; considering the one-sided limits as $\beta \to 0$ , one can identify $T=+\infty$ with $T=-\infty$ and simply write $T=\infty$ . Then temperature further increases as T changes from $T=\infty$ through negative values towards 0. Thus negative temperatures are higher than positive temperatures in this sense. (This idea corresponds to traversing a circle obtained by transforming the real number line into a circle by identifying +∞ and −∞, both for β and T; this identification is a special case of the Möbius transformation of the complex plane.) Another explanation for why negative temperatures are higher than positive ones is based on energy (and entropy) considerations. For further details, see [15] [16] [36] .

Using the Boltzmann distribution (4), we can now define in a standard way the average of any function of the states $x_i\in {\mathcal{S}}_N$ for any finite β. For example, the average energy $\langle E\rangle$ of the system ${\mathcal{S}}_N$ at a given β (or temperature T) is given by

$\langle E\rangle ={\displaystyle \underset{i=1}{\overset{M}{\sum }}}{p}_i{E}_i.$ (5)

One of the goals of this paper is to reliably approximate $\langle E\rangle$ for any finite β, which would extend the results of Chorin (see, e.g., [6] [7] [8] ) that appear reliable only for a small interval of values of β containing zero.

The Gibbs entropy of the system ${\mathcal{S}}_N$ at a given β is given by

$S=-k_B{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{p}_i\log p_i,$ (6)

which can also be viewed as an average quantity, since one can write

$S={\displaystyle \underset{i=1}{\overset{M}{\sum }}}{p}_i\left(-k_B\log p_i\right)={\displaystyle \underset{i=1}{\overset{M}{\sum }}}{p}_i{S}_i$

and view $S_i=-k_B\log p_i$ as an “entropy” of the state $x_i$ . Also note that if $\beta =0$ , then, since $p_i=1/M$ , we have $S=k_B\log M$ , which corresponds to Boltzmann’s definition of entropy.

Taking a partial derivative of the logarithm of the partition function in (4) with respect to $\beta$ , we obtain

$-\frac{\partial \log Z}{\partial \beta }=\frac{1}{Z}{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{E}_i{\text{e}}^{-\beta E_i}=\langle E\rangle .$ (7)

Using (4), expression (6) can be written as

$S={\displaystyle \underset{i=1}{\overset{M}{\sum }}}{p}_i\left(k_B\beta E_i+k_B\log Z\right)=k_B\beta \langle E\rangle +k_B\log Z,$

and differentiating it with respect to $\langle E\rangle$ and using (7) gives the well-known result

$\frac{\partial S}{\partial \langle E\rangle }=k_B\frac{\partial \beta }{\partial \langle E\rangle }\langle E\rangle +k_B\beta +k_B\frac{\partial \log Z}{\partial \beta }\frac{\partial \beta }{\partial \langle E\rangle }=k_B\beta =\frac{1}{T}.$ (8)

Finally, the Helmholtz free energy of the system is

$F=\langle E\rangle -TS={\displaystyle \underset{i=1}{\overset{M}{\sum }}}{p}_i\left(E_i+k_{B}T\log p_i\right),$

where, in view of (4) and the relationship between β and T, the term in parentheses equals $-k_{B}T\log Z$ , i.e., a constant independent of i, and hence the expression for the Helmholtz free energy gives

$F=\langle E\rangle -TS=E_i+k_{B}T\log p_i\text{for}\text{any}i=\mathrm{1,}\cdots \mathrm{,}M\mathrm{.}$ (9)

In particular, this means that the entropy (and also the Helmholtz free energy and the partition function) can be computed from β, $\langle E\rangle$ , $E_i$ , and $p_i$ for any i provided these values are known since (9) can be rewritten as

$S=k_B\left(\beta \langle E\rangle -\beta E_i-\log p_i\right)\text{for}\text{any}i=\mathrm{1,}\cdots \mathrm{,}M\mathrm{.}$ (10)

Equations (8) and (10) will be used later in Section 7, where we discuss approximating the entropy of the vortex filament system.

5. The MCMC Algorithm

In order to reliably compute various statistical quantities of interest, such as the average energy (5) or the entropy (6), one usually employs Markov chain Monte Carlo (MCMC) techniques, since the space ${\mathcal{S}}_N$ of possible configurations of the vortex filament of length N is too large for direct enumerations. As mentioned in Section 3, complete enumerations have been achieved only for small values of N up to date [33] [35] .

Here we briefly review the Metropolis-Hastings algorithm that we use in this work. We will consider two configurations, x and x’, their energies E and E’ given by (3), and their probabilities $P\left(x\right)$ and $P\left(x^{\prime }\right)$ given by (4). We define $g\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)$ to be the probability of proposing a transition from configuration x to configuration x’; $A\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)$ to be the probability of accepting a transition from x to x’; and $P\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)$ to be the probability of transitioning from x to x'. Note that $P\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)=A\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)g\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)$ and that we can easily compute the ratio $P\left(x^{\prime }\right)/P\left(x\right)={\text{e}}^{-\beta \left(E^{\prime }-E\right)}={\text{e}}^{-\beta \Delta E}$ since the difficult-to-compute partition function Z cancels out.

The detailed balance condition states that at equilibrium, each elementary process is in equilibrium with its reverse process, so we require

$P\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)P\left(x\right)=P\left(x\mathrm{<mo>\; |\; </mo>}{x}^{\prime }\right)P\left(x^{\prime }\right)\mathrm{.}$

Since $P\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)=A\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)g\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)$ , the detailed balance condition can be rewritten as

$\frac{A\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)}{A\left(x\mathrm{<mo>\; |\; </mo>}{x}^{\prime }\right)}=\frac{P\left(x^{\prime }\right)}{P\left(x\right)}\cdot \frac{g\left(x\mathrm{<mo>\; |\; </mo>}{x}^{\prime }\right)}{g\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)}={\text{e}}^{-\beta \Delta E}\cdot \frac{g\left(x\mathrm{<mo>\; |\; </mo>}{x}^{\prime }\right)}{g\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)}\mathrm{.}$

A typical way to satisfy this condition is to choose

$A\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)=\min \left\{\mathrm{1,}{\text{e}}^{-\beta \Delta E}\cdot \frac{g\left(x\mathrm{<mo>\; |\; </mo>}{x}^{\prime }\right)}{g\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)}\right\}\mathrm{,}$ (11)

and an alternative choice is [37]

$A\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)=\frac{1}{1+{\text{e}}^{\beta \Delta E}}\cdot \sqrt{\frac{g\left(x\mathrm{<mo>\; |\; </mo>}{x}^{\prime }\right)}{g\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)}}\mathrm{.}$ (12)

Note that sometimes $g\left(x^{\prime }\mathrm{<mo>\; |\; </mo>}x\right)\ne g\left(x\mathrm{<mo>\; |\; </mo>}{x}^{\prime }\right)$ and the ratio of the proposal probabilities has to be computed and retained in the algorithm. Either of the two acceptance probabilities (11) or (12) is then compared to a random number picked uniformly from $\left(\mathrm{0,1}\right)$ and based on the comparison the proposed transition is accepted or not. In this way one is able to generate a sequence of configurations $\left\{x_0\mathrm{,}{x}_1\mathrm{,}\cdots \right\}\subset {\mathcal{S}}_N$ in which each consecutive configuration is either the one resulting from an accepted transition, or a repeat of the previous configuration if the transition was rejected. The sequence can then be used to approximate the underlying probability distribution and its properties, and a good set of possible transitions is essential for fast and reliable performance.

Another requirement for the MCMC algorithm to work properly in sampling the configuration space is that it is ergodic, i.e., it is able to propose a finite sequence of transitions from any configuration $x_1\subset {\mathcal{S}}_N$ to any other configuration $x_2\subset {\mathcal{S}}_N$ .

An algorithm for proposing new configurations from old ones used in the previous works [5] [6] [7] was the pivot algorithm designed to generate a sequence of “effectively independent” SAWs. It was proposed by Lal in 1969 [25] and analyzed and popularized by Madras and Sokal in 1988 [26] . In this algorithm, given an N-step SAW, one randomly (usually uniformly) picks a point on it as a “pivot” and applies a transformation from the symmetry group of the cubic lattice to the points subsequent to the pivot, using the pivot point as the origin. The part of the SAW before the pivot stay intact. If the resulting walk is self-avoiding, it is accepted as the new element of the sequence; otherwise it is rejected and the original walk is repeated in the sequence. The symmetry group of the cubic lattice is the octahedral group O_h with 48 elements including the identity. The algorithm is ergodic and satisfies the condition of detailed balance [26] . In the works that motivated our paper [5] [6] [7] only a small range of values of β near 0 (or large in absolute value temperatures T) was explored.

In the context of our problem the pivot algorithm will not perform well in many scenarios due to the involvement of the vortex filament energy (3) in (4) and (5) and due to the fact that the pivot transformations are not local—in any given step all of the points after the pivot are expected to change position leading to a possibly large change in the filament energy. As commented in [6] , “For values of $\left|T\right|$ between 0.4 and 1 […], the results are unreliable but suggestive, and for $\left|T\right|<0.4$ they are completely unreliable.” This leaves only a small range of β values close to 0 where the authors in [6] feel confident in their results and leaves much room for improvement. Therefore, different algorithms to compute the average energy (5) more accurately across the whole spectrum of temperatures are needed.

Various modifications of the pivot algorithm have been considered in the literature for various reasons. For example, additional possible moves useful in modeling polymers are discussed in [32] and include “end flips,” “corner flips,” and “crankshaft moves.” Other modifications to speed up the algorithm have been described, for example, in [38] [39] . Algorithms for random walks with fixed endpoints have also been proposed [40] . More recently, other approaches have appeared in the polymer literature that aim to address finding the minimum energy configurations. They include “pull moves” [41] , “re-bridging moves” [42] , and “connectivity changing moves” [43] and appear to work well in the context of modeling the folding of HP polymers [44] [45] .

6. The Localized Transformations Algorithm

The rigidity of the pivot algorithm for problems that involve energy and the Boltzmann distribution away from the polymeric case ( $\beta =0$ or $T=\infty$ ) implies that a different approach is needed. In cases with $\left|\beta \right|\gg 0$ the transformations need to be local enough so that much of the shape of the current SAW remains the same, thus preserving much of the SAWs energy, and only a small part should be changed.

We have found that a very simple algorithm that we describe below works well for this problem across a wide spectrum of the values of β. The algorithm will use two types of transformations: one designed specifically to help with straightening out non-straight configurations at negative values of β, and one designed to help with compressing configurations in order to lower their energy at positive values of β. We next describe these two types of transformations and later we comment on their relationship to others existing in the literature.

For large (in absolute value) negative values of β, the pivot algorithm is very unlikely to straighten out a configuration with all but one segment in one direction (or similar, more complicated configurations), since in order to remove the “kink,” the whole part after the kink has to be temporarily transformed as well, resulting in a significant change in energy. An easy solution is to simply modify the one non-aligned segment while keeping the directions of the other segments the same. This still modifies the positions of the points after the kink, but it results in a small change in the filament energy because the relative positions of the points stay the same. We generalize this idea to the following algorithm for an N-step SAW:

First Reconstruction Algorithm:

1) Randomly select a number n_steps between 1 and N of consecutive steps to be modified.

2) Randomly select a subwalk of length n_steps to be modified.

3) Randomly reconstruct the selected subwalk.

4) Connect the reconstructed subwalk to the remaining one or two pieces by translation.

This algorithm is clearly ergodic since the reconstruction can be done on the whole SAW. Inverses of transformations are included and selected with equal probabilities, so the algorithm satisfies the detailed balance condition as well when the usual Monte Carlo acceptance probability is used. Practical considerations here include the choice of probability distributions in steps 1) and 2) (we used an exponential distribution in 1) to emphasize shorter reconstructions and a uniform distribution in step 2) and how much care should go into step 3) (e.g., random reconstruction vs. constructing a true SAW; we used a random reconstruction). It is easy to see that the configurations with one misaligned segment are easy to transform into a straight configuration with little energy change.

For large positive values of β, the filament configurations are expected to have small energy and be folded up and “balled up” to some degree in order to lower their energy. The first reconstruction algorithm is not an ideal candidate here due to the fact that the part of the filament after the reconstructed subwalk is likely to be translated during the reconnection, likely leading to self-intersections of already balled up configurations. Motivated by transformations proposed for configurations with fixed endpoints [40] , we propose a second reconstruction algorithm that reconstructs a subwalk while keeping its endpoints fixed:

Second Reconstruction Algorithm:

1) Choose a number K between 2 and N for maximum number of allowed reconstructed steps.

2) Randomly select a number n_steps between 2 and K of consecutive steps to be modified.

3) Randomly select a subwalk of length n_steps to be modified.

4) Permute the selected steps to generate a new subwalk.

5) Connect the reconstructed subwalk to the remaining one or two pieces.

Combining the two reconstruction algorithms results in an algorithm that is ergodic due to the ergodicity of the first algorithm. The second algorithm also clearly satisfies the detailed balance condition. Practical considerations here include the following. The number K can be chosen arbitrarily, but we found that $K\approx N^{2/3}$ not only roughly corresponds to the number of surface points in a perfectly “balled up” configuration, it also appears to lead to reasonably small minimum energy values. For the probability distributions in 1) and 2), we again used an exponential and a uniform distribution, respectively. The permutation in 3) could be a random permutation; we found that choosing a (uniformly) random cyclic permutation led to good results. In general, even though we never attempted to directly compute minimum energies, we found that our lowest-energy configurations in the range $-100\le \beta \le 100$ had energies below or within a fraction of a percent of the predicted minimum values as shown in Section 7.5.

We note that the first algorithm includes many of the standard simple transformations mentioned earlier, like the end flips, corner flips, and crankshaft moves. Some of the pull moves from [41] designed to find minimum-energy configurations are included in the first algorithm as well. Transformations that change the length of the SAW, such as those in the BFACF algorithm [46] [47] , are not suitable for our problem and are not part of this algorithm. For the same reason, the transformations of the join-and-cut algorithm [48] are not included either. This algorithm does not include any of the re-bridging moves [42] used in the modeling of HP polymer folding [44] [45] .

It is worth pointing out that an alternative approach to our work exists and consists in approximating the density of states using Wang-Landau sampling [49] , as done, for example, in [44] [45] . Such an approach, if successful in the context of our problem, would then allow for the computation of average energies across the whole temperature spectrum while avoiding the possibility of being “trapped” in some states due to the energy considerations. The utility of this approach is unclear in our case with an almost continuous energy spectrum, since in the context of HP polymers there are only a small number of energy levels to be explored. We hope to investigate this approach in the future and compare to our current results.

7. Numerical Results

In this section we present numerical results for the computation of the average energy $\langle E\rangle$ and the entropy S as a function of the length of the vortex filament, N, and the inverse temperature, β. These results extend the results presented in [6] [7] [8] (and the references within) to cover a much broader range of values of β than done previously. One of the goals is to demonstrate the improved quality of the results, not to extend them to the largest possible lengths of the filaments. Thus our results extend only up to $N=1000$ , but the range of values of β, for which we believe our results are reliable, is significantly extended compared to [6] [7] [8] .

In the individual subsections we present and discuss the validation of the algorithm in cases where the exact results are known, results for the average energy (5) for various lengths of the filaments, we also discuss and compute the entropy of the system, and we present some observations related to the minimum energy configurations, even though that is not the main focus of this paper.

7.1. Validation of the Algorithm

Without knowing the exact values of the average energy (5), it is impossible to assess the accuracy of the computed results. Due to the growth rate of the number of SAWs as a function of the number of steps N (see Table 1), it is only realistic to compare computed results to exact values for small values of N. To this end we completely enumerated all SAWs for $N=\mathrm{3,}\cdots \mathrm{,9}$ and obtained exact values of $\langle E\rangle$ for various values of $-100\le \beta \le 100$ .

To demonstrate the capabilities of the MCMC algorithm, we perform the following test. Given N and β, start with a straight configuration but first perform 10N iterations with $\beta =0$ to get away from it. Next, first perform 100,000 iterations with the given β to forget history and then perform 200,000 averaging iterations to approximate $\langle E\rangle$ for that β. Repeat with various values of β in the interval of interest. The computed results, together with the exact values, are shown in Figure 4. The computed approximations of $\langle E\rangle$ for each integer β are shown as dots, and the exact values, computed with increments $\Delta \beta =0.1$ , are shown as curves. Notice how the computed results agree well with the exact values, showing that the algorithm is flexible enough to handle both positive and negative values of β. We note that the pivot algorithm will not recover the straight configurations for N = 7, 8, and 9 (not shown here) even with 300,000 iterations, and there are also signs of struggle to approximate well the lowest-energy configurations (seen in the bottom left part of the plot).

7.2. Average Energy Computations for $N=100,\cdots ,1000$

In this section, we present the numerical results for computing the average energy, $\langle E\rangle$ , given in (5), using the proposed algorithm with localized transformations that comprises of the first and second reconstruction algorithms described in the previous section.

The implementation details are as follows. As in [6] , our Monte Carlo simulations are started from straight filaments of N steps, which, for a given N, are also the largest-energy configurations with energy

$E_{\max }\left(N\right)=\frac{1}{4\pi }\left[N{\displaystyle \underset{k=1}{\overset{N-1}{\sum }}}\frac{1}{k}-\left(N-1\right)\right].$ (13)

In light of (4), the straight configurations are the most likely configurations for negative temperatures ( $\beta <0$ ), and we expect $\langle E\rangle \approx E_{\text{max}}$ when β is large (in absolute value) negative. Therefore, we start our computation at one such β ( ${\beta }_{\text{first}}=-20$ for the results below) and compute an approximation of the average energy for this value of β using the straight configuration as a starting point. We then increase β by a small increment ( $\Delta \beta =0.5$ in the results below), use the last configuration from the computation with the previous value of β as a starting point, and compute $\langle E\rangle$ with the new β. We continue this way until a stopping value of β has been achieved ( ${\beta }_{\text{last}}=100$ in the results below).

The average energy for each value of β is computed by first allowing for a number of “burn-in” transformations to forget recent history and then by performing a different number of averaging transformations to compute $\langle E\rangle$ . The numbers of iterations depend on N and β and vary between 10,000 and 55,000 for burn-in and between 500,000 and 1,500,000 for averaging.

During each averaging iteration, either the first or the second reconstruction algorithm is selected randomly with probability 1/2, which could conceivably be modified based on the value of β. For the first reconstruction algorithm a ratio $r_1\in \left(\mathrm{0,1}\right]$ is passed in and a geometric probability distribution with the ratio $r_1$ is constructed for the numbers 1 through N. (If $r_1=1$ , the distribution is uniform.) A random number n_steps is drawn from the set $\left\{\mathrm{1,}\cdots \mathrm{,}N\right\}$ according to this probability distribution. Clearly, for $r_1<1$ it is more likely to reconstruct shorter subwalks than longer ones. Finally, a starting point of the subwalk to be reconstructed is chosen randomly uniformly from possible candidates. For the second reconstruction algorithm, a similar process is followed. With a ratio $r_2\in \left(\mathrm{0,1}\right]$ , first a geometric probability distribution with the ratio $r_2$ is constructed for the numbers 2 through K, where K is the integer nearest to N^2/3. Then a random number n_steps is drawn from the set $\left\{\mathrm{2,}\cdots \mathrm{,}K\right\}$ according to this probability distribution. Finally, a starting point of the subwalk to be reconstructed is chosen randomly uniformly from possible candidates. At this point, a reconstruction of the subwalk is performed according to steps 3) - 4) in the relevant reconstruction algorithm and the new filament is tested for self-avoidance and for acceptance in the MCMC algorithm. If it is self-avoiding and accepted using the Boltzmann probability criterion, the newly constructed filament becomes the new filament in the Markov chain; otherwise the current one is repeated.

In Figure 5 we show the computed results with $N=100$ through $N=1000$ . We note that in this and subsequent figures β intentionally runs from positive values on the left to negative values on the right so that temperature, $T=1/\left(k_B\beta \right)$ , increases from left to right. For each value of N, the computation was repeated six times and averaged results are graphed. For each N, the average energy increases with temperature and, as expected, it levels off at its maximum value

when $\beta <0$ and is large in absolute value. Specifically, for $\beta \in \left[-\mathrm{20,}-10\right)$ the average energies are near their maximum values (13), corresponding to, on average, straight configurations. As temperature is lowered in the computation, average energy decreases as expected in light of (4). Due to the finite lattice spacing, for each N there is a minimum energy a filament of N steps can achieve, so it is not surprising that the average energies level off as $\beta \to \infty$ . This is demonstrated in Figure 5 and its zoomed-in version in Figure 6. We do not have an explicit expression for the minimum energy similar to (13), but we provide some insights in Section 7.5. Notice that from (13) the maximum energies grow like $N\log N$ , which is illustrated in Figure 6, but the computed results suggest a linear decay in the minimum energies. In Figure 6 we also display error bars for one standard error of the averaged results to illustrate that, not surprisingly, there is more variability in the computed results closer to the infinite temperatures, but very little variability for β large positive. The same is true for β large negative, but that observation is trivial and hence not displayed here.

An important point to be emphasized here is that the vortices start folding when β (and the temperature) is still negative. Note first that for $\beta <0$ the straight configuration has the largest probability (4) of occurring, but as β approaches 0, the probability density is becoming more and more uniform. Therefore, the average energy decreases as seen in Figure 5 and a typical configuration with energy near the average value is not straight. To quantify such configurations, one can imagine fitting them into a smallest possible rectangular prism on the cubic lattice and measuring its longest dimension. This number will be between roughly N^1/3 and N for a filament of length N, or, when expressed relative to the maximal length, between N^−2/3 and 1. In Figure 7 we display these relative longest dimensions for $N=100$ and $N=200$ over the range $-12\le \beta \le 0$ .

Each displayed result is a mean of 6 configurations that are within 1% of the mean energy for each β. The error bars show ±1 standard deviation of the 6 results. We can see how the average lengths of the filaments decrease as $\beta \to 0^{-}$ , corresponding to the folding of the filaments.

7.3. Entropy Approximations

In Chorin’s work [6] [8] , entropy of the system at various temperatures is estimated using an algorithm based on an earlier work by Meirovitch [50] . The idea is to enumerate all possible very short SAWs (150 SAWs of length 3 are used in [6] ) and then use their relative frequencies in a Monte Carlo run to approximate the probabilities needed to evaluate an entropy expression similar to (6). A result for entropy per unit length of the ensemble of vortex filaments, S/N, as a function of $-1.5\le \beta \le 2.5$ , is shown in Figure 8, where $N=351$ is used (reproduced from [6] ). It is observed that “S has a maximum at $T=\infty$ ( $\beta =0$ ), as expected. The slope of S is much smaller on the positive T side than on the negative T side, as can be expected from the larger values of E for $T<0$ and from the relation $T^{-1}=\partial S/\partial E$ . Further, note that S/N varies little with N in the range where the calculation can be trusted, and thus S increases with N. The larger the filament, the larger its entropy.”

In order to validate these results and extend them to a larger range of β, we used as benchmarks the cases for which we had the complete enumerations ( $N=3,\cdots ,9$ ). The exact entropies for these cases are shown in Figure 9. As discussed in Section 4, for cases where exact enumeration is possible the entropy for $\beta =0$ is easy to evaluate as $S=\log M$ , where M is the number of SAWs of a given length. Self-avoiding walks on a cubic lattice have been enumerated for lengths up to $N=36$ [35] , and this data can be used to construct approximate expressions $\tilde{M}$ for M as a function of N for $N>36$ . While the reference [35] is missing one parameter value ( $c_1$ ), its arXiv.org version [51] provides a slightly different, but completely described formula to approximate M. For completeness, we reproduce it here,

$M\approx \tilde{M}=A{\mu }^N{N}^{\theta }\left(1+c{N}^{-\Delta }+k{\left(-1\right)}^N{N}^{-\alpha }\right)\mathrm{,}$ (14)

where $A=1.1951966888$ , $\mu =4.6840041570$ , $\theta =0.1597395125$ , $c=0.1227360755$ , $\Delta =1.4315024046$ , $k=-0.0619076482$ , and $\alpha =1.8985141134$ . We note that the relative errors resulting from this formula for $N=\mathrm{14,}\cdots \mathrm{,36}$ are all below 2 × 10⁻⁶; this results in absolute errors for the entropy, $\log \tilde{M}-\log M$ , to be of the same order. Consequently, the resulting approximation

$S/N=\left(\log M\right)/N\approx \left(\log \tilde{M}\right)/N$ (15)

is expected to be quite good and perhaps improve for longer filaments due to the N in the denominator.

Using (14) and (15) to approximate S/N for $N=351$ , we obtain $S/N\approx 1.54733$ , which appears significantly below the maximum value of about 1.67 estimated from Figure 8. In fact, as $N\to \infty$ , the quantity $\left(\log \tilde{M}\right)/N$ eventually monotonically decreases and approaches $\log \mu \approx 1.54415$ . Note that even for $N=9$ we have $S/N=\left(\log 1853886\right)/9\approx 1.60364$ . Clearly, a better algorithm for computing entropy than the one used in [6] is needed.

Improvements of the Meirovitch algorithm have appeared in the literature more recently and have been applied to magnetic systems, polymers, peptides, and liquids (see, e.g., [52] [53] [54] and references within). A possible promising approach would be to use a variant named the hypothetical scanning Monte Carlo (HSMC) method, which approximates the probability, $p\left(x\right)$ , of a given filament configuration, x, of length N in the configuration space ${\mathcal{S}}_N$ at a given $\beta$ . We provide here a brief sketch of the algorithm, which is described more fully in the references.

The approximation to $p\left(x\right)$ , denoted by $p^{\text{HSMC}}\left(x\right)$ , is computed as a product of conditional “transition” probabilities that the kth step of the filament is in the given position as in x, given that the previous $k-1$ steps are fixed as in x. To generate these probabilities, we would keep the first $k-1$ steps fixed and apply our algorithm as used to compute $\langle E\rangle$ by allowing reconstructions only between the nodes k through N. The transition probability would then be approximated by the ratio of the filaments that agree with the kth step of x versus all filaments generated in the Monte Carlo run. Clearly, this process would be computationally intensive, since $N-1$ Monte Carlo runs (progressively shorter and shorter) would have to be performed to compute a single $p^{\text{HSMC}}\left(x\right)$ . On the other hand, $p\left(x\right)$ could be approximated arbitrarily closely by allowing sufficiently long simulations if a good sampling algorithm is used.

With the computed $p^{\text{HSMC}}\left(x\right)$ the Helmholtz free energy, F, could be approximated via (9),

$F\approx F^{\text{HSMC}}=E\left(x\right)+\frac{1}{\beta }\log p^{\text{HSMC}}\left(x\right)\mathrm{,}$

where $E\left(x\right)$ is the easily computed energy of the filament x. Subsequently, the entropy could be approximated through the relationship $F=\langle E\rangle -TS$ if the average energy has been computed, as shown in (10),

$S\approx S^{\text{HSMC}}=k_B\beta \left(\langle E\rangle -F^{\text{HSMC}}\right)=k_B\left(\beta \langle E\rangle -\beta E\left(x\right)-\log p^{\text{HSMC}}\left(x\right)\right)\mathrm{.}$

It follows that the absolute error in the entropy computation, $\left|S-S^{\text{HSMC}}\right|$ , will be affected both by the accuracy of the HSMC computation of $p^{\text{HSMC}}\left(x\right)$ as well as the accuracy of the computed $\langle E\rangle$ ; more precisely, it will be roughly bounded by the sum of the relative error in approximating $p\left(x\right)$ and $\left|\beta \right|$ times the absolute error in approximating $\langle E\rangle$ (all multiplied by $k_B$ ).

Note that to compute the entropy of the system of filaments of length N for many possible values of β (such as shown in Figure 9 for example) requires $N-1$ Monte Carlo simulations for each value of β and thus a significant amount of CPU time. We have performed such computations for $N=\mathrm{3,}\cdots \mathrm{,9}$ for validation purposes, but since in the next section we will present an alternative, and much more efficient, approach, we will not present the results here.

7.4. Entropy Computations Based on Energy Values

A very efficient alternative to computing the entropy can be used if good approximations to $\langle E\rangle$ have been computed across an interval of β values and a single value for entropy is known. Recall the relationship (8), restated here for convenience:

$\frac{\partial S}{\partial \langle E\rangle }=k_B\beta \mathrm{.}$ (16)

If the β interval is partitioned with values ${\beta }_i$ , the corresponding (average) energies are denoted by $E_i$ , and the (sought) entropies by $S_i$ , then, provided S is smooth enough, (16) can be discretized as

$\frac{S_{i+1}-S_i}{E_{i+1}-E_i}=k_B{\beta }_{i+1/2}+\frac{{\left(\Delta {\beta }_i\right)}^3}{48}{\frac{{\partial }^{3}S}{\partial {\langle E\rangle }^3}|}_{\langle E\rangle =\tilde{e}}$

for some $\tilde{e}$ between $E_i$ and $E_{i+1}$ , where ${\beta }_{i+1/2}=\left({\beta }_{i+1}+{\beta }_i\right)/2$ and $\Delta {\beta }_i=\left|{\beta }_{i+1}-{\beta }_i\right|$ . This then leads to the following Euler-like method for approximating the entropy,

$S_{i+1}=S_i+k_B{\beta }_{i+1/2}\left(E_{i+1}-E_i\right)\mathrm{,}$ (17)

provided at least one of the values $S_i$ is known to start the recursion. Since, as discussed above, the entropy when $\beta =0$ can be computed for small N or approximated for larger N by $k_B\log \tilde{M}$ using (14), a starting value is readily available. Alternatively, one can use the HSMC algorithm to approximate the entropy for any one value of β, for example $\beta =0$ .

To validate this approach, we first apply it to the cases with $N=3,\cdots ,9$ and with exact (average) energy values provided at equidistant values of β with $\Delta \beta =0.5$ . The results are shown in Figure 10. The exact value of entropy at $\beta =0$ is taken as $\log M$ with M from Table 1. The Euler algorithm (17) is applied to compute the remaining entropy values. As shown in the figure, the computed values (large dots) agree well with the exact values (underlying curves). We also see that even after 100 steps in each direction away from $\beta =0$ the agreement is excellent.

For $N=100,200,\cdots ,1000$ , we use (14) and (15) to get an approximation for the entropy at $\beta =0$ , and apply algorithm (17) with the precomputed average energy values. For the values of $\langle E\rangle$ shown in Figure 5, the results for the scaled entropy, S/N, are shown in Figure 11. Here, $\Delta \beta =0.5$ as in the validation case shown in Figure 10. As expected, the scaled entropy has a maximum value when $\beta =0$ and the rate of increase for $\beta >0$ is smaller (in absolute value) than the rate of decrease for $\beta <0$ . Notice how the proximity of the ten curves suggests that entropy S increases linearly with the length of the filament, N. We note that as $\beta \to -\infty$ , we have $S\to \log 6$ since for each N there are six straight configurations with highest probability, each probability approaching 1/6. Therefore, for large enough (in absolute value) negative values of β we have ${\lim }_{N\to \infty }S/N=0$ as illustrated in Figure 11. On the other side, as $\beta \to \infty$ , the

fact that S/N appears to be approaching a finite limit suggests that the number of energy minimizing configurations grows exponentially as $N\to \infty$ or that an exponentially growing number of configurations can be found in the vicinity of the minimum-energy configurations.

To make a comparison to the results in [6] and shown in Figure 8, in Figure 12 we show a zoomed-in version of the computed entropies using an interval similar to that in Figure 8. Notice that qualitatively the results are similar. The results in Figure 8 are for $N=351$ , so our results for this N would appear between the third and the fourth curves from the top. Quantitatively we notice differences in both the slopes on the two sides of the maximum, as well as the maximum value itself.

7.5. Minimum Energy Results

Even though we do not have a direct comparison to exact values of the average energy (5) for large values of N, it appears that our algorithm allows for accurate approximation of the average energies for negative values of β and for values close to 0. How accurate the computed energies for large positive values of β are is less clear. It follows from (4) that for such values the probabilistically preferred configurations are those with lowest energies, but it is not clear what the lowest possible energy is. This is in contrast with the largest possible energy value that is clearly associated with the straight filament and its value is easily computable via (13).

Intuitively, the energy expression (3) suggests that the lowest energy will be associated with a filament folded in such a way that individual segments line up close to each other in an antiparallel way, resulting in large negative contributions to the energy, and also in such a way that there are many right angles

which contribute zero energy. From a fluid mechanics point of view this is also intuitive since antiparallel vortices with the same circulation should contribute negligible kinetic energy away from the vortices.

Even though this isn’t one of the goals of this paper, to assess how well our algorithm with localized transformations can reconstruct the low-energy configurations, we first computed the lowest possible energies for short filaments. We fully enumerated all filaments of lengths up to $N=18$ , and then based on these results we made simplifying assumptions to extend the results up to $N=21$ . These results (with 6 significant digits) are shown in Table 2.

Looking at the limited set of data, it appears to exhibit a strong linear pattern. The best fit line has the equation $E_{\text{min}}~\left(-0.604671N+1.13186\right)/\left(4\pi \right)$ and $R^2=0.996184$ . We then used our Monte Carlo algorithm and attempted to locate the minimum-energy configurations for larger values of N. Since these have not been validated by another approach, we will not report them here. These computed energies follow the linear trend and when using them to generate a best-fit line, we obtain the equation

$E_{\text{min}}~\frac{1}{4\pi }\left(-0.612418N+1.19011\right)$ (18)

with $R^2=0.998585$ . Using this slightly steeper line to generate predictions for the minimum energies for $N=\mathrm{100,}\cdots \mathrm{,1000}$ , we obtain the results in the column labeled “Predicted $E_{\text{min}}$ ” in Table 3. The next two columns of the table show the lowest energies encountered during the computations of the average energies for each N shown in Figure 5, together with their percentage proportions with respect to the predicted minima. Notice how these values track the predictions in the second column of the table.

To illustrate the idea of low-energy configurations being folded and compressed, we show in Figure 13 the lowest-energy configuration encountered during the simulation with $N=1000$ . Notice how this configuration is “compressed” into a relatively small volume, suggesting that the true energy-minimizing configurations might also be volume minimizing in some sense. This is also the case for smaller values of N (not shown here). As the figure suggests, the shown configuration is likely not a minimum as there are clearly visible regions where more “folding” could occur, likely bringing the energy further down. As finding minimum configurations wasn’t the goal of this project, we didn’t pursue it any further.

It is also interesting to point out that all the lowest-energy configurations, whether exact or found approximately in our simulations, ended up with their endpoints in the same unit square on the cubic lattice as illustrated in Figure 13. This observation has not been reported in any of the related works [4] [5] [6] [7] [8] and it has several consequences. First, any of the quantities computed in [4] [5] [6] [7] [8] based on the distance between the endpoints of the filaments (denoted by ${\mu }_{\mathrm{1,}N}$ , ${\mu }_{\mathrm{2,}N}$ , $\mu$ , $D_1$ , $\tilde{D}$ , $\gamma$ , etc. in the references) need to be computed in a different way, and therefore we will not attempt to reproduce those results here. Second, it is likely that this observation is a general attribute of energy-minimizing configurations and worth proving analytically. We have not succeeded in this effort. Third, the numerical search for energy-minimizing configurations may perhaps be done more efficiently using this observation. It limits the set of possible filament configurations that need to be considered, but it also offers a different point of view: consider such filaments as closed (also known as self-avoiding polygons) and then remove one or two segments, depending on the parity of N. We used the first idea to slightly expand the list of values in Table 2 to obtain the results for $19\le N\le 21$ , which were then matched by the results of

our numerical algorithm. Although not presented here, some preliminary work has also been done on considering these configurations as subwalks for self-avoiding polygons.

8. Conclusions

In this paper we explore vortices using a statistical mechanical approach. We consider vortices in negative, infinite, and positive temperature states, a concept previously studied in [15] [16] . This novel approach gives explanations of the role played by suction or supercritical vortices in tornadogenesis by extending the investigations of A. Chorin and his collaborators [4] [5] [6] [7] [8] in understanding turbulence. Supercritical vortices are straight, narrow, and high-energy vortices that are transient in the tornadic flow, and they are approximated in the model by self-avoiding vortex filaments (or self-avoiding walks, SAWs) on a cubic lattice. We were able to reliably compute equilibrium energies of such vortices for a wide range of temperatures (or inverse temperatures), significantly extending the range in the previous works. Having approximated the energy values, we also proposed and utilized an efficient way to compute equilibrium entropy for all temperature values.

Our results confirm that the supercritical vortices (modeled by straight filaments) correspond to the highest kinetic energy of the flow and also correspond to negative temperatures in this model. The lowest-energy configurations are folded up and compressed to a great extent. From the point of view of the flow of energy, the negative-temperature, high-energy supercritical vortices are expected to lose (some of) their energy to the surrounding flow, which can now be understood as a natural process of increasing the entropy of the system. In [8] Chorin states: “The Euler and the Navier-Stokes equations cause filaments to stretch and fold, and N increases. The energy is an increasing function of both T and N; if energy is conserved and N increases T must decrease.” Consequently, in the context of supercritical vortices in a tornadic flow, when such high-energy vortices stretch, they need to fold. As can be seen in our results in Figure 5: with a fixed energy $\langle E\rangle$ and at negative temperatures, longer vortices (larger N) are farther away from their maximum energies than shorter vortices, and thus farther away from straight configurations, which corresponds to more folds present along such vortices. This is in agreement with the results shown in Figure 7.

Our results are not intended to simulate the dynamics of supercritical vortices, only their corresponding statistical equilibrium energies. From this point of view, the results show that a more folded vortex corresponds to a flow with lower kinetic energy. Hence, if a vortex of a fixed length folds, the flow would lose some of its kinetic energy, and this loss would have to be interpreted as a transfer of the energy to the surrounding flow.

We also see that since entropy is highest for $\beta =0$ ( $T=\pm \infty$ ), a system containing a supercritical vortex, which corresponds to high negative temperatures, would be naturally driven towards the state where $\beta =0$ and thus lower the energy associated with the vortex. This corresponds to the vortex folding and energy being transferred to the surrounding flow.

In this work we proposed and utilized a simple algorithm for the computation of statistical equilibrium quantities on a cubic lattice when both an energy and a statistical temperature are involved. This was a consequence of the fact that the pivot algorithm used in [5] [6] [7] works well only for a small range of temperatures near the polymeric case, but it fails in other situations. Specifically, non-straight configurations at large (in absolute value) negative values of β (for which the Boltzmann probability distribution (4) strongly favors high-energy, straight configurations) are extremely unlikely to straighten out in the pivot algorithm MCMC simulation since intermediate steps with non-local transformations are required that greatly affect the energy. Similarly, when starting with a straight configuration (that can be viewed as corresponding to a large, in absolute value, negative β), the same need for an intermediate transformation with a large energy change affects the results at smaller (in absolute value) values of β.

At the other end of the inverse temperature spectrum, when the values of β are large positive, the Boltzmann probability distribution strongly favors lowest-energy configurations, very much folded up and compressed into a small volume. While the maximum-energy configurations are obviously the straight configurations with their energies easily computed for any N, the exact minimum-energy configurations are not known to us, nor are the actual minimum energies. However, the pivot algorithm MCMC simulation struggles to get anywhere close to the hypothesized minimum values, thus begging for a different approach to the problem. Some possibilities might include a statistical approach studied, for example, in [55] , or a quantum variational approach proposed in [56] .

The algorithm utilized in this work seems to perform well for all possible temperature values. It uses two types of transformations: one designed specifically to help with straightening out non-straight configurations in order to increase their energy at negative values of β, and one designed to help with compressing configurations in order to lower their energy at positive values of β. The latter was motivated by previous work for generating self-avoiding walks with fixed endpoints [40] . The algorithm includes many transformations used in previous works, such as the end flips, corner flips, crankshaft moves, and some of the pull moves. Transformations used in the literature to find minimum-energy configurations, such as the pull moves [41] or the re-bridging moves [42] , are not specifically included because the goal of our work is not energy minimization, but their addition is worth exploring in a future work. A new discovery, not reported in any earlier work, is that the lowest-energy configurations have the initial and terminal points on the same unit square of the cubic lattice, hence resembling a self-avoiding polygon (SAP). The meaning of this discovery is worth further exploration.

Having reliably approximated the values of equilibrium energy, we also proposed an efficient way to compute equilibrium entropy for all temperature values using an algorithm that mimics the midpoint Euler’s method for numerically approximating solutions of differential equations in that it approximates the solution to $\partial S/\partial \langle E\rangle =k_B\beta$ instead of computing the entropy directly. The algorithm requires the knowledge of $\langle E\rangle$ for a range of values of β, as well as an initial value of the entropy S for one value of β. This value (for $\beta =0$ ) was readily estimated using existing works on approximate number of possible N-step self-avoiding walks on a cubic lattice. As far as we are aware, we have not seen this approach used in any previous work.

In this work we have not addressed some of the other concepts approximated in the earlier works of Chorin et al. [4] [5] [6] [7] [8] , such as the computation of the fractal dimensions, the Flory exponents, etc. Since these concepts were based on end-to-end distance in the vortex and since the endpoints tend to end up close to each other at lower temperatures, this approach does not seem reasonable in this case. We hope to revisit this subject in a future work. The idea of using the Wang-Landau algorithm to approximate the density of states [49] might also be attempted in this case. It would be interesting to see how well this approach works in this setting, because the energy spectrum is much larger than in previous works, in which it might have a small, finite number of values (up to 83 values in [45] ). Finally, one might also consider the addition of the “pull moves” [41] or other moves designed to help with “tight,” folded configurations in order to see if their addition would result in any difference on the low-energy side of the temperature spectrum.

Acknowledgements

Sincere thanks to the Office of Undergraduate Research and Graduate Opportunity (URGO) at Augsburg University and to Dean and Amy Sundquist for providing the research opportunity and financial support for Eric Bibelnieks, Robert Laskowski, and Alek Lukanen. The authors would also like to thank Arkady Shemyakin and the anonymous referees for their useful feedback and suggestions.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References