Dynamics of Ion Channels via Non-Hermitian Quantum Mechanics

Abstract

We study dynamics and thermodynamics of ion transport in narrow, water-filled channels, considered as effective 1D Coulomb systems. The long range nature of the inter-ion interactions comes about due to the dielectric constants mismatch between the water and the surrounding medium, confining the electric filed to stay mostly within the water-filled channel. Statistical mechanics of such Coulomb systems is dominated by entropic effects which may be accurately accounted for by mapping onto an effective quantum mechanics. In presence of multivalent ions the corresponding quantum mechanics appears to be non-Hermitian. In this review we discuss a framework for semiclassical calculations for the effective non-Hermitian Hamiltonians. Non-Hermiticity elevates WKB action integrals from the real line to closed cycles on a complex Riemann surfaces where direct calculations are not attainable. We circumvent this issue by applying tools from algebraic topology, such as the Picard-Fuchs equation. We discuss how its solutions relate to the thermodynamics and correlation functions of multivalent solutions within narrow, water-filled channels.

1. Introduction

Transport of ions through narrow channels plays a big role in many biological and technological systems. Many pathogens attack cells by forming nanopores in the cell membrane by using pore-forming toxins (PFTs) [1,2]. This punches holes in the cell membrane through which ions diffuse to the outside, effectively killing the cell. Physically this is similar to artificial nanopores in, e.g., silicon [3,4]. These are heavily used in genetic sequencing techniques where high-throughput of selective transport is the most important factor [5]. Other similar examples include free-standing silicon nanowires [6,7] and water-filled nanotubes [8,9]. These systems play various different roles in biology and technology. However they all follow the same underlying physics of a quasi-1D statistical system formed by ions confined to move in a narrow water-filled tube inside a lipid membrane or solid medium [10,11,12,13,14,15,16]. What makes this system special is the large ratio between the dielectric constants of water, ${\kappa }_1\simeq 80$, and the surrounding media (e.g., for lipids or silicon oxide ${\kappa }_2\simeq 2-4$). Because of this, the electric field created by an ion within is confined to stay mostly inside the water-filled channel and does not leak into the surrounding medium. As several numerical simulations in three dimensions point out the flow field also follows almost entirely the channel direction [17,18,19]. This simple observation has profound consequences.

First, as was noticed by Parsegian [20], there is a potential barrier for an ion to enter the channel. This barrier is equal to the energy difference between an ion being inside and outside the channel. For a channel of radius a the electric field created by an ion of charge e in the middle of the channel is $E_0=2e/\left({\kappa }_1{a}^2\right)$. The corresponding field energy integrated over the channel volume is $U_0=\frac{{\kappa }_1}{8\pi }{E}_0^2\pi a^{2}L=k_{B}T\left({\lambda }_{B}L\right)/\left(2a^2\right)$, where L is the length of the channel and ${\lambda }_B=e^2/\left({\kappa }_1{k}_{B}T\right)\approx 7$Å is the Bjerrum length at ambient temperature [11]. For a typical channel with $L\approx 40$Å and $a\approx 5$Å the corresponding (self-)energy barrier exceeds ambient temperature $k_{B}T$ by a factor of 5 or 6. This means that such a channel would block the transport of ions. However, there are at least two mechanisms which can be employed to overcome this issue. One is placing charged radicals along the channel path. The other is entropic screening of the barrier by a collective effect of multiple cations and anions inside the channel. In this review we focus on this latter phenomena, while the former is addressed in References [13,14,21,22].

The second consequence of the mismatch of dielectric constants is that the mutual interactions between the ions within the channel acquire the form of the 1D Coulomb potential

$$\mathsf{\Phi }(x_i-x_j)=e{E}_0|x_i-x_j|,$$

(1)

where $x_i$ are 1D coordinates of the ions along the channel axis. As illustrated in Figure 1, the electric field lines emanating from a charge are bent to run along the channel. Only after a characteristic length $\xi$ given by the implicit relation ${\xi }^2=a^2{\kappa }_1/\left(2{\kappa }_2\right)ln(2\xi /a)$ the field lines start penetrating the lipid membrane and escaping the channel [11]. For a water-filled channel in a lipid membrane this gives $\xi \approx 7a$. Hence, for a sufficiently short channel with $L<\xi$ or (as considered in Section 3) a large concentration of salt ions where the characteristic distance between two ions is smaller than $\xi$, the interactions effectively follow the 1D Coulomb potential. The linear nature of the potential (1) leads to the curious observation that the energy barrier of transporting a charge through the channel can’t be less than $U_0$, irrespective of how many other ions are present in the channel [20]. Indeed, for the most favorable arrangement of alternating positive and negative ions, the electric field along the channel alternates between $\pm E_0$. This leads back to the value of $U_0$ for the electrostatic energy of adding a single ion to the channel in the presence of the other ions. This may seem as a predicament that collective screening can’t lower the transport barrier. Such conclusion is premature, however. The resolution of this apparent paradox is that in a system of multiple particles at a finite temperature it is the free energy (rather than the energy) which determines the transport barrier. The difference between the two is given by the entropy, i.e., it is the entropy of the ion gas within the channel which provides the screening mechanism. The nature of entropic suppression of the transport barrier can be traced to the aforementioned independence of the energy $U_0$ of the positions of individual ions. This observation implies that there is a large number of microscopic configurations which are close in energy. This is the hallmark of a state with large entropy and thus lower free energy.

Formalizing these observations is not entirely straightforward. As was first realized by Edwards and Lenard (EL) in 1962 [10] it requires mapping of the 1D statistical system onto an effective quantum mechanics with cosine potential. In fact, this is a particular case of the generic correspondence between D-dimensional statistical mechanics of the Coulomb gas and (D-1)-dimensional sine-Gordon field theory [23]. The D=2 version of this mapping is well-known in the physics of the Berezinskii-Kosterlitz-Thouless transition. The less appreciated fact is that the Hermitian potential of the form $2cos\theta =e^{i\theta }+e^{-i\theta }$ is a consequence of having a neutral plasma of monovalent ions with charge $\pm e$. In the EL mapping the $e^{\pm i\theta }$ operators shift the value of the electric field in the channel (the variable canonically conjugated to $\theta$) by a quanta $\pm 2E_0$, which corresponds to the electric field generated by a unit charge $\pm e$.

What happens in the presence of a multivalent dissociated salt, such as, e.g., $CaC{l}_2$ which produces a plasma with positive charges $+2e$ and twice as many negative charges $-e$? It is not difficult to see that the EL mapping leads to an effective Hamiltonian with the potential $\frac{1}{2}{e}^{2i\theta }+e^{-i\theta }$. Such a Hamiltonian is non-Hermitian and thus admits a complex-valued spectrum. This may present a problem for the interpretation of the original statistical mechanics of the Coulomb plasma. For example, the free energy density (a manifestly real quantity) is given by the logarithm of the partition function which therefore needs to be real and positive. Fortunately the effective non-Hermitian quantum operator obeys the so-called $\mathcal{PT}$-symmetry [24], which ensures that all eigenvalues are real or appear as complex-conjugate pairs. When calculating the partition function, which includes summing over all eigenvalues, the imaginary parts cancel and we obtain a real, physical result [25]. However, in general there exist complex eigenvalues (spontaneously broken $\mathcal{PT}$-symmetry). This translates to an oscillatory character of certain correlation functions, reflecting short-range charge density wave correlations within the channel.

To model the transport of ions through the channel in this framework we use the concept of boundary charges which was developed in Reference [11]. From now on we assume that the channel is sufficiently short so that all field lines stay inside the channel. If there are no ions inside the channel (or the sum of all charges is zero), then there is no electric field emanating from the channel. If a single ion is added in the center of the channel, then half of its electric field lines are exiting the channel on the left and the other half on the right, cf. Figure 1. This is akin to having two image boundary charges $q,q^{\prime }=\frac{1}{2}$ at the two ends of the channel (charges are measures in units of e). These charges are provided by polarization effects in the well-conducting reservoirs. There are only integer charges inside the channel. Hence, if the boundary charge at one end is q (the ion emits a fraction q of its field lines at one end), then the other boundary charge is $q^{\prime }=1-q$. Reference [11] shows that moving a unit probe charge through the channel (while allowing the other ions to equilibrate) creates boundary charges which change from zero to one. Once the boundary charges reach an integer value they may either be released from the end points and join the bulk, or enter into the channel. This makes thermodynamic properties periodic functions of q with unit period. In Section 2 we show that the boundary charge q takes the role of the quasi-momentum in the effective quantum mechanics. Hence, the bandwidth of the lowest quantum-mechanical band translates directly to the transport barrier.

This review is devoted to the mathematical apparatus needed to treat the non-Hermitian operators appearing in the physics of multivalent 1D plasmas. However, we want to stress that these methods can be applied more broadly to a wide range of non-Hermitian systems. In particular we focus on semiclassical methods applicable for relatively large concentrations of the dissociated salts. Our central observation is that the corresponding (complex) semiclassical trajectories may be viewed as closed cycles on Riemann surfaces of non-zero genus. The action integrals along such cycles are given by solutions of the Picard-Fuchs differential equation, allowing for their analytic evaluation. As a result one obtains asymptotically exact thermodynamic and correlation functions of the 1D multivalent Coulomb plasmas. Of particular interest is the transport barrier, given by the width of the lowest Bloch band (i.e., energy difference between anti-periodic and periodic ground-states of the Schrödinger equation). We obtain analytic results for the transport barriers for various combinations of ion valencies as functions of salt concentration and temperature.

The structure of this paper is as follows: in Section 2 we discuss the EL mapping of statistical mechanics to an effective quantum mechanics with a cosine potential and its generalizations to the non-Hermitian cases. Section 3 is devoted to the semiclassical treatment of the corresponding non-Hermitian operators using the Picard-Fuchs equation. At the end of that section we go beyond the usual semiclassical formulas and describe how to obtain second- and higher-order corrections with little computational effort. We provide a brief summary and discussions in Section 4.

2. Thermodynamic Description and Equivalent Quantum Mechanics

In this section we discuss the relationship between statistical mechanics of the ion channel and (non-Hermitian) quantum mechanics. We start with a thermodynamic description of the ion channel in terms of the grand-canonical partition function. Then we review how to map the partition function onto a Feynman propagator and derive a Hamilton operator from there. This mapping was pioneered by Edwards and Lenard [10] and subsequently used in several works as starting point [11,25,26,27]. If the system consists of cations and anions with the same valency and concentration, then the resulting Hamilton operator is Hermitian. However, if the positive and negative charges have different valency, for example solutions of the divalent salts $MgC{l}_2$ or $CaC{l}_2$, non-Hermitian terms appear. Hence, the spectrum of the resulting operator also contains complex eigenvalues. We discuss how reality and positivity of the partition function is ensured. In the end we comment on the case if charge neutrality is violated.

2.1. Derivation of the Hamilton Operator

As discussed in Section 1 charged ions inside the channel interact with the effective 1-dimensional Coulomb potential $\mathsf{\Phi }\left(x\right)=-e{E}_0\left|x\right|$, where $E_0=2e/{\kappa }_1{a}^2$ is the electric field strength generated by a single ion with charge e inside a channel of radius a and dielectric constant ${\kappa }_1$ [11]. The total interaction energy of all ions in the channel is given by

$$U=\frac{1}{2}{\int \int }_0^{L}dxd{x}^{\prime }\rho \left(x\right)\mathsf{\Phi }(x-x^{\prime })\rho \left(x^{\prime }\right).$$

(2)

Here we write the charge density for point charges in terms of $\delta$-functions,

$$\rho \left(x\right)=\sum _{j=1}^{N_1+N_2}{\sigma }_j\delta (x-x_j)+q\left(\delta \left(x\right)-\delta (x-L)\right),$$

(3)

where ${\sigma }_j=n_1$ for $1\le j\le N_1$ and ${\sigma }_j=-n_2$ for $N_1+1\le j\le N_1+N_2$. This charge density represents $N_1$ cations with valency $n_1$ and $N_2$ anions with valency $-n_2$, and the two fractional boundary charges $\pm q$ at $x=0,L$. The channel is open and can exchange particles with two 3D bulk reservoirs at the ends. Therefore the thermodynamic properties are given by the grandcanonical partition function,

$$\mathcal{Z}=\sum _{N_1,N_2=0}^{\infty }\frac{f_1^{N_1}{f}_2^{N_2}}{N_1!N_2!}\prod _{j=1}^{N_1+N_2}{\int }_0^{L}d{x}_j{e}^{-U/k_{B}T},$$

(4)

where $f_{1,2}$ are the fugacities of the two charge species. As shown in References [10,11] and in Appendix A, the partition function can be converted into a functional integral by introducing an auxiliary field $\theta \left(x\right)$ as conjugate to the charge density $\rho \left(x\right)$. Through this process all integrals over the variables $x_j$ decouple, bringing them to the form ${\sum }_N{[f\int dx{e}^{i\sigma \theta \left(x\right)}]}^N/N!=exp\{f\int dx{e}^{i\sigma \theta \left(x\right)}\}$. The interaction potential (2), being inverse of the 1D Laplace operator, leads to an additional term $exp\{(k_{B}T/e{E}_0)\int dx\theta {\partial }_x^2\theta \}$. As a result the partition function (4) is identically written in terms of the Feynman path integral with an “imaginary time” x, describing quantum mechanics with the Hamiltonian

$$\widehat{H}={(i{\partial }_{\theta }-q)}^2-\left({\alpha }_1{e}^{i{n}_1\theta }+{\alpha }_2{e}^{-i{n}_2\theta }\right),$$

(5)

where ${\alpha }_{1,2}=f_{1,2}{k}_{B}T/e{E}_0$ are dimensionless ion concentrations. The Feynman integral is the expectation value of the evolution operator over the imaginary “time” L,

$${\mathcal{Z}}_L=〈q\left|\mathcal{X}{e}^{-\frac{e{E}_0}{k_{B}T}{\int }_0^{L}dx\widehat{H}}\right|q〉=\sum _m{\left|\langle q|m\rangle \right|}^2{e}^{-\frac{e{E}_{0}L}{k_{B}T}{\epsilon }_m\left(q\right)},$$

(6)

where $\mathcal{X}$ is the x-ordering operator. Here ${\left\{{\epsilon }_m\left(q\right)\right\}}_m$ is the spectrum of the effective Hamiltonian $\widehat{H}$, and $|m\rangle ={\psi }_m\left(\theta \right)$ are its eigenvectors in the Hilbert space of periodic functions, ${\psi }_m\left(\theta \right)={\psi }_m(\theta +2\pi )$. The matrix elements are $\langle q|m\rangle ={\int }_0^{2\pi }d\theta e^{iq\theta }{\psi }_m\left(\theta \right)$. The boundary charge q plays the role of the Bloch quasi-momentum and the spectrum is periodic in q with unit period.

Note that for ${\alpha }_1={\alpha }_2$ and $n_1=n_2$ the potential in Equation (5) reduces to the cosine function and the Hamiltonian becomes the well-known Mathieu Hamiltonian [10]. However, if these conditions are violated the potential is non-Hermitian [25]. We discuss implications of this in the following section.

2.2. Physical Observables

The partition function in Equation (6) gives the thermodynamic properties of the ion gas. However, to be physically meaningful the partition function needs to be real and positive, while the spectrum of the non-Hermitian Hamiltonian (5) may contain non-real eigenvalues. This issue is resolved because the Hamiltonian obeys a symmetry akin to $\mathcal{PT}$-symmetry. The combined action of the “parity operator” $\mathcal{P}:\theta \to -\theta$ and “time reversal” $\mathcal{T}:i\to -i$ leaves the Hamiltonian in Equation (5) unchanged. Bender et al. [24] proved that all eigenvalues of a $\mathcal{PT}$-symmetric Hamiltonian are either real or appear in complex conjugated pairs. Hence, summing over all eigenvalues in Equation (6) gives a real result. In [25] is was shown that for positive values of concentrations ${\alpha }_{1,2}>0$ the lowest energy band ${\epsilon }_0\left(q\right)$ is entirely real-valued, ensuring positivity of the partition function. The higher bands ${\epsilon }_m\left(q\right)$ are in general complex-valued.

Hence we obtain a physically meaningful partition function, and can connect it to thermodynamic observables. The pressure of the Coulomb gas is its free energy per unit length

$$P=k_{B}T\frac{\partial ln{\mathcal{Z}}_L}{\partial L}\stackrel{L\to \infty }{⟶}-e{E}_0{\epsilon }_0\left(q\right),$$

(7)

which for a long channel is determined by the eigenvalue with the smallest real part, ${\epsilon }_0\left(q\right)$. In equilibrium the system minimizes its free energy by choosing an appropriate boundary charge q. In [25,26] this minimum was found to generally be the non-polarized state of the channel, i.e., $q=0$. Adiabatic charge transfer through the channel is associated with the boundary charge q sweeping through its full period. As a result, the (free) energy barrier for ion transport is

$$U_0=e{E}_{0}L{(\Delta \epsilon )}_0,$$

(8)

where ${(\Delta \epsilon )}_0$ is the width of the lowest Bloch band. Therefore the ground state energy and the width of the lowest Bloch band of the Hamiltonian (5) give the leading thermodynamic and transport properties of the $(n_1,n_2)$ Coulomb gas. In Section 3 we discuss analytic results for the eigenvalues and the bandwidth.

2.3. Charge Non-Neutrality

In [10] it was shown that for arbitrary values of ${\alpha }_{1,2}$ the Hamiltonian (5) is always isospectral to a similar charge-neutral Hamiltonian. This can be seen by shifting the coordinate as $\theta \to \theta +{\theta }_0$. Upon such transformation the dimensionless concentrations ${\alpha }_{1,2}$ renormalize as ${\alpha }_1\to {\alpha }_1{e}^{i{n}_1{\theta }_0}$ and ${\alpha }_2\to {\alpha }_2{e}^{-i{n}_2{\theta }_0}$. Notice that the combination ${\alpha }_1^{n_2}{\alpha }_2^{n_1}$ remains invariant. Hence, the family of Hamiltonians (5) with

$${\alpha }_1^{n_2}{\alpha }_2^{n_1}=\mathrm{const}$$

(9)

is isospectral [10,25]. Therefore one may choose one representative from each isospectral family. A convenient choice is taking the representative with charge neutrality in the bulk reservoirs, i.e., $n_1{\alpha }_1=n_2{\alpha }_2\equiv \alpha$. The physical reason for this symmetry is that the interior region of the channel always preserves charge neutrality due to the large self-energy of charges. The edge regions screen charge imbalances of the reservoirs. Therefore, irrespective of the relative fugacities of cations and anions in the reservoirs, the thermodynamics of the long channel are equivalent to the one in contact with neutral reservoirs with an appropriate salt concentration $\alpha$. This brings the Hamiltonian (5) to the form

$$\widehat{H}=\alpha \left[{\widehat{p}}^2-\left(\frac{1}{n_1}{e}^{i{n}_1\theta }+\frac{1}{n_2}{e}^{-i{n}_2\theta }\right)\right],$$

(10)

where we define the momentum operator as

$$\widehat{p}={\alpha }^{-1/2}(-i{\partial }_{\theta }+q);[\theta ,\widehat{p}]=i{\alpha }^{-1/2}.$$

(11)

The commutation relation shows that ${\alpha }^{-1/2}$ plays the role of the effective Planck constant. Hence, a large concentration of charges corresponds to the semiclassical limit of the Hamiltonian (10). We further rescale the eigenvalues $\epsilon$ as

$$u\equiv \frac{n_1{n}_2}{n_1+n_2}\frac{\epsilon }{\alpha }.$$

(12)

This keeps the classical minimum of the potential at $u=-1$, irrespective of the concentration $\alpha$ and the valencies $n_1,n_2$. In Section 3 we discuss the spectral properties of the Hamiltonian (10) in the semiclassical limit.

3. Large Charge Concentration

In Section 2 we mapped the grand-canonical partition function of the Coulomb gas onto an equivalent quantum system. The resulting Hamiltonian, Equation (10), contains one free parameter $\alpha$ which is proportional to the concentration of charged ions. In this section we analyze the spectral problem of this Hamiltonian in the limit of large $\alpha$. As argued after Equation (11), this is the semiclassical limit of the equivalent quantum problem. We use the main semiclassical results, Bohr-Sommerfeld quantization and Gamow’s formula, to calculate the eigenvalues and bandwidths of the Hamiltonian for several different cases of valencies $(n_1,n_2)$. In the case of equal valencies, $n_1=n_2$, the Hamiltonian (10) is the well-known Mathieu Hamiltonian which we discuss in Section 3.1. It’s spectral properties were calculated using several different approaches [10,11,25,26,27,28]. In this review we focus on an approach based on integration on a complex Riemann surface [25,26,28]. We choose this method because it can also be applied to the cases with different valencies, $n_1\ne n_2$, see Section 3.2. In that situation the Hamiltonian is non-Hermitian, and the required action integrals are not attainable by straightforward integration. Instead we show how to relate them to integrals along closed cycles on a Riemann surface. Then we use powerful tools from algebraic topology to derive a differential equation for the action integrals. This is known as the Picard-Fuchs equation. The required actions are a combination of the solutions of this differential equation. Through this procedure we bypass the use of direct integration methods. From the actions we obtain the eigenvalues and the bandwidths, which are directly related to the ion pressure and transport barrier for ions in the channel. In Section 3.3 we go one step further. We use the same concepts to calculate the second-order corrections in the WKB series. Most importantly we show that these can be expressed in terms of the already-calculated action and its derivatives, and therefore can be obtained with minimal computational effort. This gives an improved semiclassical approximation of the eigenvalues. Relating this to the pressure in the ion channel we find that beyond the ideal-gas pressure and the Debye-Hueckel correction there is another correction which only depends on the geometry of the channel but not on the concentration of ions. We compare these results to numerical calculations.

3.1. Equal Valency

As mentioned in Section 2 the Hamiltonian in Equation (5) is Hermitian if the valencies of the two charges are equal, $n_1=n_2$. Indeed, in this case it reduces to the well-known Mathieu Hamiltonian,

$$\mathcal{H}=\alpha \left[{\widehat{p}}^2-2cos\theta \right].$$

(13)

In literature there exist several studies of the Coulomb gas with charges of equal valency. In [10] it was first noted that the Coulomb gas is mapped onto the Mathieu equation. In [27] the authors perform a semiclassical calculation on this equation via direct integration. From this they obtain the required actions and analytic approximations of the eigenvalues and bandwidths. [11] provides additional qualitative arguments which lead to the same results. However, as mentioned above, in this section we will follow the Riemann surface methods developed in [25] because in that framework one can also study the case of unequal valencies $n_1\ne n_2$ in Section 3.2, and these concepts form the basis of our considerations for higher-order corrections in Section 3.3.

3.1.1. Construction of the Riemann Surface

In the semiclassical ansatz we look for wave functions of the form $\psi =e^{i{\alpha }^{1/2}S}$, where S is the action of the classical problem with the normalized Hamiltonian (13). The semiclassical trajectories satisfy the classical Hamilton equations of motion and thus conserve the (complex) energy u in Equation (12),

$$2u=p^2-2cos\theta .$$

(14)

In this normalization $u=\mp 1$ corresponds to the bottom (top) of the cosine potential. Our approach to calculate the action integrals $S={\oint }_{\gamma }p(\theta ,u)d\theta$ is based on complex algebraic topology. First we set $z=e^{i\theta }$ and consider $(z,p)$ as complex variables. Energy conservation, Equation (14), defines a family of complex algebraic curves parametrized by u and satisfying

$${\mathcal{E}}_u:\mathcal{F}(p,z)=p^{2}z-(z^2+2uz+1)=0.$$

(15)

For $u\ne \pm 1$ it can be checked that $(\partial \mathcal{F}/\partial z,\partial \mathcal{F}/\partial p)$ does not vanish on ${\mathcal{E}}_u$, so each ${\mathcal{E}}_u$ is nonsingular. Then $\mathcal{F}(p,z)$ implicitly defines a locally holomorphic map $p=p\left(z\right)$. The exceptions to this occur at $z=0,\infty ,z_{\pm }$, where $z_{\pm }=-u\pm i\sqrt{1-u^2}$ are the roots of $p^2=0$ (i.e., classical turning points). In a vicinity of these four branch points $p\left(z\right)$ behaves as

$$\begin{array}{ccc}\hfill & p\sim z^{-1/2},\hfill & \hfill (z\sim 0)\\ \hfill & p\sim z^{1/2},\hfill & \hfill (z\sim \infty )\\ \hfill & p\sim {(z-z_{\pm })}^{1/2},\hfill & \hfill (z\sim z_{\pm })\end{array}$$

(16)

respectively, i.e., $p\left(z\right)$ is locally double-valued. Note that we added the point at infinity to have an even number of branch points. This compactifies the complex plane and makes it topologically equivalent to a Riemann sphere, cf. Figure 2.

To avoid dealing with $p\left(z\right)$ as a double-valued function we introduce a second copy of the complex z-plane and the corresponding Riemann sphere. On both sheets we define two branch cuts connecting the four branch points, between $0,\infty$ and the turning points $z_{\pm }$ respectively. $p\left(z\right)$ is analytically continued across the branch cuts, i.e., when crossing a branch cut we jump from the first sheet to the second and vice versa. Identifying the branch cuts as edges we can deform the two Riemann spheres into tubes and glue them together to form a torus. This construction is visualized in Figure 2. Thus, the complex algebraic curve ${\mathcal{E}}_u$ in Equation (15) defines a torus which is a compact Riemann surface of genus $g=1$. (Generically, every compact Riemann surface is topologically equivalent to a sphere with some number of handles g, or a (multi-)torus with g holes, called the genus of the surface).

3.1.2. Integrals on the Riemann Surface and the Picard-Fuchs Equation

The action integrals can be understood as integrals over closed cycles $\gamma$, $S\left(u\right)={\oint }_{\gamma }\lambda \left(u\right)$, where

$$\lambda \left(u\right)=p\left(\theta \right)d\theta =p\left(z\right)\frac{dz}{iz}=\frac{{(z^2+2uz+1)}^{1/2}}{i{z}^{3/2}}dz$$

(17)

is the action 1-form which, by construction, is holomorphic on the Riemann surface.

To visualize the relevant trajectories we momentarily return to $\theta$ and consider it as complex. In this representation one has square-root branch cuts along the real axis, connecting the classical turning points along the classically allowed region. The integration trajectories run just above or below the real axis and connect the turning points. After combining them to form closed cycles one can push these cycles off the real axis and away from the turning points without altering the integrals (by Cauchy’s theorem). We call these the classical cycle ${\gamma }_0$ and the instanton cycle ${\gamma }_1$, as shown in Figure 3. Translating these two cycles to the complex z-plane yields the contours in the right panel of that figure.

Cauchy’s theorem is also valid on the Riemann surface since the action form (17) is, by construction, holomorphic on the torus. Therefore all closed cycles can be deformed without changing the integrals, and can be expressed as a combination of an integer number of these two basis cycles. This leads to our key idea how to calculate the action integrals: for this we employ a central theorem of algebraic topology, de Rham’s theorem. It states that on a Riemann surface there are exactly as many linearly independent holomorphic 1-forms to integrate upon as there are independent closed cycles to integrate along. This is valid up to exact forms, i.e., 1-forms which integrate to 0 along any closed cycle, and boundaries, i.e., closed curves which can be continuously deformed to a point. Hence, there are exactly two independent holomorphic 1-forms on the Riemann surface. Any set of three 1-forms is linearly dependent modulo an exact form which integrates to 0 upon integration along any closed cycle. (A full explanation of the mathematical concepts is beyond the scope of this review. A detailed discussion of relevant and related concepts is in [29], basic definitions and additional background are in [30,31]. All concepts can also be found online at [32]. A simplified derivation specifically for complex-valued Riemann surfaces is in chapter 2 of [28]).

Equipped with this we look at a set which contains the action 1-form (17) and its first two derivatives with respect to energy u, $\{\lambda \left(u\right),{\lambda }^{\prime }\left(u\right),{\lambda }^{″}\left(u\right)\}$. Taking derivatives does not change the structure of branch points, therefore these are three 1-forms which are all defined on the same Riemann surface. Hence, we know that there must exist a linear combination of these which is an exact form. Reference [25] explains in detail how to find the linear combination and the exact form as

$$\left((u^2-1){\partial }_u^2+\frac{1}{4}\right)\lambda \left(u\right)=\frac{d}{dz}\left[\frac{i}{2}\frac{1-z^2}{z^{1/2}{(z^2+2uz+1)}^{1/2}}\right]dz.$$

(18)

It is evident from Stokes’ theorem that the right-hand-side integrates to 0 along any closed cycle on the Riemann surface. Hence, we obtain

$${\oint }_{\gamma }\left((u^2-1){\partial }_u^2+\frac{1}{4}\right)\lambda \left(u\right)=(u^2-1)S^{″}\left(u\right)+\frac{1}{4}S\left(u\right)=0.$$

(19)

This differential equation for the action $S\left(u\right)$ is called the Picard-Fuchs Equation [29]. Integration is performed along a closed cycle $\gamma$, which can be the classical or the instanton cycle, ${\gamma }_{0,1}$ in Figure 3. Therefore both the classical and instanton actions $S_{0,1}\left(u\right)$ are solutions of the Picard-Fuchs Equation (19). This equation is a second-order ordinary differential equation, therefore it admits two independent solutions. These can be found in the form $F_0\left(u^2\right)$ and $u{F}_1\left(u^2\right)$, where

$$\begin{array}{cc}\hfill F_0\left(u^2\right)& ={}_2{F}_1\left(-\frac{1}{4},-\frac{1}{4};\frac{1}{2};u^2\right),\hfill \\ \hfill F_1\left(u^2\right)& ={}_2{F}_1\left(+\frac{1}{4},+\frac{1}{4};\frac{3}{2};u^2\right),\hfill \end{array}$$

(20)

are hypergeometric functions [33,34]. These solutions form a basis out of which $S_{0,1}\left(u\right)$ must be composed, so we write

$$\begin{array}{cc}\hfill S_0\left(u\right)& =C_{00}{F}_0\left(u^2\right)+C_{01}u{F}_1\left(u^2\right),\hfill \\ \hfill S_1\left(u\right)& =C_{10}{F}_0\left(u^2\right)+C_{11}u{F}_1\left(u^2\right).\hfill \end{array}$$

(21)

To find the correct coefficients $C_{jk}$, $j,k=0,1$ it is sufficient to evaluate the periods at one specific value of u. Employing the fact that the hypergeometric functions (20) are normalized and analytic at $u=0$, i.e., $F_k\left(u^2\right)=1+\mathcal{O}\left(u^2\right)$, one notices that $S_j\left(u\right)=C_{j0}+u{C}_{j1}+\mathcal{O}\left(u^2\right)$. Thus, to identify $C_{jk}$ we expand the integrand $\lambda \left(u\right)$ to first order in u and evaluate the integrals $S_j\left(u\right)$ at $u=0$. Straightforward calculation yields

$$\begin{array}{cc}\hfill & C_{00}=e^{-i\pi /2}{C}_{10}=8{\pi }^{-1/2}\mathsf{\Gamma }{(3/4)}^2,\hfill \\ \hfill & C_{01}=e^{+i\pi /2}{C}_{11}={\pi }^{-1/2}\mathsf{\Gamma }{(1/4)}^2.\hfill \end{array}$$

(22)

The relations between $C_{0k}$ and $C_{1k}$ are not accidental. They originate from the fact that the cycle ${\gamma }_1$ transforms into ${\gamma }_0$ by substitution $z^{\prime }=e^{-i\pi }z$ and $u^{\prime }=e^{i\pi }u$, and vice versa. This gives a global symmetry between the two periods,

$$\begin{array}{c}\hfill S_0\left(u\right)=e^{-i\pi /2}{S}_1\left(e^{i\pi }u\right).\end{array}$$

(23)

Equations (20)–(23) fully determine the classical and instanton actions $S_{0,1}\left(u\right)$. We now proceed to relate them to physical observables.

3.1.3. Semiclassical Results

We seek semiclassical results for the sequence of low-energy bands terminating at $u=-1$. Therefore we quantize the classical action $S_0\left(u\right)$ according to the Bohr-Sommerfeld rule to determine the normalized energies $u_m$ as solutions of the equation

$$S_0\left(u_m\right)=2\pi {\alpha }^{-1/2}(m+1/2),m=0,1,\dots$$

(24)

We see that the cycle ${\gamma }_0$ contracts to a point when the energy goes to the bottom of the potential, $u\to -1$. This corresponds to vanishing of the classical action, $S_0(u=-1)=0$. To obtain an approximate analytic expression for the lowest energy levels ${\epsilon }_m=2\alpha u_m$ we expand the classical action to first order near the bottom of the potential,

$$S_0\left(u\right)=2\pi (u+1).$$

(25)

Equations (24) and (25) combined imply ${\epsilon }_m=-2\alpha +2{\alpha }^{1/2}\left(m+1/2\right)$. As a result the pressure (7) of a monovalent gas is

$$P=-e{E}_0{\epsilon }_0=2k_{B}Tf-\sqrt{k_{B}Te{E}_{0}f}.$$

(26)

The two terms here are the pressure of the ideal gas with fugacity f and the mean-field Debye-Hueckel interaction correction [22].

The instanton action $S_1\left(u\right)$ determines the bandwidth ${(\Delta u)}_m$ according to Gamow’s formula,

$${(\Delta u)}_m=\frac{\omega }{\pi \sqrt{\alpha }}{e}^{i{\alpha }^{1/2}{S}_1\left(u_m\right)/2}.$$

(27)

Here $\omega =2$ is the frequency of the harmonic-oscillator approximation of the potential near the classical minimum. We expand the instanton action near the classical minimum and at the quantized energies $u_m=-1+{\alpha }^{-1/2}(m+1/2)$ to obtain

$$S_1\left(u_m\right)=16i+2i\left(m+\frac{1}{2}\right)ln\left(\frac{m+1/2}{32e{\alpha }^{1/2}}\right).$$

(28)

Applying this to Gamow’s Formula (27) leads to

$${(\Delta \epsilon )}_m=2\alpha {(\Delta u)}_m=\frac{4}{\pi }{\left(\frac{32e}{m+1/2}\right)}^{m+1/2}{e}^{-8{\alpha }^{1/2}+(m/2+3/4)ln\alpha },$$

(29)

This coincides with the known asymptotic results for the Mathieu Equation [27,35,36]. As explained below Equation (3), adiabatic charge transport is associated with a change of the boundary charge q (i.e., quasi-momentum) across the interval $0<q<1$ (i.e., the Brillouin zone). Therefore the free energy transport barrier is given by the width of the lowest Bloch band, ${(\Delta \epsilon )}_0$. One notices that increasing the concentration of salt ions leads to an exponential entropic suppression of the transport barrier, ${(\Delta \epsilon )}_0\propto {\alpha }^{3/4}{e}^{-8\sqrt{\alpha }}$.

3.2. Multivalent Ions

So far we worked with the Hermitian example of the Mathieu Hamiltonian, i.e., when both ion species are monovalent, $n_1=n_2=1$. With that we could validate the Riemann surface method by comparing the results to literature. In this section we discuss four different cases with multivalent ions (assuming $n_1>n_2$ without loss of generality). In such a scenario the Hamiltonian (5) is non-Hermitian. This leads to complex values in the spectrum, which we present in Section 3.2.1. Furthermore, in classical motion the coordinate and momentum acquire complex values. This results in a phase space $(\theta ,p)$ with two complex dimensions (instead of two real dimensions). The classical (instanton) action is obtained by integrating the momentum $p\left(\theta \right)$ along the trajectory which connects two turning points and solves the classical equations of motion with real (imaginary) time. However, solving the equations of motion in complex phase space $(\theta ,p)$ is non-trivial, if at all attainable. Therefore we go from an integral along the trajectory to an integral along a closed cycle in the plane of complex $z=e^{i\theta }$ which encloses the trajectory, similar to the mapping in Figure 3. With that we connect the non-Hermitian problem to the method that we validated in the previous section. We discuss this calculation for four different combinations of charge valencies in Section 3.2.2. In Section 3.2.3 we connect the results to the classical and instanton actions and physical observables.

3.2.1. Spectrum of the Non-Hermitian Hamiltonian

Non-Hermiticity of the Hamiltonian (10) has a significant effect on its spectrum. Namely, not all eigenvalues are real. In Figure 4 we show numerical results for the eigenvalues at large concentration $\alpha$, for four different combinations of the integers $(n_1,n_2)$. Most importantly all non-real eigenvalues appear as complex conjugate pairs. This is a consequence of the $\mathcal{PT}$-symmetry of the Hamiltonian and crucial to obtain a physically meaningful partition function, as discussed in Section 2. Furthermore we see sequences of narrow bands which emerge from $u=-\nu$ with ${\nu }^{n_1+n_2}=1$. These sequences approximately follow the lines connecting $u=-\nu$ and $u=1$, but avoid the special point $u=1$. At some point all of these branches merge. Beyond this the nature of the spectrum changes drastically, instead of narrow bands and large gaps we see wide bands separated by small gaps. This feature is similar to the case of a periodic Hermitian potential: as long as the energy lies below the maximum of the potential there are narrow bands, while for energies exceeding the maximum there are wide bands. Hence, we associate the point where the spectral branches meet with the top of the potential. (It is important to bear in mind that for a complex-valued potential there is no proper definition of a “maximum”.) The energy variable u is normalized so that in the Hermitian $(1,1)$ case this point lies at $u=1$. In the non-Hermitian cases we observe $u\approx 0.96$ for $(2,1)$, $u\approx 1.09$ for $(3,1)$, $u\approx 1.20$ for $(4,1)$, and $u\approx 0.84$ for $(3,2)$. These values are independent of $\alpha$, so this must be a consequence of the underlying classical mechanics.

To calculate the statistical partition function in Equation (6) the most important eigenvalues are those with small real part. Therefore we will focus on the narrow bands and treat them in semiclassical approximation.

3.2.2. Riemann Surface and Picard-Fuchs Equation

We use the rescaled energy variable u in Equation (12), substitute $z=e^{i\theta }$ in the Hamiltonian (10), and write the classical energy-momentum relation as

$$u\frac{n_1+n_2}{n_1{n}_2}=p^2-\left(\frac{1}{n_1}{e}^{i{n}_1\theta }+\frac{1}{n_2}{e}^{-i{n}_2\theta }\right).$$

(30)

The generalization for the complex algebraic curve in Equation (15) is the family of curves

$${\mathcal{E}}_u:\mathcal{F}(p,z)=n_1{n}_2{p}^2{z}^{n_2}-\left(n_2{z}^{n_1+n_2}+(n_1+n_2)u{z}^{n_2}+n_1\right)=0.$$

(31)

This defines implicitly a double-valued function $p\left(z\right)$. It is easy to see that $(\partial \mathcal{F}/\partial z,\partial \mathcal{F}/\partial p)$ does not vanish on ${\mathcal{E}}_u$ unless $u=-e^{\frac{2\pi im}{n_1+n_2}}$ for an integer m. For the non-singular values of u the function $p\left(z\right)$ is locally holomorphic except for the points $z=0,\infty ,z_j$, where $z_j,j=1,\dots ,n_1+n_2$ are the roots of $p^2=0$. The $z_j$ are the turning points of classical motion in complex coordinates. Near these special points $p\left(z\right)$ behaves as

$$\begin{array}{ccc}\hfill & p\sim z^{-n_2/2},\hfill & \hfill (z\sim 0)\\ \hfill & p\sim z^{n_1/2},\hfill & \hfill (z\sim \infty )\\ \hfill & p\sim {(z-z_j)}^{1/2}.\hfill & \hfill (z\sim z_j)\end{array}$$

(32)

The $z_j$ are $n_1+n_2$ branch points. If $n_2$ ($n_1$) is odd, then 0 (∞) is an additional branch point; for even $n_2$ ($n_1$) there is a normal pole at 0 (∞). Hence, there are $n_1+n_2+1$ branch points on the Riemann sphere if one of the integers is odd, and $n_1+n_2+2$ branch points if both are odd. (Here we ignore the case that $n_1,n_2$ are both even, because if both integers can be divided by the same number n we can define $z^{\prime }=e^{in\theta }$ to obtain a simpler algebraic curve.) In all cases there is an even number of branch points which can be connected pairwise to form branch cuts. For $(n_1,n_2)=(2,1)$ we obtain four branch points and two branch cuts and a Riemann surface of genus 1, as in Figure 2. For $(n_1,n_2)=(3,1),(4,1),(3,2)$ the asymptotic expansions (32) give six branch points. Consequently there are three branch cuts in the complex plane. Through a similar construction as in Figure 2 one obtains a Riemann surface which is topologically equivalent to a figure “8”, i.e., a figure with two holes and genus 2 [26,28]. In the following we consider these four cases because there are no naturally occurring ions with larger charge. However, mathematically the algebraic curves for higher values of the integers can be constructed in the same way, yielding Riemann surfaces with larger genus.

In Figure 5 we show the structure of branch points in the z-plane for these four cases. On a Riemann surface with genus $g=1\left(2\right)$ there are two (four) independent closed cycles [29]. In Figure 5 we define three cycles for the $(2,1)$ case, and five cycles for $(4,1)$ and $(3,2)$. This is done for convenience and symmetry reasons. The superfluous cycle can be expressed by the other cycles. For $(2,1)$ the linear combination ${\gamma }_0-{\gamma }_1-{\gamma }_2$ does not contain any of the branch points and is contractible to a point. For $(4,1)$ the trivial cycle is ${\gamma }_0-{\gamma }_1+{\gamma }_2+{\gamma }_3-{\gamma }_4\cong 0$, and for $(3,2)$ we see that ${\gamma }_0+{\gamma }_1-{\gamma }_2-{\gamma }_3+{\gamma }_4\cong 0$. We choose to include the additional cycle because it gives an easy representation for the symmetry relation between the corresponding actions $S_j\left(u\right)$, akin to Equation (23). By substituting $z^{\prime }=e^{-i\varphi }z$ and $u^{\prime }=e^{i\varphi }u$ the cycles transform ${\gamma }_j\to {\gamma }_{j+1}$. For the $(2,1)$ case the resulting symmetry relation is

$$S_0\left(u\right)=e^{\pi i/3}{S}_1\left(e^{-2\pi i/3}u\right)=e^{-\pi i/3}{S}_2\left(e^{2\pi i/3}u\right).$$

(33)

The analogous symmetry relations for the genus-2 cases are shown in Reference [26].

To calculate the actions $S\left(u\right)={\oint }_{\gamma }\lambda \left(u\right)$ we continue in the same manner as in Section 3.1. The 1-form (cf. Equation (17)) with general $n_1,n_2$ is

$$\lambda \left(u\right)=p\left(\theta \right)d\theta =p\left(z\right)\frac{dz}{iz}=\frac{{\left(n_2{z}^{n_1+n_2}+(n_1+n_2)u{z}^{n_2}+n_1\right)}^{1/2}}{i\sqrt{n_1{n}_2}{z}^{1+n_2/2}}dz.$$

(34)

On a Riemann surface of genus $g=1\left(2\right)$ there are two (four) independent closed cycles. According to the de Rham theorem, this is equal to the number of linearly independent 1-forms, modulo exact forms. Therefore a set of the 1-form (34) and its first few derivatives, ${\left\{{\partial }_u^k\lambda \left(u\right)\right\}}_{k=0}^K$, is linearly dependent if it contains the first $K=2\left(4\right)$ derivatives. We build a linear combination of these which equals an exact form (for details see [26]). The integral of the exact form along a closed cycle gives zero. What is left is a linear combination of the action and its first derivatives, cf. Equation (19). In the $(2,1)$ case we find this Picard-Fuchs equation as

$$(u^3+1)S_j^{″}\left(u\right)+\frac{u}{4}{S}_j\left(u\right)=0.$$

(35)

This is a second-order differential equation. The Picard-Fuchs equations for the genus-2 cases are fourth-order ODEs which can be found in Reference [26]. Equation (35) admits two solutions $F_0\left(u^3\right)$ and $u{F}_1\left(u^3\right)$ which are given in terms of hypergeometric functions [26,34],

$$\begin{array}{ccc}\hfill F_0\left(u^3\right)& =& {}_2{F}_1\left(-\frac{1}{6},-\frac{1}{6};\frac{2}{3};-u^3\right),\hfill \\ \hfill F_1\left(u^3\right)& =& {}_2{F}_1\left(+\frac{1}{6},+\frac{1}{6};\frac{4}{3};-u^3\right).\hfill \end{array}$$

(36)

The actions are a linear combination of these, $S_j\left(u\right)=C_{j0}{F}_0\left(u^3\right)+C_{j1}u{F}_1\left(u^3\right)$. Expanding the hypergeometric functions near the origin, $F_{0,1}\left(u^3\right)=1+\mathcal{O}\left(u^3\right)$, one notices that $S_j\left(u\right)=C_{j0}+u{C}_{j1}+\mathcal{O}\left(u^3\right)$ as $u\to 0$. The constants $C_{0k}$ are therefore given by $C_{00}=S_0\left(0\right)$ and $C_{01}=S_0^{\prime }\left(0\right)$. Straightforward integration and the symmetry relation (33) yield

$$\begin{array}{ccc}\hfill C_{00}=C_{10}{e}^{\pi i/3}=C_{20}{e}^{-\pi i/3}& =& \frac{2^{11/6}3{\pi }^{3/2}}{\mathsf{\Gamma }\left(\frac{1}{6}\right)\mathsf{\Gamma }\left(\frac{1}{3}\right)},\hfill \\ \hfill C_{01}=C_{11}{e}^{-\pi i/3}=C_{21}{e}^{\pi i/3}& =& \frac{3^{1/2}\mathsf{\Gamma }\left(\frac{1}{6}\right)\mathsf{\Gamma }\left(\frac{1}{3}\right)}{2^{11/6}{\pi }^{1/2}}.\hfill \end{array}$$

(37)

The actions $S_j\left(u\right)$ for $(n_1,n_2)=(2,1)$ are fully given by Equations (33), (36), and (37). The analogous expressions for the genus-2 cases with $(n_1,n_2)=(3,1),(4,1),(3,2)$ are given in Reference [26]. In the next section we discuss how to obtain semiclassical results for the physical observables.

3.2.3. Semiclassical Results in the Non-Hermitian Cases

In this section we calculate the eigenenergies and bandwidths of the non-Hermitian Hamiltonian in Equation (10) with the Bohr-Sommerfeld quantization condition and Gamow’s formula. To utilize these standard semiclassical results we need to calculate the classical and the instanton actions, $S_{cl,inst}\left(u\right)={\oint }_{{\gamma }_{cl,inst}}\lambda \left(u\right)$. The crucial part hereby is identifying the correct cycle of integration. In Section 3.1, when discussing the case of a Hermitian Hamiltonian, we identified these with trajectories which connect the classical turning points through the classically allowed or forbidden region respectively, cf. Figure 3. In the non-Hermitian case this is not so clear, because there exist more than two turning points, and in the space with complex coordinate, momentum, and energy the concept of classically allowed or forbidden regions doesn’t apply. Instead, to identify the correct actions $S_{cl,inst}\left(u\right)$ we look at the analytic behavior of these actions near special values of the energy u.

The Bohr-Sommerfeld condition requires that the classical action goes to zero at the classical minimum of the potential. This happens when two turning points collide which causes the corresponding cycle of integration to collapse to a point. We can easily check that in all four cases in Figure 5 the cycle ${\gamma }_0$ collapses to a point as $u\to -1$. The corresponding action goes to zero, $S_0(-1)=0$. Therefore we identify $S_0\left(u\right)$ as the classical action which quantizes into the branch of eigenstates that terminates at $u=-1$. For $(n_1,n_2)=(2,1)$ it follows immediately from the symmetry relation (33) that at the singular point $u=e^{i\pi /3}$ ($e^{-i\pi /3}$) the cycle ${\gamma }_1$ (${\gamma }_2$) collapses to a point and the action $S_1\left(u\right)$ ($S_2\left(u\right)$) goes to zero. It should be thus identified with the classical action for the spectral branch terminating at $u=e^{i\pi /3}$ ($e^{-i\pi /3}$). In the same manner the analogous symmetry relations for the genus-2 cases in Reference [26] allow us to identify the classical actions for all the spectral branches in Figure 4. Quantizing these classical actions according to the Bohr-Sommerfeld rule,

$$S_j\left(u_m^{\left(j\right)}\right)=2\pi {\alpha }^{-1/2}(m+1/2),m=0,1,\dots ,$$

(38)

one finds the semiclassical energies $u_m^{\left(j\right)}$ determining the $q=0$ edges of the narrow bands in the complex plane. These results are compared with numerical data in Figure 6. The excellent agreement holds all the way up to the point where all spectral branches coalesce. Beyond this point the semiclassical approximation breaks down, which manifests in e.g., the appearance of wide Bloch bands.

All graphs exhibit spectral branches along the lines where one of the actions $S_j\left(u\right)$ is real, while the narrow bands lie at the points determined by the Bohr-Sommerfeld condition (38). For $(2,1)$ and $(3,1)$ there exists a total of three spectral sequences, for $(4,1)$ and $(3,2)$ five sequences due to a higher number of special energies. In the $(4,1)$ case the two complex-valued branches intersect at $u\approx 0.90+0.32i$. Beyond this point the two sequences merge into one, for which the quantization condition is neither determined by $S_1$ nor $S_2$ individually, but instead by the sum $S_1+S_2$ (shown in green). For $(3,2)$ the two lines for the complex-conjugate pair $S_2$ and $S_3$ collide at $u\approx 0.84$, the other pair collides at $u\approx 0.98$ where the semiclassical approximation breaks down. A closer look at the state at $u\approx 0.89$ reveals that this cannot be explained by the quantization of $S_0$ along the real axis. However, it meets the Bohr-Sommerfeld condition (38) for $S_2+S_3$ with $m=17$. Thus we may conclude that the spectral branches can be derived from the Bohr-Sommerfeld condition for one of the actions, or upon intersection of two branches by the sum of the two actions of these branches.

To calculate the width of these bands with Gamow’s formula,

$${(\Delta u)}_m=\frac{\omega }{\pi \sqrt{\alpha }}{e}^{i{\alpha }^{1/2}{S}_{inst}\left(u_m\right)/2},$$

(39)

we need to identify the instanton actions. The classical frequency $\omega$ is determined from the harmonic oscillator approximation, i.e., by expanding the potential around $\theta =0$. In Hermitian quantum mechanics the instanton trajectory connects the two classical turning points through the classically forbidden region, cf. Figure 3. Hence, we identify the instanton cycle as the other possible cycle that connects the same two turning points. This is a combination of all other integration cycles ${\gamma }_i$. The instanton actions that correspond to the classical actions $S_0\left(u\right)$ are

$$\begin{array}{ccc}\hfill S_{inst}\left(u\right)=& -S_1\left(u\right)+S_2\left(u\right),& (2,1);\hfill \\ \hfill S_{inst}\left(u\right)=& -S_1\left(u\right)-S_2\left(u\right)+S_3\left(u\right),& (3,1);\hfill \\ \hfill S_{inst}\left(u\right)=& -S_1\left(u\right)-S_2\left(u\right)+S_3\left(u\right)+S_4\left(u\right),& (4,1);\hfill \\ \hfill S_{inst}\left(u\right)=& -S_1\left(u\right)+S_2\left(u\right)-S_3\left(u\right)+S_4\left(u\right),& (3,2).\hfill \end{array}$$

(40)

From the symmetry relation (33) between the actions and its analogons for the genus-2 cases it is easy to check that these combinations are purely imaginary, which makes the bandwidth in Equation (39) real, as required.

More can be said when considering the analytic structure of the classical and instanton action in a vicinity of $u=-1$. Therefore we use a concept called monodromy [29,32], which is visualized in Figure 7. We choose some $u\gtrsim -1$ and allow u to wind around $-1$ (i.e., $(u+1)\to (u+1)e^{2\pi i}$). The two branch points inside the cycle ${\gamma }_0$ in Figure 5 are exchanged by this transformation via a counter-clockwise half-turn; the branch cut in effect rotates by ${180}^{\circ }$. For ${\gamma }_0$ this has no effect, the cut turns within it. Not so for ${\gamma }_1$: if this cycle is never to intersect the branch points, it is continuously deformed and as a result of this monodromy transformation we obtain ${\gamma }_1\to {\gamma }_1+{\gamma }_0$, thus $S_1$ picks up a contribution of $S_0$. This effect is visualized in Figure 7. While we have returned to the initial value of u, the period $S_1$ does not return to its original value and thus can’t be analytic. This occurs for every monodromy cycle near $u=-1$. The only function which monotonically increases as the phase of its argument grows is the complex logarithm. Thus, $S_1$ must have a logarithmic dependence on $1+u$. One can check that

$$S_1\left(u\right)=Q_1\left(u\right)-\frac{i}{2\pi }{S}_0\left(u\right)ln(1+u)$$

(41)

yields the correct behavior, where $Q_1\left(u\right)$ and $S_0\left(u\right)$ are analytic functions of $(1+u)$. The same applies to the other cycle which is connected to the same branch cut. Therefore the instanton action $S_{inst}$ in Equation (40) picks up a contribution of $-2S_0$. Hence, we can derive the Bohr-Sommerfeld quantization condition (38) from the requirement that the monodromy transformation leaves the bandwidth (39) unchanged.

A comparison of the results for the bandwidth with numerical simulations is shown in Figure 8 for the four non-Hermitian cases and the Hermitian $(1,1)$ case. All cases show good agreement with the numerical data already for moderate values of the parameter $\alpha$. (Note however, that for the genus-2 cases Gamow’s formula had to be multiplied by an overall factor of $3/2$ (in $(3,1)$ case) or 2 (in $(4,1)$ and $(3,2)$ cases), respectively. The origin of this preexponential factor is beyond the scope of this paper.)

To summarize, we find that in all cases the bandwidth is of the form

$${(\Delta \epsilon )}_m=A\times {\left(\frac{k}{m+1/2}\right)}^{(m+1/2)}\times exp\left(-b\sqrt{\alpha }+(m/2+3/4)ln\alpha \right).$$

(42)

The pressure, which is calculated from the lowest eigenvalue, contains the ideal gas pressure and the Debye-Hueckel correction,

$$P=C{k}_{B}Tf-c\sqrt{k_{B}Te{E}_{0}f}.$$

(43)

Here A, k and b, and C and c, are numerical factors that can be calculated directly by expanding $S_0$ and $S_{inst}$:

$({\mathit{n}}_{\mathbf{1}},{\mathit{n}}_{\mathbf{2}})$	$\mathit{A}$	$\mathit{k}$	$\mathit{b}$	$\mathit{C}$	$\mathit{c}$
$(1,1)$	$4/\pi$	$32e$	8	2	1
$(2,1)$	$2\sqrt{6}/\pi$	$36\sqrt{6}e$	$3\sqrt{6}$	$3/2$	$\sqrt{3/2}$
$(3,1)$	$4\sqrt{2}/\pi$	$6.35$	$7.06$	$4/3$	$\sqrt{2}$
$(4,1)$	$5\sqrt{5/2}/2\pi$	$1303.46$	$6.90$	$5/4$	$\sqrt{5/2}$
$(3,2)$	$5\sqrt{5/2}/3\pi$	$6740.06$	$5.65$	$5/6$	$\sqrt{5/2}$

These values quantify the thermodynamic properties of the ion channels for all five different combinations of charged ions which give a Riemann surface of genus 1 or 2. With a maximum valency of 4 these are also the physically relevant cases. Most importantly we show that the Coulomb gas with unequal valency $n_1\ne n_2$ has the same qualitative behavior as the standard gas with ions of equal valency, $n_1=n_2$. In all cases the pressure consists of the ideal gas pressure and the Debye-Hueckel correction, see Equation (43). Crucially for transport through the ion channel, in all cases the bandwidth shows exponential decay with the square-root of the fugacity $\alpha$ and has a universal pre-exponential factor of ${\alpha }^{3/4}$. However the factor b in the exponent shrinks when the valency is increased, meaning that the transport barrier falls off slower with increased charge concentration when transporting ions with larger valency.

3.3. Higher-Order Corrections from Exact Wkb Method

The approximations for the eigenvalues of the non-Hermitian Hamiltonian can be improved further by considering second- and higher-order terms in the WKB series. The inspiration comes from the exact WKB method which was studied extensively in the context of resurgence theory [37,38]. We use this to get a better approximation for the eigenvalues, and with that the pressure of the Coulomb gas, at moderate values of the charge concentration $\alpha \gtrsim 1$. The key is that the $q=0$ band edge, which gives the pressure in equilibrium, is determined by an infinite series in ${\alpha }^{-1}$ (i.e., ${\hslash }^2$ in usual quantum mechanics),

$$\sum _{n=0}^{\infty }\frac{{(-1)}^n}{{\alpha }^n}{\oint }_{{\gamma }_{cl}}{\rho }_{2n}(\theta ,u_m)d\theta =\frac{2\pi (m+1/2)}{\sqrt{\alpha }}.$$

(44)

${\rho }_0(\theta ,u)=p(\theta ,u)$ is the classical momentum, and the other terms can be found through a recursive relation [37]. Equation (44) is sometimes also referred to as the generalized Bohr-Sommerfeld quantization condition. Reference [38] shows a calculation of the exact WKB series at all orders for a class of Hermitian genus-1 cases which include the cosine potential, i.e., the $(1,1)$ case in our notation. Here we follow the ideas in [39] and chapter 5 of [28] which give a general procedure to calculate the terms order-by-order for any potential, and can also be applied to non-Hermitian Hamiltonians.

It is evident that truncation of Equation (44) at the $n=0$ term leads to the usual Bohr-Sommerfeld quantization condition. To improve upon this we include the $n=1$ term. The integrand is given by

$${\rho }_2(\theta ,u)d\theta =\left(\frac{{\partial }_{\theta }^2\left({\rho }_0{(\theta ,u)}^2\right)}{48{\rho }_0{(\theta ,u)}^3}+\frac{5}{24}{\partial }_{\theta }\frac{{\rho }_0^{\prime }(\theta ,u)}{{\rho }_0{(\theta ,u)}^2}\right)d\theta ,$$

(45)

where the prime denotes a derivative with respect to $\theta$ [37]. The second term is an exact form which integrates to zero. We drop this exact form, use the expression (30) for the classical momentum $p={\rho }_0$, and perform the coordinate transformation $z=e^{i\theta }$ to write the second-order 1-form as

$${\tilde{\rho }}_2(z,u)dz=\frac{-n_1{z}^{n_1}-n_2{z}^{-n_2}}{48{\left(u\frac{n_1+n_2}{n_1{n}_2}+\frac{1}{n_1}{z}^{n_1}+\frac{1}{n_2}{z}^{-n_2}\right)}^{3/2}iz}dz.$$

(46)

A comparison with Equation (34) shows that the second-order 1-form ${\tilde{\rho }}_2(z,u)dz$ has the same branch points as the action 1-form $\lambda \left(u\right)$. Therefore it is defined on the same Riemann surface. As discussed in the preceding sections, on the Riemann surfaces of genus $g=1\left(2\right)$ there exist two (four) linearly independent 1-forms, up to an exact form. We take ${\left\{{\partial }_u^k\lambda \left(u\right)\right\}}_{k=0}^K$ as this maximal independent set with $K=1\left(3\right)$. This forms a basis for the space of all 1-forms. Hence, the second-order correction can be written as a linear combination of these basis 1-forms, modulo an exact form. We find this linear combination in the same way as in the derivation of the Picard-Fuchs Equations (19) and (35) and integrate it along the classical cycle ${\gamma }_{cl}$ to get

$${\oint }_{{\gamma }_{cl}}{\tilde{\rho }}_2(z,u)dz=-a\left(S_0^{\prime }\left(u\right)+2u{S}_0^{″}\left(u\right)\right),\begin{array}{cccccc}(n_1,n_2)& (1,1)& (2,1)& (3,1)& (4,1)& (3,1)\\ a& 1/48& 1/18& 3/32& 2/15& 3/10\end{array}.$$

(47)

These expressions fully define the second-order corrections in terms of the classical action and its derivatives with respect to u. These are easily obtained from the previous results, Equations (20)–(22), (36) and (37) (see Reference [26] for the genus-2 cases). Note that in the genus-1 cases the second derivative $S_0^{″}\left(u\right)$ can be replaced with $S_0\left(u\right)$ by using the Picard-Fuchs Equations (19) and (35).

Here we want to stress that calculation of the second-order (and any higher) correction is only as computationally demanding as deriving the Picard-Fuchs equation. It does not require solving the differential equation and matching boundary conditions because the correct classical action was already identified. Therefore this can also be used as a simple method to simply calculate the higher-order WKB terms if the classical action was obtained in a different manner. The improvement in the approximation of the lowest eigenvalue is shown in Figure 9.

With the second-order result we can calculate the eigenvalues u up to order ${\alpha }^{-1}$. Therefore we expand the classical action $S_0\left(u\right)$ for $u\gtrsim -1$ to order ${(u+1)}^2$ and solve for u. Taking the lowest eigenvalue $u_0$ and applying this to the formula for the pressure (7) gives

$$P=c_0{k}_{B}Tf-c_1\sqrt{e{E}_0{k}_{B}Tf}-c_{2}e{E}_0,$$

(48)

with the following constants:

$({\mathit{n}}_{\mathbf{1}},{\mathit{n}}_{\mathbf{2}})$	${\mathit{c}}_{\mathbf{0}}$	${\mathit{c}}_{\mathbf{1}}$	${\mathit{c}}_{\mathbf{2}}$
$(1,1)$	2	1	1/16
$(2,1)$	$3/2$	$\sqrt{3/2}$	1/9
$(3,1)$	$4/3$	$\sqrt{2}$	19/144
$(4,1)$	$5/4$	$\sqrt{5/2}$	1/8
$(3,2)$	$5/6$	$\sqrt{5/2}$	13/36

This gives the ideal gas pressure and the Debye-Hueckel correction from the usual Bohr-Sommerfeld condition. The second-order WKB term gives an additional correction which is independent of the fugacity but only depends on the geometric properties of the channel which are included in the definition of $E_0$.

4. Summary of Semiclassical Results

In this review we discussed analytic calculations of the thermodynamic properties of an ion channel at large charge concentrations, with an extension to moderate concentrations. We started with discussing a standard mapping of a statistical system onto an effective quantum system [10,23]. When performing this mapping there is no guarantee that the resulting effective Hamiltonian is Hermitian and has a purely real spectrum. Physically one needs to obtain a real and positive partition function. This is e.g., guaranteed if the Hamiltonian obeys $\mathcal{PT}$-symmetry and its lowest eigenvalue is purely real.

Translation between the quantum results and thermodynamic observables is straightforward. Most importantly, the pressure (i.e., free energy density) is given by the quantum mechanical ground-state energy. The adiabatic transport barrier is the width of the lowest Bloch band. The complex energies of excited states, c.f. Figure 4 and Figure 6, describe higher-order correlation functions. Their imaginary part is responsible for spatial oscillations, while the real part yields an overall exponential decay. Such decaying oscillatory correlation functions reflect short-range charge density wave ionic order within the channel. As seen in Figure 4 and Figure 6, the onset of complex eigenvalues happens at lower energies for ions with larger valencies, which implies stronger charge density fluctuations. In all cases we observe that an increase of the charge concentration leads to an exponential reduction of the transport barrier, however this decay is slower if the ion valencies are large. This is visualized in Figure 8.

The approximation with the effective 1D Coulomb potential, Equation (1), works best at large ion concentration. Electric field lines leak out of the channel after a characteristic length $\xi$ which is given by ${\xi }^2=a^2{\kappa }_1/\left(2{\kappa }_2\right)ln(2\xi /a)$, where a is the radius of the channel and ${\kappa }_1,{\kappa }_2$ are the dielectric constants of water and the surrounding medium. Therefore the 1D Coulomb potential best approximates the situation where the characteristic distance between the ions is small. This is the case of large charge concentration, which is also the case when then semiclassical approximation is applicable.

Here we discuss a method how to perform semiclassical calculations without the need to solve the classical equations of motion and without direct integration. This is particularly useful in the non-Hermitian cases when the solutions to the equations of motions are hardly attainable. Instead we derive and solve the Picard-Fuchs differential equation, which is a tool from algebraic topology. The power of the Picard-Fuchs equation is that it is a coordinate-free expression, i.e., one does not need to know the classical trajectories. In the last part we extend our calculations to second- and higher-order terms in the WKB series. These provide a clearly improved approximation for the eigenvalues especially at moderate charge concentrations, see Figure 9.

The applicability of the Picard-Fuchs method extends far beyond the case of ion channels. It can be a powerful tool for Hermitian and non-Hermitian systems alike, as it can be applied to generic Hamiltonians. Especially the extension to second- and higher-order terms in the WKB series requires very little computational effort once the classical action has been calculated. Mappings of a generic statistical system onto an effective quantum system can lead to a non-Hermitian Hamiltonian for which semiclassical calculations with direct integration are difficult. We believe that the Picard-Fuchs method can be especially useful in these cases, as it allows us to circumvent the complications associated with direct integration like solving equations of motion with complex coordinates.