Quantum Eigenstates of Curved and Varying Cross-Sectional Waveguides

: A simple one-dimensional differential equation in the centerline coordinate of an arbitrarily curved quantum waveguide with a varying cross section is derived using a combination of differential geometry and perturbation theory. The model can tackle curved quantum waveguides with a cross-sectional shape and dimensions that vary along the axis. The present analysis generalizes previous models that are restricted to either straight waveguides with a varying cross-section or curved waveguides, where the shape and dimensions of the cross section are ﬁxed. We carry out full 2D wave simulations on a number of complex waveguide geometries and demonstrate excellent agreement with the eigenstates and energies obtained using our present 1D model. It is shown that the computational beneﬁt in using the present 1D model to calculate both 2D and 3D wave solutions is signiﬁcant and allows for the fast optimization of complex quantum waveguide design. The derived 1D model renders direct access as to how quantum waveguide eigenstates depend on varying cross-sectional dimensions, the waveguide curvature, and rotation of the cross-sectional frame. In particular, a gauge transformation reveals that the individual effects of curvature, thickness variation, and frame rotation correspond to separate terms in a geometric potential only. Generalization of the present formalism to electromagnetics and acoustics, accounting appropriately for the relevant boundary conditions, is anticipated.


Introduction
With the vast possibilities to manufacture complex topological structures in quantum technology, it becomes increasingly important to address how physical properties change due to shape, size etc. Nanowire technology provides a rich platform to discover novel physical properties related to curvature effects such as flexible electronics [1], battery [2] and nanoelectromechanical sensors and generators [3][4][5][6][7]. A selection of recent nanoarchitecture work [8] includes the influence of a varying thickness of 2D WS 2 nanolayers [9], the role of topology of the thermopower of six-terminal Andreev interferometers [10], the voltage induced by moving vortices as a function of transport and magnetic fields for rolled-up nanostructured nanotubes [11], topological transitions in superconducting structures [12,13], etc. Recently, the present authors [14] derived the first four-band spinless k · p model in curved coordinates using Kane's model. In Ref. [15], a detailed discussion of the combined influence of curvature, torsion, and cross-sectional rotation of a nanostructure for electron eigenstates and symmetry properties in nanowires was presented based on a one-dimensional effective mass equation. Evidently, the use of curved coordinates is effective to recast computationally intensive 3D wave problems of complex nanoarchitectures into lower-dimensional problems; examples include the well-known cases of a torus, a helix, and a Möbius structure [16][17][18][19][20][21][22][23][24].
The use of analytical or numerical methods to determine 3D waveguide eigenstates and eigenfrequencies in acoustics and electromagnetics has a long history and is of high importance for applications. Webster's horn model in acoustics provides a simple, accurate one-dimensional approximation for the estimation of the sound field in waveguides of varying thickness that otherwise requires a full 3D analysis [25,26]. Stevenson [27,28] derived a set of coupled ordinary differential equations for the varying thickness waveguide. In the latter case, even a small number of coupled equations was shown to provide an accurate description of the groundstate, the first excited states and their frequencies.
Both methods [25,27,28], however, assumed a straight centerline (cylinder structures) and therefore cannot be applied to a more general class of nanowire geometries characterized by, simultaneously, a curved axis and a varying cross section.
In this work, we present a new method using differential geometry and perturbation theory to determine eigenstates to the Schrödinger equation of a quantum-mechanical particle confined to a curved, varying-thickness waveguide with small cross-sectional dimensions relative to the waveguide length. The cross-sectional geometry can be arbitrary and may rotate along the waveguide axis. We examine the influence of all geometry parameters and compare with accurate full-wave numerical solutions. For the waveguide structures analyzed, excellent accuracy is demonstrated for the first several eigenstates. The developed model provides insight into the influence of waveguide curvature, cross-sectional thickness variation and frame rotation for eigenstate energies and symmetry properties. In particular, the present method sheds new light on resonance phenomena as a result of coupling between geometry parameters that is not attainable from a standard complex full-wave problem. Besides the insight that can be gained from the 1D equations, there is also significant computational advantages, as discussed in Appendix A.

The General Equations
The goal is to express the Laplace operator in tubular coordinates in a thin neighbourhood of a curve. The normal cross sections of the structure we consider (the neighbourhood) are assumed to have the same shape, but the size and orientation in the normal plane are allowed to be arbitrary, see Figure 1, where the shapes of the normal intersections are disks. r h t p q Figure 1. The physical domain. We consider a neighbourhood of a curve r (in black). The tangent vectors t (in red) are orthogonal to the normal planes (in pink). They intersect in this case the neighbourhood in disks of varying radius h centered at the curve. The vectors p and q, in green and blue respectively, give an orthonormal basis for the normal plane.
Nevertheless, it is advantageous to first consider the case where the orientation of the cross sections relative to a minimal rotating frame is constant.

Constant Orientation
We assume that we have a curve r parametrised by arc length with tangent t = r . We first choose a minimal rotating frame, i.e., an orthonormal frame t, p, q (see Figure 1) along the curve satisfying, Then, the curvature of the curve is κ = κ 2 1 + κ 2 2 and the torsion is τ = . We now consider the following parametrisation of a neighbourhood of r, see Figure 1: where (v, w) belongs to a fixed domain Ω 0 in the vw-plane and L is the length of the curve. If we fix s then we obtain a scaled (by h(s)) and rotated copy of Ω 0 , in Figure 1 it is the unit disk. The scale is allowed to vary from point to point, this is encoded in the function h and the overall scaling is given by the positive number . We will expand the Laplace operator in powers of and as we consider as a small number we will disregard terms in of degree one or higher. The partial derivatives are We let and find the metric tensor by taking the inner products between the partial derivatives, The determinant is and the inverse is We let where r = √ v 2 + w 2 . We furthermore have that Denoting the coordinates (s, v, w) by (u 1 , u 2 , u 3 ) we can now write the Laplace operator as As usual [15], we let χ = √ Fψ and have that such that the Helmholtz equation reads Later, we shall exemplify the present analysis to the Schrödinger equation for a quantum-mechanical particle. For now, we focus on the Helmholtz equation. To zeroth order in , we have Unless h = 0, i.e., the thickness is constant, we cannot separate the variables. Instead we will use perturbation theory. We have The last three terms are considered as perturbations and first order perturbation theory yields where R and R 2 are the expectation values of R and R 2 , respectively.
where ψ 0 is a transversal eigenstate. This procedure is easy in the case of non-degenerate transverse eigenstates such as for the transverse groundstate. However, in case of degenerate transverse eigenstates, standard degenerate perturbation theory must be used. We arrive at We can now separate the variables, χ = χ 0 (v, w)φ(s), and if χ 0 is a solution to the 2D Helmholtz equation, i.e.,

Variable Orientation
We now allow the orientation of the cross section to vary relative to the minimal rotating frame. More precisely, consider a general frame t, n 1 , We then have where ω = ϑ . By following steps similar to the analysis in Refs. [15,19], the resulting 1D equation, accounting for a varying thickness and curvature effects, becomes, where L is the θ-component of the orbital angular moment operator, i.e., and θ is the polar angle. Let us perform a gauge transformation, φ → ϕ, Then Hence, Equation (18) is equivalent to A similar transformation was used in Ref. [15] in the case without a thickness variation. We see that all geometric effects: curvature, thickness variation, and frame rotation appear in the form of a geometric potential only. It is hard to say something general about the significance, except that the last term λ 0 /( h) 2 dominates when → 0 and unless h is constant the eigenstates will concentrate at the maxima of h. In Appendix B, we examine examples where either the first or the last term is most important.

Non Arc Length Parametrisation
So far we have assumed that the curve is parametrized by arc length, but even though it is always possible theoretically, it is more often than not impractical. Thus, we will here give the 1D approximative equation for a general parametrisation. Observe that where prime ( ) now denotes differentation with respect to a general parameter u. Then, Equation (18) becomes (22) and Equation (21) becomes It must be pointed out that the above analysis and reduction of a 3D problem to an effective 1D model does not impose restrictions on the form of the boundary conditions at the waveguide ends. Hence, the boundary conditions in the u coordinate can be arbitrary.

2D Problems
The analysis above obviously carries over to 2D problems. The only difference is that we now have a unique frame and the terms involving ω and L disappears. The domain Ω 0 is now an interval on the v-line.
In Appendix A, we will validate the method by comparing 1D calculations with full 2D calculations.
It can also be relevant for 3D problems where one variable can be separated away exactly, cf. Section 3.

Schrödinger Equation
The Schrödinger equation reads where m and E are the particle's mass and energy, respectively, and ψ is the wavefunction. If −∇ 2 ψ = λψ and L 0 is the unit length then we have

Transmission Studies through a Straight Waveguide with a Varying Thickness
In this section, we use the developed 1D model to compute transmission and reflection coefficients of a particle confined to a varying thickness quantum waveguide with a straight axis defined as the x axis. The waveguide cross-sectional dimension in one direction (z) is assumed to be much larger than the dimensions along the x, y directions. In this case, the 3D problem reduces to a 2D problem in the x − y plane. In the region confined by the waveguide, the Schrödinger equation reads and the wavefunction satisfies ψ = 0 on the surface of the waveguide corresponding to infinite barriers. We will attempt to compute the transmission and reflection coefficients for the case where the particle approaches the central region from the far left, and the waveguide thickness is assumed to vary in a continuous (and piecewise two times differentiable) manner. The waveguide thickness is considered constant to the left and the right of the central region (for simplicity, the same thickness h 0 and a zero curvature are assumed in the left and right regions). In this case the wavefunction part φ and its derivative dφ ds must be continuous everywhere, and we may write φ(s) = exp (iαs) + B exp (−iαs) , to the left of the central region, to the right of the central region, where α = 2mĒ The latter inequality guarantees that wave solutions are propagating in the left and right regions. We can now obtain the transmission T and reflection R coefficients as Evidently, since the wave problem is lossless, we must require In order to determine the coefficients B and C, Equation (15) must be solved in the central region by matching to the wave solutions in the left and right regions. This has been done numerically.
We shall consider the central region to have a trigonometrically varying thickness, i.e., where L m defines the size of the central region and n is an integer. In Figure 2 (upper plot), we show the waveguide thickness and its first-and second-order derivatives as a function of the centerline coordinate u for the case where h 0 = 1 nm, b = 2 Å, L m = 10 nm, and n = 1, = 1.
In Figure 2 (lower plot), the transmission and reflection curves as a function of the energy of the quantum-mechanical particle are shown.   In Figure 3 (lower plot), the transmission and reflection curves as a function of the energy of the quantum-mechanical particle are shown. The mass of the particle is the free-electron mass.
Note that a threshold energy is required for the particle to be transmitted to the right. This comes about as the thickness in the central region is lower than the thickness to the far left and right. Hence, the particle will essentially only propagate to the right if its energy is higher than the energy associated with the smallest thickness h min = h 0 − b in the central region. Otherwise, if its energy is lower than h 2 2m π 2 4 2 h 2 min , transmission can only take place by virtue of tunneling, i.e., with a small probability.
In Figure 4 (upper plot), we show the waveguide thickness and its first-and second-order derivatives as a function of the centerline coordinate u for the case where h 0 = 1 nm, b = 2 Å, L m = 10 nm, and n = 2, = 1.
In Figure 4 (lower plot), the transmission and reflection curves as a function of the energy of the quantum-mechanical particle are shown.
In this case, where h(s) follows a full period of a sine curve, the thickness is smaller than the thickness h 0 to the far left and right, for one half of the central region. Hence, again, a threshold is needed for the quantum-mechanical particle to be transmitted to the far right.  We point to that the 1D model can be used in a completely analogous way to calculate transmission properties of an arbitrarily curved quantum waveguide.

3D Eigenstate Problems
We will first consider 3D curves with constant curvature, e.g., a straight line, a helix, and a circle. Moreover, we will consider quantum waveguides with two thickness variations: a linear (with varying slope) and a trigonometric (with varying amplitude, frequency, and phase).
Note that, for all curves with constant curvature, the 1D equation Equation (18) can be treated on equal footings, since in this case, the only effect of curvature is a shift in the energy E by − mκ 2 2h 2 L 2 0 . We will, in all cases, consider a circular cross section. We choose the unit disk as our 2D parameter space and the first step is to calculate the first eigenvalue λ 0 and the expectation values R , R 2 , L , and L 2 . Using polar coordinates (r, θ) in the unit disk the ground state does not depend on the angle θ, so we have immediately that L = L 2 = 0. As a function of the radius r we have a Bessel function J 0 (j 0,1 r), where j 0,1 ≈ 2.4 is the first zero of J 0 and λ 0 = j 2 0,1 ≈ 5.78. We have and hence RJ 0 (j 0,1 r) = rj 0,1 J 0 (j 0,1 r) = −rj 0,1 J 1 (j 0,1 r) ,

Line with Circular Cross Section
We will consider trigonometric variations of the thickness along the line segment (of length L 0 = 10 nm and average thickness 1 nm). More specifically, we have = 0.5 nm, α = 0.5, and consider the cases ψ 0 = 0, π/2 and β ∈ [π/32, 5π]. The result can be seen in Figure 5. As mentioned before, any curve with constant curvature has essentially the same spectrum. The only difference from the straight line is a shift of the whole spectrum by − mκ 2 2h 2 L 2 0 .

Elliptic Helix with Circular Cross Section
By an elliptic helix, we mean a geodesic on an elliptic cylinder. That is, it can be parametrised as where L 0 is the required length and L = T 0 √ 1 + c 2 a 2 sin 2 t + b 2 cos 2 t dt. Notice that the third coordinate is arc length on the ellipse. We consider the case a, b > 0 and find that the speed and curvature is given by If a > b, then the maximal curvature is L L 0 a b 2 (1+c 2 ) . We let T ∈ [π/32, 5π] and have length L 0 = 100 nm and average thickness 1 nm. The energies of the first five states are shown in Figure 6. We clearly have an effect, but not nearly as strong as when we vary the thickness.

Closed Ellipse with Circular Cross Section
We now consider a closed ellipse parametrised as where L is a scaling such that the length is L 0 , i.e., L = 2π 0 a 2 sin 2 t + b 2 cos 2 t dt .
With length L 0 = 10 nm, constant width 1 nm, letting a = 1, and considering b ∈ [0.25, 1] we obtain the energies of the first five states shown in Figure 7. 2 v cos 2 π 2 w dv dw = 1, Consider a trigonometric variation, Equation (34), of the thickness. As a closed structure is considered, we need to have βL = 2πm. In Table 1, Table 1. The first five eigenvalues for the tennis ball curve with b = a/20, length 25 nm, average width 1 nm ( = 0.5 nm), no turn, and trigonometric variation of width h = 1 + α 2 sin(βs − ψ 0 ), α = 0.5, ψ 0 = 0. The first five eigenvalues in the case of no turn, α = 0.5, ψ 0 = 0, L = 25 nm, and = 0.5 nm are listed. The cases of a quarter turn and a half turn yield the same eigenvalues to four digits precision. This is due to the fact that the term ω 2 ( L 2 − L 2 ) is much smaller than the term

Conclusions
A numerically effective 1D model describing eigenstates and energies of an arbitrarily complex quantum waveguide geometry is presented. A combination of differential geometry and perturbation theory is used to derive a 1D model describing the added effects of waveguide curvature, varying cross-sectional dimensions and frame rotation along the waveguide. A detailed comparison of the present 1D model with full wave simulations is carried out and excellent agreement is demonstrated for a number of complex waveguide geometries, many of which cannot be tackled using earlier 1D models. The only significant assumption of the present 1D model is that the cross-sectional dimensions are significantly smaller than the waveguide length, although experience has shown this restriction is not a severe one. A main new result is that the individual effects of curvature, thickness variation, and frame orientation enter as separate terms in a geometric potential only.
In Figure A3, we plot the difference between the 1D and 2D calculations as a function of for different values of α.
When is very small, the difference is seen to increase instead of decrease. This is due to an increase in the discretisation error from solving the exact 2D equation.
For the 1D-equation, we used B-splines of degree 3 and 64 knot intervals. That leads to a system matrix of size 65 × 65. For the 2D calculations we used standard finite element analysis on a mesh of size 256 × 256 and with basis functions of degree 3. That leads to a system matrix of size more than 10 6 × 10 6 (in 3D we would have a size greater than 2 × 10 9 × 2 × 10 9 . The reason we need higher resolution for the 2D calculation is that we are interested in the differences λ n − λ 0 and using the 1D model we calculate this quantity directly, while in the 2D case we calculate λ n and as the difference is small we need a very high precision. Even if we had used a mesh of size 65 in all directions we would have a system matrix of size more than 8 × 10 4 × 8 × 10 4 in 2D and 4 × 10 7 × 4 × 10 7 in 3D. Numerical analyses of the logarithmic spiral in cases with a linearly or a trigonometrically varying thickness were also carried out. Compared to the exponentially varying thickness case, the error analysis gives approximately the same good result for the trigonometric case, while even smaller errors are found for the linearly varying thickness case. Similar conclusions are found for the cases with straight and circular-arc waveguides. We now consider the case R = 1, µ = 0.45π, with θ ∈ [0, 3π] for a linearly varying thickness h = 1 + α (s − β), with α = ±0.8/L and β = L/2, where L = R tan µ (e 3π cot µ − 1) is the length of the arc. The two cases are depicted in Figure A4.    Table A2. The colours represent the fifth eigenfunction calculated based on the 2D equation.
In Table A2, we have listed the first five eigenvalues. In order to compare the eigenfunctions obtained from the 2D and 1D calculations, respectively, we calculate the probability density along the curve. The determinant of the metric tensor is G(s, v) = 2 h(s) 2 F(s, v) hence for the 2D calculation the probability density along the curve is In Figure A5, we have plotted the probability density for the first and fifth state using the 1D and 2D equations. We have only shown the cases of the first and the fifth eigenfunction, but the good agreement between probability densities obtained from the 1D and 2D equations applies to the other cases as well.

Appendix B. Significance of the Different Terms
We consider the case of a line with length L = 1 and no rotation, that is we have the 1D equation Equation (15) with κ = 0. In order to examine the significance of the different terms we will consider the eigenvalue problem where the potential V is given by (A4) Observe that case (A4) gives the full 1D equation. It is clear that the different geometric potentials lead to pronounced differences in the eigenvalues. We furthermore consider = 0.03 and h(s) = 1 + e −β(s−1/2) with β = 20 and 1000. In Table A3 we have listed the first five eigenvalues And in Figure A6 we have plotted the probability distributions for the first three states together with the full geometric potential V 4 .