Precise determination of lattice phase shifts and mixing angles

We introduce a general and accurate method for determining lattice phase shifts and mixing angles, which is applicable to arbitrary, non-cubic lattices. Our method combines angular momentum projection, spherical wall boundaries and an adjustable auxiliary potential. This allows us to construct radial lattice wave functions and to determine phase shifts at arbitrary energies. For coupled partial waves, we use a complex-valued auxiliary potential that breaks time-reversal invariance. We benchmark our method using a system of two spin-1/2 particles interacting through a finite-range potential with a strong tensor component. We are able to extract phase shifts and mixing angles for all angular momenta and energies, with precision greater than that of extant methods. We discuss a wide range of applications from nuclear lattice simulations to optical lattice experiments.


Introduction
Lattice methods are widely used in studies of quantum few-and many-body problems in nuclear, hadronic, and condensed matter systems, see e.g. Refs. [1][2][3][4][5]. A necessary step in such studies is the computation of scattering phase shifts and mixing angles from an underlying microscopic lattice Hamiltonian. Remarkably, the same problem arises in the context of experiments on optical lattices. Several groups have pioneered the use of ultracold atoms in optical lattices produced by standing laser waves, to emulate the properties of condensed matter systems and quantum field theories [6][7][8][9][10]. The basic concept is to tune the interactions of the atoms, both with each other and with the optical lattice, to reproduce the single-particle properties and particleparticle interactions of the "target theory". Such studies often require a more general setup than a simple cubic lattice, for instance in the case of the hexagonal Hubbard model [11], which closely resembles the physics of graphene [12] and carbon nanotubes [13]. Clearly, a robust and accurate method for computing scattering parameters on arbitrary lattices is needed.
For the scattering of particles on a cubic lattice, Lüscher's finite-volume method [14] uses periodic boundary conditions to infer elastic scattering phase shifts from energy eigenvalues. The method has been widely used in lattice QCD simulations with applications to different angular momenta [15][16][17][18][19] as well as partial-wave mixing [20], see Ref. [2] for a recent review. An important advantage of Lüscher's method is that periodic boundary conditions are typically already used in lattice calculations of nuclear, hadronic, ultracold atomic, and condensed matter systems. Since no additional boundary constraints are needed, the method is easily applied to a wide class of systems.
However, Lüscher's method requires that the finitevolume energy levels can be accurately determined, with errors small compared to the separation between adjacent energy levels. This is not practical in cases such as nucleus-nucleus scattering, where the separation between finite-volume energy levels is many orders of magnitude smaller than the total energy of the system. Fortunately, this problem has been solved using an alternative approach called the adiabatic projection method [21][22][23][24]. There, initial cluster states are evolved using Euclidean time projection and used to calculate an effective two-cluster Hamiltonian (or transfer matrix). In the limit of large projection time, the spectral properties of the effective two-cluster Hamiltonian coincide with those of the original underlying theory. This method has been applied to nuclei and ultracold atoms, while applications to lattice QCD simulations of relativistic hadronic systems are currently being investigated.
Since the adiabatic projection method reduces all scattering systems to an effective two-cluster lattice Hamiltonian, additional boundary conditions can be applied to the effective lattice Hamiltonian in order to compute scattering properties. This opens the door to methods more accurate than Lüscher's by removing the effects of the periodic boundary conditions, which are otherwise a significant source of rotational symmetry breaking. One promising approach is to place the particles in a harmonic oscillator potential and extract phase shifts from the energy eigenvalues [25,26]. Another prominent example is the method used in Refs. [5,27], whereby a "spherical wall" is imposed on the relative separation between the two scattering particles. Phase shifts are then determined using the constraint that the wave function vanish at the wall boundary. This method has been applied to the two-nucleon problem in lattice effective field theory (EFT) [28][29][30][31] and to lattice simulations of nucleus-nucleus scattering using the adiabatic projection method [21][22][23][24].
In spite of such progress in lattice scattering theory, all methods are still lacking in precision, especially when partial-wave mixing and high angular momenta are concerned. In previous work, numerical approximations were used for the study of coupled-channel systems [5]. We now describe an extension of the spherical wall method, which enables an efficient and precise determination of two-particle scattering parameters for arbitrary energies and angular momenta. We use angular momentum projection and solve the lattice radial equation with spherical wall boundaries, supplemented by an "auxiliary potential". We test our method on a lattice model with strong tensor interactions that induce appreciable partial-wave mixing. We expect our method to be applicable in theoretical lattice studies of nuclear, hadronic, ultracold atomic, and condensed matter systems, as well as in the experimental design of optical lattices. While we discuss only non-relativistic wave mechanics in our examples here, the extension to relativistic systems simply entails replacing the non-relativistic dispersion relation with the relativistic one.

Benchmark system
We begin with the eigenvalue equation where r is the relative displacement, σ i , with i = 1, 2, are the spins of the two scattering nucleons with m N ≡ 2µ = 938.92 MeV. Following Ref. [5], we take with C = −2.00 MeV and R 0 = 2.00 × 10 −2 MeV −1 .
We only consider states of total intrinsic spin S = 1.
The radial equation is where L is the orbital angular momentum and J the total angular momentum. The "effective" potential is for uncoupled channels, and for coupled ones. In the continuum, phase shifts and mixing angles are obtained by solving Eq. (3) using the potentials (4) and (5) with appropriate boundary conditions.
As rotational symmetry is broken by the lattice, the energy eigenstates of Eq. (1) belong to the irreducible representations (irreps) A 1 , A 2 , E, T 1 and T 2 of the cubic group SO(3, Z) rather than the full SO(3) rotational group [5,32,33]. For cubic periodic boundary conditions, as in Lüscher's method [14], the cubic symmetry remains exact, thus our solutions can still be classified by cubic irreps. Nevertheless, the rotational symmetry breaking due to the boundaries makes it difficult to identify states of high angular momentum and to extract scattering parameters. In order to remove these effects, we impose a hard spherical wall of radius R W , where θ is the Heaviside step function and Λ is a (large) positive constant, intended to sufficiently suppress the wave function beyond R W (we set Λ = 10 8 MeV). We take R W to exceed the range of the interaction, such that the boundary is placed in the asymptotic (noninteracting) region. We also take 2R W to be less than the difference of the box size and the interaction range, which ensures that cubic boundary effects remain negligible.

Angular momentum decomposition
Let | r ⊗ |S z denote a two-body quantum state with separation r and z-component of total intrinsic spin S z . We define radial lattice coordinates (ρ, ϕ) by grouping equidistant mesh points, as shown in Fig. 1. To construct radial wave functions, we project onto states with total angular momentum (J, J z ) in the continuum limit, using where the Y L,Lz are spherical harmonics with orbital angular momentum (L, L z ). The C J,Jz L,Lz,S,Sz are Clebsch-Gordan coefficients. The parentheses around J, J z and L on the left hand side signify that these quantum numbers are not exactly good quantum numbers. Note that Eq. (7) is applicable to arbitrary geometries. Here, n runs over all lattice points and the "radial shell" is given by the integer m. Then, ρ m is the distance from the origin in units of the lattice spacing a, and δ ρm,| n| picks out all lattice points for which ρ m = | n|. It may be  7). The colors show the magnitude of the matrix elements. To study unphysical mixings, we remove the tensor component of V J (r). The resulting Hamiltonian matrix should ideally be block-diagonal in the S-, Dand G-waves etc. Clearly, the matrix elements that cause unphysical mixings are suppressed by several orders of magnitude. In each block, the row and column indices represent the radial coordinates of the mesh points. For higher partial waves, entire "radial shells" ρ m vanish due to the angular dependence of the wave function, and such redundant rows and columns have been removed. practical (especially for non-cubic lattices) to relax this condition to include all lattice points with |ρ m −| n|| < δ for small, positive δ. On the lattice, the |m J,Jz L form a complete (but non-orthonormal) basis. We therefore compute the norm matrix of these states before solving for the eigenstates of the lattice Hamiltonian.
We find that rotational symmetry breaking is almost entirely due to the non-zero lattice spacing a. As we take a → 0 at fixed R W , rotational symmetry is exactly restored. The degree of mixing between different total angular momenta J and J is a useful indicator of rotational symmetry breaking. Such effects can be interpreted as arising from the non-orthogonality of wave functions in different partial waves when their inner product is computed as a sum over discrete lattice points. The degree of mixing is difficult to estimate a priori, as it depends strongly on the details of the interaction.
Given a simple cubic lattice with a cubic-invariant interaction, unphysical J-mixing only occurs between cubic irreps of the same type. If the objective is to describe a rotationally invariant system on the lattice, then we may simply drop all unphysical couplings between channels with different J. We find that rotational For this purpose, we use a simple cubic lattice with a = 100 MeV −1 and R W = 10.02a. In the radial basis (7), the Hamiltonian matrix becomes nearly blockdiagonal, with each block corresponding to a specific J.
The non-block-diagonal elements induce unphysical Jmixing. In Table 1, we examine the lowest energy levels with and without J-mixing matrix elements. When Jmixing is included, we solve directly for the eigenstates of the lattice Hamiltonian without a spherical harmonic projection. In Fig. 2, we show the Hamiltonian matrix elements in the projected basis defined in Eq. (7). In order to focus entirely on unphysical mixings caused by rotational symmetry breaking, we have neglected the tensor component of V J (r) in Fig. 2. The magnitude of such unphysical mixing matrix elements is found to be greatly suppressed.

Auxiliary potential
We first consider uncoupled channels, where V vanishes beyond an "inner" radius R I . A hard wall at R W gives access to discrete energy eigenvalues only, and a very large box is needed at low energies. To resolve these issues, we define an "outer" radius R O , between R I and R W , as shown in Fig. 1. We also introduce a Gaussian "auxiliary" potential in region III, with R O ≤ r ≤ R W , where the separation between R O and R W is chosen such that V aux is negligible at R O . Note that V aux vanishes in regions I and II. The energy eigenvalues can now be adjusted continuously as a function of V 0 . In Fig. 1, we show V J (r) for V 0 = −25 MeV.
In order to extract phase shifts, we express ψ(r) in region II as for R I ≤ r ≤ R O , where h + J (kr) and h − J (kr) are spherical Bessel functions, and k = √ 2µE. The constants A and B can be determined e.g. by a least-squares fit in region II. We note that with S ≡ exp(2iδ J ), from which δ J can be obtained. For coupled channels, ψ has two components with L = J ± 1. Given Eq. (5), both satisfy the spherical Bessel equation in region II, and are therefore of the form (9). If we denote A ≡ (A J−1 , A J+1 ) T and B ≡ (B J−1 , B J+1 ) T , the S-matrix couples channels with L = J ± 1. In the Stapp parameterization [34], where J is the mixing angle. When solving S from Eq. (10) as in the uncoupled case, we encounter a subtle problem. For a simple hard wall boundary, only one independent solution per lattice energy eigenvalue is obtained. In order to determine S unambiguously, two linearly independent vectors A and B are needed. In Ref. [5], this problem was circumvented by taking two eigenfunctions with approximately the same energy and neglecting their energy difference. However, such a procedure introduces significant uncertainties.
As the potential (5) is real and Hermitian, an exact time-reversal symmetry results. We now add to V J (r) an imaginary component, where U aux (r) is an arbitrary, real-valued function with support in region III only. This leaves V J (r) Hermitian and the energy eigenvalues real, while the time-reversal symmetry is broken. Also, ψ and ψ * are now linearly independent and satisfy Eq. (3) in regions I and II with identical energy eigenvalues. In addition to Eq. (10), we have the conjugate expression, and the S-matrix from (10) and (13). Phase shifts and mixing angles can then be obtained from Eq. (11). Note that the inverse in Eq. (14) cannot be computed without U aux (r), since in that case A = −B * . We use for R O ≤ r 0 ≤ R W , where r 0 is a radial mesh point in region III and U 0 is an arbitrary real constant. We find that the distortion of the energy eigenvalues and radial wave function introduced by this choice is minimal. The same methodology we have applied here for coupled partial waves can also be applied to more general problems with different scattering constituents.

Numerical results
We benchmark our method numerically with the interaction (2) using a cubic lattice with a = 100 MeV −1 (π/a = 314 MeV), box size 35a, and we take R I = 9.02a, R O = 12.02a, and R W = 15.02a. For all channels, we use the real auxiliary potential (8), while for coupled channels we add the complex auxiliary potential (15) with U 0 = 20.0 MeV and r 0 R W .
In Fig. 3, we show our lattice phase shifts and mixing angles. We compare with continuum results, obtained by solving the Lippmann-Schwinger equation for each channel. All our lattice results agree well with the continuum ones, from threshold to a relative center-ofmass momentum of p CM ≡ k = 140 MeV. We note the marked improvement over Ref. [5] for the same benchmark system.

Application to arbitrary lattices
While Lüscher's method has been extended to asymmetric rectangular boxes [35], no standard method yet exists for an arbitrary lattice. Our method can be used to characterize particle-particle interactions on arbitrary lattices, in any number of spatial dimensions. This is significant for optical lattices, as the lattice geometry is then engineered to reproduce the single-particle energies of a given condensed matter or quantum field theoretical system. Anisotropic lattices exhibit more breaking of rotational invariance than a simple cubic lattice does. This is often an essential feature, e.g. in the crossover from a three-dimensional system to a layered two-dimensional one. In Fig. 4, we show the 1 S 0 phase shift on an anisotropic rectangular lattice, where the spacing along the z axis, a z , exceeds those along the x and y axes, denoted collectively by a. The unit cell volume is 100 3 MeV −3 in all cases. While we find good agreement with the continuum up to a z 1.4a, this breaks down when a z becomes comparable to the range of the interaction, with increasing deviation at high p CM . Such a crossover to two-dimensional behavior can be characterized in terms of mixing between the 1 S 0 and 1 D 2 (J z = 0) partial waves, an effect of rotational symmetry breaking. The low-energy particleparticle interactions of any lattice system can be similarly described.

Summary and discussion
We have described a general and systematic method for the calculation of scattering parameters on arbitrary lattices, which we have benchmarked using a lattice model of a finite-range interaction with a strong tensor component. Extensions to more general interactions are straightforward. The Coulomb interaction can be accounted for by replacing the spherical Bessel functions by Coulomb functions, and by defining the distance between particles as the minimum distance on a periodic lattice. The spherical wall then removes unphysical boundary effects. When combined with the adiabatic projection method, the techniques we have discussed can be applied to any scattering system in nuclear, hadronic, ultracold atomic, or condensed matter physics. We expect our method to be applicable to optical lattice experiments, in addition to its immediate usefulness for lattice studies in nuclear, hadronic, and condensed matter theory. In fact, the method proposed here has already been used to significantly improve the adibatic projection method, as detailed in Ref. [36].