Phase Unwrapping and One-Dimensional Sign Problems

Sign problems in path integrals arise when different field configurations contribute with different signs or phases. Phase unwrapping describes a family of signal processing techniques in which phase differences between elements of a time series are integrated to construct non-compact unwrapped phase differences. By combining phase unwrapping with a cumulant expansion, path integrals with sign problems arising from phase fluctuations can be systematically approximated as linear combinations of path integrals without sign problems. This work explores phase unwrapping in zero-plus-one-dimensional complex scalar field theory. Results with improved signal-to-noise ratios for the spectrum of scalar field theory can be obtained from unwrapped phases, but the size of cumulant expansion truncation errors is found to be undesirably sensitive to the parameters of the phase unwrapping algorithm employed. It is argued that this numerical sensitivity arises from discretization artifacts that become large when phases fluctuate close to singularities of a complex logarithm in the definition of the unwrapped phase.


I. INTRODUCTION
If the properties of quantum states with large baryon number could be calculated directly from the Standard Model, open questions could be answered regarding the boundaries of the periodic table, the composition of neutron stars, and the interpretation of low-energy experimental searches for beyond the Standard Model interactions in nuclear targets. Many electroweak effects can be calculated accurately in perturbation theory, but at low energies relevant to nuclear systems the effects of the strong nuclear force described by Quantum Chromodynamics (QCD) can only be accurately computed non-perturbatively. Lattice Quantum Field Theory (LQFT) provides a method for nonperturbatively calculating path integrals in many QFTs in which ultraviolet divergences have been regularized by replacing spacetime with a discrete lattice of points and infrared divergences have been regularized by restricting spacetime to a finite volume. Renormalized QFT observables are obtained from the continuum and infinite volume limits of LQFT results. Monte Carlo (MC) methods can be used to calculate path integrals representing observables in Lattice QCD (LQCD) and other LQFTs by stochastically sampling field configurations from appropriately chosen probability distributions and averaging observables over quantum fluctuations. If the probability distribution used for MC sampling is proportional to the contribution of each field configuration to the thermal partition function, then equilibrium thermodynamic observables can be computed from MC ensemble averages.
Sign problems arise in MC calculations when contributions of different field configurations to path integrals have different signs or more generally when path integral contributions are complex and have different phases. For example, the partition function for QCD at non-zero baryon chemical potential has a sign problem and cannot be used to define a probability distribution for MC simulations of non-zero baryon density systems. When a partition function has a sign problem, one can instead MC sample according to a different probability distribution and then attempt to reweight the contribution of each field configuration by the ratio of the desired complex weight to the positive weight used for MC importance sampling. In reweighting approaches to the baryon chemical potential sign problem, the Signal-to-Noise (StN) ratio of the reweighting factor, that is the average reweighting factor divided by the square root of its variance, vanishes exponentially quickly as the spacetime volume is taken to infinity [1][2][3][4][5][6][7][8]. Sign problems can arise for particular observables even when the partition functions does not have a sign problem if different contributions to the path integrals representing these observables have different phases. Ensemble averages of the observable may then have small StN ratios analogous to the StN ratios of complex reweighting factors for non-zero baryon density partition functions.
Baryon and nuclear correlation functions have StN ratios that decrease exponentially at a rate predicted by the moment analysis of Parisi [9] and Lepage [10] when the total baryon number integrated over spacetime is increased.
Many LQCD calculations of baryon and nuclear correlation functions rely on a Golden Window of intermediate source/sink separations where signals are consistent with single-state behavior but this StN problem is not too severe [11][12][13][14][15][16][17]. The baryon correlation function StN problem arises from phase fluctuations that lead to sign problems in correlation functions and is absent when phase fluctuations are ignored [18]. The average values of baryon and nuclear correlation functions decay exponentially faster in the limit of large source/sink separation than their average magnitudes, which demonstrates that exponentially precise cancellations must occur between different contributions to MC ensemble averages and signals an exponentially severe StN problem. LQCD simulations of mesons and baryons with the quark masses tuned to reproduce experimental observables are being used today to predict increasingly complex observables directly from the Standard Model. However, LQCD simulations of nuclei are still generally limited to simulations of small (baryon number B = 1 − 5) nuclei with quark masses tuned to unphysically heavy values . Although physical quark mass simulations of small nuclei are becoming computationally feasible, current methods for calculating nuclear correlation functions require computational resources that grow exponentially with baryon number. These problems motivate the study of alternative approaches to calculating path integrals with sign problems arising from phase fluctuations.
A zero-plus-one-dimensional ((0 + 1)D) complex scalar field StN problem arising from phase fluctuations in correlation functions with non-zero U (1) charge is studied in this work as a toy model of the baryon correlation function StN problem in LQCD. 1 Scalar field phase fluctuations are found to qualitatively resemble the LQCD baryon correlation function phase fluctuations described in Ref. [18] and in particular are shown to be wrapped normally distributed and have exponential StN problems in an analytically tractable approximation where magnitude fluctuations are neglected and phase fluctuations are assumed to be small. Analytically integrating over phase fluctuations using a change of variables similar to the dual lattice variables employed in Ref. [43] enables calculations of correlation functions that avoid sign and StN problems in (0 + 1)D complex scalar field theory, but it remains challenging to extend similar analytic integration methods to more complicated theories such as LQCD [7,[44][45][46][47][48][49]. Instead, a new method is explored in this work in which phase differences are "unwrapped," or numerically integrated over a series of spacetime separations. The resulting unwrapped phases are non-compact random variables rather than circular random variables defined modulo 2π. Moments of unwrapped phase differences can be calculated from positive-definite path integrals that do not have sign problems and do not generically require computational resources that increase exponentially with increasing U (1) charge. Correlation functions can be calculated from moments of unwrapped phases using cumulant expansion techniques similar to those of Ref. [50], although beyond leading order in the cumulant expansion sign and StN problems can re-emerge from differences of cumulants. Phase unwrapping in conjunction with this cumulant expansion allows generic complex correlation functions with non-zero U (1) charge to be represented by series of path integrals without sign problems; however, it is shown below that finding a robust numerical implementation of phase unwrapping is challenging even in (0 + 1)D complex scalar field theory.
The phase unwrapping techniques used here to map series of random phases to series of non-compact random variables are analogous to phase unwrapping techniques used in signal processing, radar interferometry, X-ray crystallography, Magnetic Resonance Imaging (MRI), and other areas of science and engineering [51][52][53][54]. The idea of phase unwrapping correlation functions in LQFT was briefly mentioned in Ref. [18] but its numerical implementation faces challenges arising from large phase fluctuations and was not pursued in detail. This work presents a detailed study of phase unwrapping in (0 + 1)D complex scalar field theory in order to explore its potential applicability in physically relevant higher-dimensional LQFTs. The 1D phase unwrapping algorithms studied here can be immediately applied to three-momentum-projected correlation functions in LQCD and other higher-dimensional LQFTs; however, 1D phase unwrapping algorithms generically suffer from numerical instabilities and do not provide a robust solution to LQFT sign and StN problems. These numerical instabilities are argued to arise from an accumulation of phase unwrapping ambiguities related to large phase jumps that occur with non-negligible probability even on lattices with very fine levels of discretization in (0 + 1)D complex scalar field theory. Multi-dimensional phase unwrapping algorithms are known to avoid analogous numerical instabilities, and more robust phase unwrapping algorithms might be achieved in future investigations of phase unwrapping in multi-dimensional LQFTs.
The remainder of this work is organized as follows. LQFT for (0 + 1)D complex scalar fields is reviewed in Sec. II. After discussing sign and StN problems in free complex scalar field theory in Sec. II A, analytic integration over phase fluctuations and its effect on sign and StN problems is discussed in Sec. II B and is used to derive wrapped phase statistics in Sec. II C that confirm the phase fluctuation origin of the sign and StN problems. Phase unwrapping is introduced in Sec. II D. A cumulant expansion method for relating wrapped and unwrapped phase distributions is introduced in Sec. II E. Numerical studies comparing 1D phase unwrapping and the cumulant expansion, analytic phase integration, and standard MC methods are presented in Sec. III. Conclusions and outlook on future multidimensional phase unwrapping applications are presented in Sec. IV.

II. COMPLEX SCALAR FIELD STATISTICS
A. Sign and Signal-to-Noise Problems Consider a complex scalar field ϕ(t) defined on a uniform lattice of points t = 0, . . . , L − 1 representing a discretized (0 + 1)D flat Euclidean spacetime. Units where the lattice spacing is set to unity are used throughout. Periodic boundary conditions (PBCs) are imposed on ϕ, and L is assumed to be even for simplicity. The free complex scalar field action is then given by where L dependence of the action and other quantities below is left implicit, ϕ(L) ≡ ϕ(0) by PBCs, and the last expression introduces a discrete Fourier transform The partition function can be evaluated by Gaussian integration in polar coordinates ϕ(n) = |ϕ(n)|e iθ(n) , The scalar field propagator is given by In the limit L → ∞ the sum can be replaced by an integral that can be evaluated by contour integral methods [55].
Corrections to this form are exponentially suppressed in L and can be evaluated analytically if the LQFT action in Eq. (1) is replaced by its continuum counterpart. This provides a spectral representation for the propagator, where the L → ∞ complex scalar propagator pole E and residue Z 1;0,1 are given by The free complex scalar field action Eq. (1) has a U (1) symmetry This U (1) symmetry can be used to classify sectors of states in the LQFT Hilbert space that do not mix under (Euclidean) time evolution. The charge of the vacuum is defined to be Q = 0. Field products of the form In operator language, the O Q,2P for all Q, P ∈ Z with P ≥ 0 form a complete basis of interpolating operators for Hilbert space states connected to the Euclidean vacuum by polynomial functions of field operators. Two-point correlation functions involving these operators can be evaluated in terms of the scalar propagator by Wick's theorem as Comparing to the general spectral representation guaranteed by unitarity, the energies of states produced by O Q,2P are seen to be E n ∈ {|Q|E, (|Q| + 2)E, . . . , (|Q| + 2P )E}, and the overlaps Z n;Q,2P of the n-th energy eigenstate can be determined straightforwardly. The full spectrum of the theory includes energies E m = mE for all integers m ≥ 0. Each eigenvalue E m appears in (m + 1) different charge sectors Q = −m, −m + 2, . . . , m − 2, m, signaling an (m + 1)-fold degeneracy of energy eigenstates. Since the action in Eq. (1) is real, e −S is a positive-definite function that can be interpreted as a probability distribution MC methods can be used to produce stochastic samples of complex scalar fields distributed by Eq. (12) using for example the Metropolis algorithm. 2 After generating a MC ensemble of field configurations ϕ i with i = 1, . . . , N and N 1, the scalar field propagator can be determined by approximating the path integral with the ensemble average At fixed t and asymptotically large N , convergence of the sample mean to the true average with 1/ √ N scaling is guaranteed (neglecting MC autocorrelations) by the Berry-Esseen theorem since all moments of P(ϕ) are finite. The 1/ √ N prefactor of Var(G(t)) describes the size of statistical errors at asymptotically large N according to the central limit theorem but may only provide a rough guide to the error scaling for arbitrary P(ϕ) and finite N . It is noteworthy that G(t) is real by Eq. (10), but individual samples ϕ i (t)ϕ i * (0) and G(t) at finite N are complex. At large N the . . , 4 in (0 + 1)D complex scalar field theory. Analytic results valid in the L → ∞ limit from Eq. (6) are shown as red lines. The right plot shows the variance in these ground-state energies. The red lines show the theoretically predicted e 2E Q,0 t scaling. Error bars denote sixty-eight percent confidence intervals determined by bootstrap resampling correlation functions calculated from L source points on N = 5, 000 MC field configurations of complex scalar fields with M 2 = 0.00625 and L = 512 generated using the Metropolis algorithm. This ensemble is denoted C0 below, see Sec. III and Appendix C for more details.
real part of G(t) converges to G(t) as in Eq. (13) and the imaginary part of G(t) converges to zero with similar 1/ √ N scaling. Ensemble average estimators can similarly be constructed for general correlation functions, where O i Q,2P is defined by Eq. (8) with ϕ replaced by ϕ i and 1/ √ N error scaling follows from the Berry-Esseen and central limit theorems. Ground-state energies can be determined from the large-t behavior of the effective masses derived from ensemble average correlation functions, Following standard Parisi-Lepage arguments [9,10], the variance of G Q,2P can be described by a linear combination of correlation functions where G −Q,2P = G Q,2P has been used following Eq. (10). The variance of G Q,2P is related to the variance of G Q,2P by 1/ √ N in the large N limit by the central limit theorem, giving at large N Correlation functions describing sectors with U (1) charge Q = 0 face an exponentially severe StN problem where the exponent is proportional to the U (1) charge of the system. This result is confirmed numerically in MC results shown in Fig. 1.
In LQCD, baryon correlation functions similarly face an exponentially severe StN problem whose exponent is proportional to U (1) B baryon number charge. MC studies indicate that this StN problem is related to the sign problem caused by correlation function phase fluctuations [18]. Analogous features can be seen analytically in free complex scalar field theory correlation functions. A magnitude-phase decomposition of the complex scalar field and an analogous decomposition of the scalar boson propagator and correlation functions can be inserted into the path integral representation of the propagator, Eq. (4), to give Fluctuations of the scalar field phase give scalar boson propagators a sign problem. The path integrand in Eq. (20) is not positive-definite and cannot be interpreted as a probability distribution in MC simulations. 3 Applying a similar decomposition to generic correlation functions gives Free complex scalar field correlation functions have a sign problem if they describe states with U (1) charge Q = 0. Identical considerations apply to the StN problem in Eq. (17), demonstrating that sign and StN problems for complex scalar field correlation functions both arise in the presence of non-zero U (1) charge. The average of correlation function magnitude, depends non-analytically on the Fourier modes ϕ(n) but can be calculated simply in MC studies of (0 + 1)D complex scalar field theory. As shown in Fig and e iΘ Q has both an expectation value and a StN problem of O(e −QEt ) as shown in Fig. 2. Additional StN degradation is present in calculations of excited-state energies. Correlation functions G 0,2P , which include both vacuum and Q = 0 excited state contributions, have qualitatively similar behavior to G 0,1 = |ϕ(t)ϕ(0)| in Fig. 2 with O(1) signal and root-mean-square variance independent of t. After subtracting vacuum contributions, connected Q = 0 correlation functions G conn 0,2P are given by  where ∼ denotes proportionality at large t. Connected Q = 0 correlation functions are O(e −2Et ) in expectation but O(1) sample-by-sample and have a StN problem identical to Q = 2 charged correlation functions. Eq. (24) provides a simple example of a sign problem: expectation values of O(1) random variables must cancel to exponentially increasing precision in order to achieve constant precision in calculations of G conn 0,2P at increasing t. Interpreting the vacuum-subtracted correlation functions G conn 0,2P as belonging to the E n = 0 "charge sector," excited state spectroscopy in free (0 + 1)D scalar field theory can be interpreted as possessing a sign problem and StN problem associated with exponentially small differences of averages that is associated with the presence of non-zero energy "charge." This appearance of a sign problem and StN problem for interpolating operators overlapping on to excited states of the Q = 0 vacuum sector is analogous to the exponentially worse StN problem faced by excited-state correlation functions than ground-state correlation functions in LQCD when excited-state correlation functions are constructed from nonpositive linear combinations of correlation functions with the same quantum numbers.

B. Dual Lattice Variable Sign Problem Solution
Generalizing Eq. (1) to include an arbitrary U (1) invariant potential V (|ϕ|) and inserting the magnitude-phase decomposition of Eq. (18) gives with The partition function can be factored into magnitude and phase contributions as It is shown in Appendix A that the integral over phase fluctuations in Eq. (27) can be performed analytically using a change of variables analogous to the dual lattice variables used in Ref. [43]. Analytically integrating over phase fluctuations provides a positive-definite representation for correlation functions, where the I q are modified Bessel functions, P 0 is a probability distribution suitable for MC sampling of scalar field magnitudes after phase fluctuations have been analytically integrated out, and Z 0 is a normalization coefficient defined by´D|ϕ|P 0 (|ϕ|) = 1. The explicit sum over scalar field phase winding numbers q is included in Eq. (28) in order to avoid "topological freezing" arising with stochastic sampling over winding numbers; see Appendix A for details. Significant contributions to Eq. (28) arise for q = −Q, . . . , +Q but topological charge sectors with |q| > |Q| make sub-dominant contributions that rapidly converge to zero and allow the sum over topological charge sector to be truncated in practical calculations. Given a finite MC ensemble of scalar field magnitude |ϕ i |, i = 1, . . . , N sampled from Eq. (29), correlation functions can be estimated from the corresponding ensemble averages where G dual Q,2P denotes ensemble average calculations of G Q,2P in this dual variables approach. The variance of correlation functions after integrating over dual lattice phase variables is given by where the scaling estimate arises from counting each positive-definite factor of I |Q| /I 0 as O(e −|Q|E ) for t L and 1/ √ N corrections have been neglected. This suggests that charged scalar correlation functions with analytically integrated dual phase variables avoid both the sign problem for importance sampling correlation functions and the O(e −QEt ) StN degradation associated with correlation functions from standard MC methods where the phase is stochastically sampled. One still expects StN degradation arising from numerically estimating the average product of an increasingly large number of variables as t is increased, but this residual StN problem is not associated with estimating a signal that vanishes in the large-t limit and should therefore be much less severe. 4 For numerical verification, the free scalar boson propagator for ensemble C 0 is compared to the propagator determined by MC sampling of the dual representation Eq. (29) with identical parameters M 2 = 0.00625 and L = 512 in Fig. 2. The dual representation provides calculations of the propagator and effective mass with slower StN degradation than the standard representation.
The integrand in Eq. (27) includes products of I |q| functions reminiscent of transfer matrix products appearing in symmetry-projected path integral constructions of Ref. [58][59][60]. Times between the scalar source and sink are associated with factors of I |Q+q| Bessel functions orthogonal to the I |q| Bessel functions associated with the partition function and with the propagator for times outside the source and sink. This suggests that integration over the phase projects the transfer matrix to sectors of definite scalar U (1) charge. Formally, group theoretical projectors are constructed for path integrals by integration over all elements of a symmetry group. Integration over decoupled e i∆ factors associated with each link is equivalent to integration over the group of local U (1) transformations ϕ(t) → e iα(t) ϕ(t), and so dual lattice phase integration acts as a projector to sectors of the Hilbert space with definite U (1) charge. This projection provides the essential mechanism by which dual lattice phase integration avoids the U (1) charged correlation function sign and StN problem.

C. Wrapped Phase Statistics
The magnitude-phase decomposition of the partition function in Eq. (27) shows that for a given scalar field magnitude the phase differences ∆(t) are independent in the L → ∞ limit where the PBC constraint L−1 t=0 ∆(t) = 2πw can be neglected. The L → ∞ distribution for ∆(t) is given from Eq. (27) in terms of κ(t) by This distribution is known as a von Mises distribution and is well-studied in circular statistics [63,64]. The resulting probability distribution describing phase differences Θ(t) = θ(t) − θ(0) as sums of independent von Mises random variables can be expressed as I |n| (κ(t )) I 0 (κ(t )) .
It is difficult to calculate many properties of this probability distribution analytically. The remainder of Sec. II studies a simpler approximation to Eq. (33) where StN ratios can be calculated analytically for correlation function estimators using phase unwrapping introduced in Sec. II D.
A simpler approximation to Eq. (33) can be derived under the assumptions For fine discretizations with M 2 1, the gradient term provides the dominant contribution to the action and Eq. (34) should approximately hold for generic neighborhoods of generic field configurations. Note however that Eq. (34) is not exact in any limit of complex scalar field theory. The non-trivial consequences of relaxing Eq. (34) are explored below by comparing numerical MC results to analytic expressions derived to leading order in Eq. (34), see in particular Figs. 4 and 13. Throughout the remainder of this section ≈ will be used to denote equality to leading order in the small quantities indicated in Eq. (34). In this approximation, phase differences between adjacent lattice sites are identically distributed as well as independent since 4 Multi-level hierarchical integration [57] can be used to exponentially reduce the StN problem associated with sampling products of increasingly many factors in Eqs. (A10)-(A13) as in Refs. [58][59][60][61][62]. This approach has been explored, and for instance a two-level hierarchical integration scheme for calculating correlation functions from Eqs. The assumption ∆ 1 can be used to further simplify Eq. (32). Expanding the cosine to second order in ∆, restoring invariance under ∆ → ∆ + 2πk shifts through explicit summation, and adjusting the overall normalization to enforcé d∆ P(∆) = 1 exactly at this order gives where the second line can be obtained using Poisson summation. Eq. (36) defines the wrapped normal probability describing a normally distributed random variable defined modulo 2π. The wrapped normal and von Mises distributions both approach normal distributions in the limit of small width κ → ∞ and uniform distribution in the limit of large width κ → 0, but the distributions differ at intermediate κ.
The wrapped normal characteristic function is identical to the normal characteristic function, The characteristic function of Θ can be described as a product of characteristic functions of ∆, The probability distribution of Θ is given by a Fourier transform of this characteristic function, Under the assumption of Eq. (34), the scalar boson propagator is given by Comparing this to the large-time spectral representation Eq. (5), the correct ground-state energy and overlap factor are reproduced if The expectation value of the ensemble average correlation function can be calculated in this approximation as Its variance is given by and its StN ratio is It is noteworthy that the full StN problem for the scalar propagator arises at leading order in Eq. (34) where magnitude fluctuations are neglected and phase differences are wrapped normally distributed. Determination of the scalar propagator pole mass from MC sampling phases distributed according to Eq. (39) is equivalent to parameter inference for a wrapped normal distribution with variance 1/κ ≈ 2E. Avoiding large finite sample size errors in wrapped normal parameter inference requires [63] 1 indicating the window of time in which reliable parameter inference is possible has size scaling only as log N . As shown in Fig. 4, Eq. (39) roughly captures the t dependence of phase difference distributions for MC ensemble C 0 but does not provide a precise fit to MC results. The empirical distribution of ∆ is better described by a heavy-tailed wrapped stable distribution than by a wrapped normal distribution. Similar heavy-tailed phase derivatives were seen to arise for baryon correlation functions in Ref. [18], where it was conjectured that these heavy tails arose from nonperturbative strong interaction physics. The appearance of heavy tails in free scalar field theory suggests that they have a more generic origin. The von Mises distribution describing phase derivatives for a fixed field magnitude does not have heavy tails, and so the heavy tails present in MC distributions must arise from integration over magnitude fluctuations. Large phase jumps leading to deviations from Eq. (39) and their relation to magnitude fluctuations are discussed further below. It is also noteworthy that correlation functions for higher charge sector can be computed under the assumptions of Eq. (34) by Eq. (39) as which does not reproduce the linear spectrum of free field theory in Eq. (10). These deficiencies are addressed in Sec. III with numerical MC studies not relying on the assumptions of Eq. (34).

D. Unwrapped Phase Statistics
The results of the last section demonstrate that exponential StN degradation appears in calculations of the average cosines of wrapped normal phase differences. This suggests that to avoid the sign and StN problem, one needs to avoid numerical sampling of circular random variables. For (0 + 1)D complex scalar field theory phase fluctuations can be integrated analytically and the resulting path integrals describing charged scalar correlation functions have positivedefinite weights and a less severe StN problem. Similar methods can be applied to more complicated scalar field theories in more dimensions [43] and also have a long history of application to lattice gauge theory [65]. The search for a more general dual representation of LQCD where properties of finite density matter can be computed with path integrals with positive-definite weights is an active area of ongoing research, see for instance Refs. [7,[44][45][46][47][48][49]. It is also possible to look for path integral deformations or changes of variables that reduce the severity of phase fluctuations. Methods based on Lefschetz thimbles and more general classes of complex path integral deformations have successfully transformed path integrals with phase fluctuations in a variety of LQFTs into (sums of) path integrals where the phase is exactly fixed or at least fluctuating less severely than in the original theory [66][67][68][69][70][71][72][73][74][75][76][77][78][79]. Still, it is an open challenge to find an efficient representation for computing finite-density observables or multi-baryon correlation functions in LQCD that avoids sign and StN problems.
The problems inherent to numerical sampling of circular random variables could be avoided if one could instead numerically sample a non-compact real random variable. Intuitively, one may imagine stochastically sampling a real random variable representing the angular displacement accumulated by the scalar field phase in the interval [0, t] including any 2π revolutions around the unit circle. This alternative has been explored in other areas of science and engineering where circular random variables appear. A variety of "phase unwrapping" techniques have been developed to extract non-compact variables representing angular displacement from numerical samples of compact phases, see Refs. [51][52][53][54] for reviews. An unwrapped phase difference describes a difference between phases at opposite ends of a parametrized path plus 2π times the "winding number" counting the number of full revolutions of the unit circle accumulated along the path. A formal definition of the unwrapped phase dating back to the 1975 homomorphic signal processing of Oppenheim and Schafer [80] can be used algebraically to compute the unwrapped phase of a complex polynomial [54,[81][82][83][84][85][86], while numerical techniques can be used to approximately compute the unwrapped phase from sufficiently finely sampled time series of phases [80,87,88] and higher-dimensional arrays of wrapped phases [89][90][91][92][93][94][95].
A continuous function denoted Arg that describes accumulated phase differences along a 1D path is defined in Appendix B and when applied to correlation functions gives Since Θ(t) = arg(C(t)) is not continuous at branch cut crossings, the fundamental theorem of calculus cannot be directly applied to Eq. (47) if there is a branch cut crossing in the interval [0, t]. By deforming the integration contour to replace branch cut crossings by integrals encircling a neighborhood of the origin (see Appendix B), Eq. (47) can be transformed into an integral over a domain where θ(t) is analytic plus 2π times the total number of oriented branch cut crossings to give The unwrapped phase difference Θ(t) associated with a LQFT propagator therefore differs from the principle-valued or "wrapped" phase difference by 2π times an integer winding number ν(t) that counts the number of oriented branch cut crossings of the propagator in [0, t].
In LQFT as well as in other applications of complex time series, wrapped phases θ(t) are determined directly from "data" by applying arg to complex random variables. The phase unwrapping problem is to determine winding numbers ν(t) that make the unwrapped phase θ(t) a continuous function of t across the branch cuts of θ(t). For a complex time series that samples a smooth function with sufficiently fine resolution, one expects that branch cut discontinuities of θ(t) can be identified and winding numbers can be assigned to keep θ(t) continuous across these branch cuts of θ(t). It is explicitly demonstrated by Itoh in Ref. [88] that the assumption is sufficient to uniquely define winding numbers for a time series of wrapped phases In LQFT, one might hope that Eq. (49) is valid in generic field configurations when the lattice spacing is much smaller than all physical correlation lengths. However, Eq. (47) shows that points with |ϕ(t)| = 0 have infinite d θ/dt even in the continuum. As demonstrated by example in Fig. 5, near-zeros of |C(t)| can occur for (0 + 1)D free complex scalar field theory with M 2 = 0.01. The wrapped phase of the same correlation function is also shown in Fig. 5 along with results for three different (0 + 1)D phase unwrapping schemes: 1. Numerical integration of d θ/dt according to a linear discretization of Eq. (47).
3. Algebraic phase unwrapping of a linear polynomial interpolation of ϕ using the numerically stabilized Strum sequence method of Kitahara and Yamada [54].  49) is violated. Near-zeros therefore produce an accumulating O(π) sensitivity in the unwrapped phase which increases with increasing t unless Eq. (49) holds at all points. The accumulation of errors problem will present numerical difficulties for the simple phase unwrapping schemes explored in Sec. III. These difficulties become more tractable in higher-dimensional phase unwrapping problems [52,89,90], essentially because redundancies in the multi-dimensional gradient of a smooth function provide additional information that can be used to guide unwrapping across regions where Eq. (49) is violated. The theory of multi-dimensional phase unwrapping is much richer than the 1D theory explored below and besides a brief mention in Sec. IV is deferred to future work. The remainder of this section addresses the StN behavior of the unwrapped phase and how properties of its distribution can be used to usefully estimate the average phase cosine. The small fluctuation assumptions of Eq. (34) are used to leading order and phase differences are therefore wrapped normally distributed. It is also assumed that the unwrapped phase θ differs from θ by 2π times an integer winding number ν and therefore that where the wrapping operator W restricts the unwrapped phase to the interval (−π, π]. Wrapped normal phase differences, Eq. (36), can be generated by applying W to a normally distributed unwrapped phase difference, By construction the average cosine of the wrapped and unwrapped phases are identical The sample mean cosine of an ensemble of unwrapped phases could be used to estimate the ground-state energy with identical results and identical StN degradation as calculations based on the sample mean of the wrapped phase cosine. The boson mass can be estimated more efficiently from a MC ensemble of normally distributed unwrapped phase differences by where Θ → − Θ symmetry has been assumed on the basis of Θ → −Θ symmetry (which follows in the infinite statistics N → ∞ limit from reality of correlation functions guaranteed by unitarity). The corresponding estimate of the correlation function is Under the present assumptions this provides an accurate estimate of the correlation function as N → ∞, The variance of G (2) in the N → ∞ limit can be computed similarly, The correlation function computed from normally distributed unwrapped phases therefore has a StN ratio Eq. (58) demonstrates that normally distributed unwrapped phases provide correlation function estimates whose StN ratios decrease polynomially as t −1 rather than exponentially as e −Et as the spacetime volume t containing non-zero U (1) charge is increased.

E. Unwrapped Characteristic Function and Cumulant Expansion
For field configurations violating the small fluctuation assumptions of Eq. (34), as demonstrated to occur even in free field theory in Figs. 4-5, it is necessary to construct an estimator for cos Θ from the unwrapped phase that does not depend on assumptions about the distribution of Θ. The constraint that the winding numbers of the unwrapped phase are integer values, W [ θ] = θ, can be interpreted as a statement that the characteristic functions of the wrapped and unwrapped phase differences agree at every integer, For non-integer n the wrapped and unwrapped characteristic functions can differ. By constraining the unwrapped characteristic function with results not limited to integer n, winding number information present in the unwrapped phase can be incorporated. Once the unwrapped phase characteristic function is fit to numerical results by some method, the mean cosine of the (wrapped or unwrapped) phase is given by evaluating the resultant fit function at n = 1, Cumulant expansion methods similar to those explored in Refs. [18,50,[96][97][98][99][100] can be used to estimate Φ Θ (1) with systematic uncertainties whose size can be assessed by varying the truncation order of the expansion. For a generic complex random variable z with characteristic function Φ z (k) = e ikz , cumulants can be defined as the coefficients of a Taylor series for ln(Φ z ), Equivalent expansions can be constructed that perform a dual expansion in cumulants of the real and imaginary parts of z. The cumulants appearing in Eq. (61) can be related to the moments of z by comparing Taylor series expansions for the exponentials in Eq. (61), If a phase unwrapping algorithm only adds ±1 or 0 to ν(t) at each lattice site as in Eq. (50), then there is a maximum possible value for the unwrapped phase difference |Θ| ≤ 2πt + π and all moments ofΘ are guaranteed to be finite. The finiteness of these moments is sufficient to guarantee that the cumulant expansion converges. Eq. (61) therefore provides a convergent expansion for the characteristic function in the N → ∞ limit. Noting that scalar field propagators are given in terms of the field's log-magnitude and unwrapped phase by an estimator for the scalar boson mass can be defined by In the limits n max → ∞ and N → ∞, Eq. (64) converges to the scalar mass, or to the ground-state energy for twopoint correlation functions in general LQFTs. This consistency requires that W [ θ] = θ and finiteness of the moments of Θ but is otherwise independent of the particular choice of phase unwrapping algorithm used to define Θ. The leading contributions to Eq. (64) are since κ 1 (Θ) and the covariance of R and Θ guaranteed to vanish by Θ → −Θ symmetry. In general LQFTs the magnitude and phase might make very different contributions to the effective mass, and so an arbitrary hierarchy is possible between odd cumulant contributions only involving the magnitude and even cumulant contributions that also involve the phase. In particular, κ 2 ( Θ) dominates κ 1 (R) for free complex scalar field theory and the leading contribution to the cumulant effective mass above is where the ellipses denotes contributes from κ n (R+i Θ) with n ≥ 3. All omitted contributions with n ≥ 3 would vanish in the infinite statistics N → ∞ limit if R and Θ were exactly normally distributed and independent. For distributions with finite moments, contributions from n ≥ 3 provide sub-dominant corrections that will be small for approximately normally distributed R and Θ. The size of these contributions can be assessed in practice by comparing results for E (nmax) with multiple truncation points n max and systematic uncertainties can be assigned based off sensitivity of E (nmax) to the truncation point. Since terms with odd n have vanishing phase contributions by Θ → −Θ symmetry, the convergence pattern of E (nmax) should be expected to strongly depend on whether n max is even or odd and be comparatively smooth as a function of even (odd) cumulant number n max = 2, 4, 6, . . . (n max = 3, 5, 7, . . . ). An analogous expansion to Eq. (64) could be defined for the wrapped phase; however, the P(Θ) approaches a uniform distribution at large times and large cumulants will make sizable contributions to the wrapped analog of Eq. (64). Much faster convergence is expected for Eq. (64) if θ(t) differs from θ(t) by non-zero winding numbers and P( Θ) is approximately normal. Any unwrapping algorithm with W [ θ] = θ will define a Θ such that Eq. (64) is a consistent estimator in the n max → ∞ limit, but different algorithms may have different convergence rates.
Estimators for correlation functions including cumulant expansions of unwrapped phases are constructed by generalizing Eq. (55) as Despite the vanishing of all cumulants with n ≥ 3 in the N → ∞ limit of an exactly normal unwrapped phase distribution, the statistical uncertainties of these higher cumulants increase with increasing n. For large n, the variance of the n-th cumulant will be dominated by the variance of the n-th moment. The large moment n and large statistical ensemble size N behavior of cumulant expansion contributions κ n /n! is therefore determined by the statistical behavior of For normally distributed unwrapped phases, these sample moments have expectation value The variance of Θ can be calculated straightforwardly from Eq. (69) and leads to StN behavior for large moments given by Cancellations between contributions from different moments could lead to a smaller variance for E (nmax) and G (nmax) than that of Θ 2n /(2n)!, but exponentially precise cancellations that would be required to avoid 2 −n StN degradation are not expected to occur at finite N . This suggests that E (nmax) and G (nmax) have StN ratios proportional to √ N 2 −n as in Eq. (70). This suggests that even under the assumptions of Eq. (34), construction of a complete solution to the sign problem using phase unwrapping and the cumulant expansion still requires an extrapolation n max → ∞ where N must be taken exponentially large in n max to remove all truncation errors at fixed statistical precision.
The appearance of such an exponentially hard extrapolation should be expected: phase unwrapping is applicable to generic LQFT correlation functions and the sign problem has been demonstrated to be NP-hard for some quantum systems by Troyer and Wiese [101]. For LQFTs including LQCD, observations of the ubiquity of (complex-)lognormally distributed correlation functions [18,50,[102][103][104][105][106][107] suggest that useful results might be obtained using modest n max despite the exponential difficulty of extrapolating to n max → ∞. Understanding the size of truncation errors in practical calculations and systematic limitations of this method will likely require specific studies for particular LQFTs of interest.
Generic correlation functions in (0 + 1)D complex scalar field theory can be analyzed similarly to the scalar boson propagators above. The wrapped phase for a general correlation function G Q,2P only depends on its U (1) charge Q and not on P and is denoted by The unwrapped phase can similarly be chosen to be independent of P and defined by 5 Note that because Arg is a nonlinear function, Θ Q is not simply related to Θ despite the identity W [Q Θ] = Θ Q . In particular Q Θ is sensitive to branch cut crossings of the variable QΘ with principle domain (−πQ, πQ] rather than branch cut crossings of Θ Q = W [QΘ] with principle domain (−π, π]. With Θ defined by Eq. (50), Q Θ will include jumps of 2πQ at branch cut crossings of Θ rather than jumps of 2π at branch cut crossings of Θ Q and is therefore not equal to Θ Q . With a consistent phase unwrapping of Θ Q,2P ∈ (−π, π], correlation functions and ground-state energies can be estimated with the cumulant expansions For any phase unwrapping algorithm satisfying W [ Θ Q ] = Θ Q,2P , these provide unbiased estimators for correlation functions and effective masses in the dual limit N → ∞ and n max → ∞. In general charge sectors, Θ Q is wrapped normally distributed under the assumptions of Eq. (34) and Θ Q can be consistently defined to be normally distributed with variance 1/κ Q chosen to reproduce the ground-state energy charge Q sector. In analogy to Θ, the correct groundstate energy E Q ≡ E Q,0 is reproduced if the varince of Θ is taken to be 1/κ Q ≈ 2E Q . The StN results of Eq. (58) can therefore be applied to G The moment analysis of Eq. (70) can be applied to Θ Q , and suggests that StN ratios for G decrease polynomially with increasing E Q t but exponentially with increasing n max .
The avoidance of exponential StN degradation with increasing E Q t at fixed order in the cumulant expansion can also be understood in the language of sign problems. Integration over phase fluctuations making non-positive-definite contributions to path integrals is replaced by calculation of the moments of the unwrapped phase, that vanish for odd n by unitarity and are path integrals of positive-definite quantities without sign problems or phase fluctuations for even n. A sign problem can re-emerge beyond leading order in the cumulant expansion from linear combinations of positive-definite moments that enter the cumulant expansion with opposite signs. In particular if E (2) approximates E poorly, then the sum of cumulants could be O(e −E Q t ) while the individual cumulant contributions are O(1) and the full sign problem could re-emerge at full strength.

III. ONE-DIMENSIONAL PHASE UNWRAPPING
Phase unwrapping is shown analytically above to remove exponential StN degradation at fixed order in a cumulant expansion under the assumptions of Eq. (34). Relaxing these assumptions, the probability distribution of phase fluctuations becomes more complicated and an unwrapped phase distribution satisfying W [ θ] = θ cannot be easily found analytically. Numerical MC simulations are used in this section to analyze the accuracy and precision of cumulant expansions involving the log-magnitude and unwrapped phase in (0 + 1)D scalar field theory without the small fluctuation assumptions of Eq. (34).

A. Numerical Phase Unwrapping Schemes
For a field defined on a discrete lattice of points, the unwrapped phase is not uniquely defined without further assumptions that for instance could be based on a discretized definition of smoothness. Precisely defining this smoothness assumption is essential for constructing a numerical phase unwrapping algorithm. The assumptions |∂ t θ| < π and W [ θ] = θ lead uniquely to the path unwrapping algorithm of Eq. (50). This section employs this phase unwrapping algorithm and two variations with alternative smoothness criteria that enforce smoothness on larger distances that a single lattice spacing,

Single-point integration: θ(t) is determined by demanding
as in Eq. (50). This technique assumes a finely sampled lattice, such that the probability density of phase jumps near π is vanishing.

Windowed integration, with window w: θ(t) is determined by demanding
This technique more robustly handles large phase jumps by considering the average of previously unwrapped phases. This locally may allow unwrapped phase jump magnitude to exceed π, but in such a way that global fluctuation is reduced. When w = 1, this reduces to single-point integration.
3. Gaussian-weighted integration, with width σ: θ(t) is determined by demanding The normalization N is fixed by t−1 t =0 N e −(t −t) 2 /(2σ 2 ) = 1. This technique allows smoothly interpolating between integer window sizes by providing a non-integer tunable parameter. When σ ∼ w, we expect the two techniques to perform similarly.
1,0 are identical to those shown in Fig. 6. The right plot shows the average inverse StN of these groundstate energy measurements for a time region t = 10 → 20 as a function of Q for various cumulant expansion truncation orders. Gaussian-weighted integration with σ = 1.41 is used to calculate the unwrapped phase.
Only phase unwrapping algorithms satisfying W [ θ] = θ are considered, since this condition guarantees that the truncated cumulant expansion converges to the exact correlation function in the dual limit of infinite truncation order and infinite statistics.
A time-reversal symmetric integration path is used in which one of the forward integration techniques above is applied to determine the unwrapped phase in [0, L/2] and the corresponding reverse integration technique is applied to determine the unwrapped phase in [L/2 + 1, L]. This symmetric integration path has the advantage of beginning with regions closest to the source where phase gradients are smallest and the probability of unwrapping errors due to physical fluctuations violating |∂ t θ| < π is correspondingly smallest. The t = 0 phase θ(0) = θ(0) is used as an initial condition for unwrapping, although this is irrelevant for correlation functions since they only involve phase differences θ(t) − θ(0). With this scheme the unwrapped phase is discontinuous at L/2 if the wrapped phase is associated with a q = 0 field configuration with non-zero U (1) winding number. Since the same point t = 0 is used for the correlation function source and the initial unwrapping point in this scheme, the unwrapped phase should be separately calculated for each correlation function in a MC ensemble when multiple source points are used with each field configuration. As discussed above, the phase differences Qθ(t) − Qθ(0) for each charge sector Q also need to be unwrapped individually since the unwrapped phase is a nonlinear function of the wrapped phase.

B. Complex Scalar Field MC Ensembles
Nineteen different choices of the (0 + 1)D complex scalar field theory parameters M 2 and λ indicated in Table I are employed to generate a variety of free and interacting MC ensembles. Free-field ensembles A 0 , B 0 , and C 0 are generated with parameters M 2 = 0.1, 0.025, and 0.00625 and serve as toy models for LQFTs with coarse, moderate, and fine lattice spacings. In lattice units, the free-field correlation lengths defined by ξ = 1/E are ξ A0 = 3.175, ξ B0 = 6.331, ξ C0 = 12.652, where the L → ∞ approximation Eq. (6) to the free scalar boson mass E has been used. For each choice of M 2 , the length of the lattice has been rescaled to L = 128, 256, and 512 respectively to enforce M L = 128 √ 0.1 ≈ 40.48 and maintain a roughly constant temporal extent in units of the free-field correlation length. Two additional free-field  ensembles D 0 and E 0 are generated with finer lattice spacing to explore lattice spacing dependence. In addition to the free-field ensembles, a variety of interacting scalar field theory ensembles are generated. Quartic self-interactions are used that correspond to potentials with a variety of couplings λ indicated in Table I. Repulsive couplings λ > 0 are necessary for the action to be bounded from below S(ϕ) ≥ 0 and for the path integral representing the thermal partition function to converge. With λ > 0, the partition function is well defined for M 2 < 0. In higher dimensions, this corresponds to a phase of complex scalar field theory where the U (1) global symmetry is spontaneously broken. In (0 + 1)D at finite L the correlation length is finite in the M 2 < 0 phase but much larger than in the M 2 > 0 phase with the same λ and |M 2 |. Ensembles A ± n , B ± n and C ± n describe interacting scalar field theories with the same |M 2 | as A 0 , B 0 , and C 0 with the sign of M 2 indicated by a superscript and two different values of λ shown in Table I denoted by subscripts n = 1, 2. Ensembles with different |M 2 | but the same λ subscript correspond to choices of λ that keep λL/|M 2 | fixed. Two additional negative M 2 ensembles D − 1 and E − 1 , corresponding to D 0 and E 0 , are also generated for a detailed study of lattice spacing dependence in the interacting case.
To interpret calculations at different M 2 and L as having a fixed physical correlation length and varying lattice spacing a for a (0+1)D field of mass dimension [ϕ] = −1/2, the dimensionless parameters used in the MC calculations should be interpreted as (M a) 2 , L/a, and a 3 λ. This scaling is obtained if the dimensionless parameter a 3 λ is chosen for calculations at different (aM ) 2 and L/a such that λL/M 2 is held fixed. In appropriately rescaled units, the spectrum Ea obtained at different (M a) 2 , L/a, and a 3 λ but fixed λL/M 2 will differ by O(λ) in small-λ perturbation theory. It is tempting to interpret this as a renormalization condition that permits quantitative comparison of results at different M a and L/a; however, (0+1)D complex scalar field theory is non-renormalizable and an infinite number of renormalization conditions need to be imposed to consistently match results for the spectrum of LQFTs with different parameters to all orders in λ. Spectral results from LQFT calculations with different parameters but fixed λL/M 2 should approximately agree for λ 1 and L/M 2 1 but will differ in general by non-universal corrections. This non-universality is an artifact of working in (0 + 1)D and arises even at λ = 0 where it can be understood as arising from higher-derivative operators in a Symanzik-improved lattice action [108]. Since there is no universal continuum limit for (0 + 1)D scalar field theory, no attempt is made to employ non-perturbative renormalization conditions to match observables between ensembles with different parameters. In the M L ≈ 40.48 MC ensembles considered here, the two different values of (λL/M 2 ) n held fixed among A ± n , B ± n and C ± n correspond to (λL/|M 2 |) 1 = 16 and (λL/|M 2 |) 2 = 32 respectively.
In addition to MC ensembles generated using the standard action in Eq. (25) with the parameter choices described above, ensembles are also generated using the analytically phase-integrated dual form of the theory given in Eq. (29)- (30). These dual ensembles only involve MC sampling over the magnitude of the scalar field, and as demonstrated analytically above they have no sign problem and a much milder (though still exponentially severe) StN problem with increasing t. It is verified below that these dual ensembles lead to more precise calculations of ground-state energies E Q,0 in charge sectors Q = 1, . . . , 4. These results are used as a precise check on the accuracy of results obtained using the standard ensembles with and without phase unwrapping. Results show small but statistically significant differences between the low-lying energy levels of ensembles with different parameters but equal λL/M 2 . More details on free-field consistency checks and autocorrelation times can be found in Appendix C.

C. 1D Phase Unwrapping Results
Where the phase varies smoothly, a nearest-neighbor unwrapping scheme accurately captures the variation in phase across the ensemble. Close to large phase jumps with |∂ t θ| > π/2, winding number assignment can vary depending on the phase unwrapping algorithm. Different algorithms will lead to distributions that broaden more or less quickly with t and therefore different low-order cumulant expansion truncation errors. Finding an unwrapping scheme with fast convergence in the cumulant expansion amounts to choosing an algorithm for winding number assignment in neighborhoods of large phase jumps that appropriately tunes the unwrapped phase variance growth controlling E (2) Q,0 . It is empirically found that the single-point integration phase unwrapping method described in Sec. III A that enforces | θ(t) − θ(t − 1)| < π gives poor results for the scalar boson mass across all ensembles. Results do not markedly improve at finer lattice spacing. Statistical precision is generally good for correlation functions and groundstate energies estimated from cumulant expansions truncated at low orders with qualitatively similar t scalings to those shown for Gaussian-weighted integration of C 0 in Figs. 6-7. However, the truncation errors of second-and third-order results for E (n) are large, sometimes an order of magnitude larger than both the statistical uncertainties and central values of standard ensemble average estimates of E. Truncation errors decrease at higher orders in the The green points show E (2) obtained using windowed integration with window sizes w shown on the horizontal axis. Dark green error bars on these points show 68% confidence intervals including statistical uncertainties. The dark purple points and error bars show E (2) and its statistical uncertainties obtained using Gaussian integration. The Gaussian widths σ used are proportional to the window size w with w = 1.65σ for C − 2 (left) and w = 1.62σ for C + 2 (right) empirically found to provide agreement between windowed and Gaussian integration. The lighter purple error bars on both Gaussian and windowed unwrapping points show the extent of the variation in central values of E (2) , E (4) , and E (6) and demonstrate that results tend to converge towards dual ensemble results as nmax is increased and also that the size of truncation errors is sensitive to the unwrapping algorithm parameters used. The red bands show dual ensemble results and statistical uncertainties for comparison. expansion, generally with a pattern of visible decreases at even orders that are sensitive to the shape of the phase distribution, but statistical errors increase dramatically.
The windowed and Gaussian-weighted integration methods that enforce smoothness on scales larger than the lattice scale provide estimates for correlation functions and energies with much smaller truncation errors than single-point integration. The statistical precision of results at various orders in the cumulant expansion is roughly independent of the choice of smearing scale (the window size w or Gaussian width σ used when calculating winding numbers) but the central values of low-order results depend sensitively on the smearing scale. A representative demonstration of this tendency is shown for E (n) for ensembles C + 2 and C − 2 in Figs. [8][9]. Results with Gaussian-weighted integration tend to match results with windowed integration up to a single O(1) constant of proportionality between w and σ. Gaussian-weighted integration can be tuned to interpolate between integer-valued window sizes in this way. An empirical condition relating the phase unwrapping smearing scale to the correlation length or another cost function penalizing large truncation errors could be used to self-consistently define an optimal smearing scale, but it is difficult to assess the systematic errors of optimized estimates without sacrificing precision by going to higher orders in the cumulant expansion.
Empirically, window sizes tuned to reproduce the correlation length in each charge sector w Q ∼ ξ Q = 1/E Q,0 tend to give accurate results for E Q,0 at low orders in the cumulant expansion. Results for the ground-state energies in Q = 1, . . . , 4 charge sectors obtained with optimally tuned Gaussian integration phase unwrapping for the free field ensemble C 0 are summarized in Figs. 6-7. It is noteworthy that second-order truncated cumulant expansion energy estimates E have StN ratios that noticeably decrease with increasing t and with increasing Q. This StN decrease shows less curvature on the log-log scale in Fig. 7 than the exponential StN decrease of E Q,0 , and numerical results are consistent with constant StN at second-order and increasingly high order polynomial StN degradation at increasingly high cumulant expansion truncation order.
The cumulant expansion tends to converge from above or below depending on whether the smearing scale is tuned to be larger or smaller than the physical correlation length. A heuristic explanation for these observations is that unwrapping with an overly small smearing scale is overly sensitive to short-distance fluctuations and erroneously adds winding numbers while unwrapping with an overly large smearing scale penalizes diffusive motion away from physically uncorrelated points and leads to under-broadening. The upper plots use the optimal window sizes for C + 2 and C − 2 indicated in Fig. 8 and show little variation with truncation order at all t. The lower plots use sub-optimal window sizes that include significant truncation errors with nmax = 2 and smaller truncation errors with larger nmax. The red band indicates the dual variables estimate of the scalar mass and its uncertainties. show the variance of the effective mass in the dual lattice variable ensemble and demonstrate exponential variance growth that is significantly less severe than the standard effective mass. The nmax = 2 estimate with phase unwrapping has even less severe variance growth and becomes more precise than the dual variable estimate at large t. Phase unwrapped cumulant effective masses with nmax = 2, 4, 6 show variance growth with downward curvature on the logarithmic scale shown that is consistent with polynomial variance growth, though it is difficult to robustly distinguish high-order polynomial from exponential variance growth numerically.  time-independent feature when estimating the effective mass. In either case, truncation errors are reduced by going to higher order in the cumulant expansion at the cost of decreased statistical precision.
The StN behavior of phase unwrapped ensembles using optimally tuned smearing parameters is shown in Fig. 11. There is very little StN degradation in E (2) . Standard ensemble average correlation functions show exponential StN degradation with the expected O(e −E Q,0 t ) scaling, while dual variable correlation functions show much more mild but likely still exponential StN scaling. For the largest source/sink separations, the dual estimate precision grows to become worse than the precision of E (2) . The limiting factors on the accuracy of low-order results in the cumulant expansion extracted with optimally-tuned phase unwrapping are the systematic uncertainties associated with truncation errors and phase unwrapping parameter tuning, not statistical precision. The precision of E Q,0 results are consistent with precise results from the dual ensembles. Systematic truncation uncertainties determined in this way are significantly larger than statistical uncertainties. The combination of large truncation errors in E (2) Q,0 and significant variance growth with increasing n max prevents phase unwrapped results from providing precise and accurate results for the spectrum of (0 + 1)D complex scalar field theory.
As the lattice spacing is taken much smaller than the physical correlation length, phase differences between neighboring lattices become smaller on average. Under the fixed scalar field magnitude assumption, this probability can be calculated using the von Mises distribution derived exactly for ∂ t θ in Eq. (32) to be P(|∂ t θ| > π − ε) = 2 I 0 (κ)ˆπ π−ε d∆ 2π e κ cos(∆) .
Under small fluctuation assumptions and neglecting excited states, κ ≈ 1/(2E) as in Eq. (41) and κ therefore becomes large as the correlation length becomes large in lattice units. The probability in Eq. (81) vanishes rapidly as κ → ∞ with ε > 0, and so one may expect the probability of large phase jumps in a MC ensemble to vanish as M 2 → 0. However, there is some non-negligible probability that |ϕ| fluctuates to become arbitrarily small even at very small lattice spacing; for example, this occurs due to nearly coincident zero-crossings of the real and imaginary parts of ϕ as they fluctuate from one sign to the other as shown in Fig. 5. The distribution of ∂ t θ in LQFT MC ensembles is given by marginalizing over κ, and non-trivial correlations between the magnitude and phase could lead to significant departures from the expectations of Eq. (81). Such departures are seen in Fig. 13, where the probability of jumps larger that ε = π/2 appears to vanish as M 2 → 0 much more slowly than predicted by Eq. (81). Similar scaling is found in free field theory, interacting field theory with M 2 > 0, and somewhat surprisingly also in the M 2 < 0 regime where the magnitude typically fluctuates about local minima where κ in Eq. (81) is non-zero. The expected number of large phase jumps per field configuration L × P(|∂ t θ| > π/2) is empirically observed to grow as |M 2 | is decreased and L is increased to hold |M |L fixed, suggesting that there is never a physically-relevant regime that is likely to be free of large phase jumps and the phase unwrapping ambiguities associated with them. The result that the number of large phase jumps per configuration grows faster than the number of sites grows as the lattice spacing is decreased is particularly troubling because of an accumulation of errors problem arising in 1D phase unwrapping. If there is a single point t jump with a large phase jump, then a phase unwrapping algorithm might erroneously assign the wrong winding number when the phase unwrapping algorithm reaches t jump . With the forward integration schemes described above, this erroneous winding number at t jump will lead to an error of 2π in the unwrapped phase at all t ≥ t jump . This accumulation of errors problem means that large phase differences between nearest neighbor lattice sites, which might be considered lattice artifacts, lead to errors of size (2π) n in contributions of a MC correlation function to the n-th moment of the unwrapped phase that do not disappear as the lattice spacing is reduced. Further studies are necessary to understand whether this scaling is an artifact of the non-universality of (0 + 1)D complex scalar field theory or a feature that persists in 1D unwrapping of momentum-projected correlation functions of renormalizable LQFTs.

IV. CONCLUSIONS
In (0 + 1)D complex scalar field theory, phase fluctuations distinguish correlation functions in Q = 0 charge sectors from vacuum sector correlation functions. These phase fluctuations result in sign problems for the path integrals representing correlation functions, even though the vacuum sector partition function is positive-definite. A method for avoiding (0 + 1)D scalar field sign problems is introduced that relies on numerically integrating time series of phase differences at a range of source/sink separation using phase unwrapping techniques developed for signal processing and a variety of engineering applications. The non-zero moments of the unwrapped phase distribution can be computed with positive-definite path integrals without sign problems. A cumulant expansion involving moments of correlation function log-magnitudes and unwrapped phases can be used to reproduce the spectrum of (0 + 1)D complex scalar field theory. The numerical results presented here include large systematic truncation errors at low orders in the cumulant expansion and decreased precision as well as the re-emergence of a mild StN problem at higher orders. It is argued that the large truncation errors arise from isolated large phase jumps that lead to errors in the n-th moment proportional to (2π) n at all subsequent times. This accumulation of error problem makes results using a cumulant expansion of the unwrapped phase numerically sensitive to the presence of large phase jumps. Numerical MC studies suggest that in (0 + 1)D scalar field theory the probability of having one or more large phase jumps per lattice extent grows as the M 2 → 0 limit is taken to increase the correlation length, and this accumulation of errors problem leads to large systematic errors even at very fine lattice spacing. This may be due to the non-renormalizability of (0 + 1)D scalar field theory and these investigations should be extended to renormalizable field theories to better understand this issue. The appearance of heavy-tailed phase derivative distributions in free field theory as well as in LQCD baryon correlation functions [18] suggests that these problematic large phase jumps are present in physically relevant LQFTs and are possibly generic features of correlation functions with phase fluctuations. If heavy-tailed phase differences are generic features of LQFT, then high moments of the unwrapped phase sensitive to the tails of the distribution may be noisy and convergence of the cumulant expansion may be slow. Leading-order cumulant expansion results using appropriately tuned phase unwrapping algorithms provide precise approximations to correlation functions that avoid sign or StN problems, but robust applications of phase unwrapping in multi-dimensional LQFTs will require a solution to the 1D accumulation of errors problem and perhaps alternative methods of including corrections from noisy higher-order terms.
Other applications of phase unwrapping provide encouraging results demonstrating that the 1D accumulation of errors problem becomes more tractable in higher dimensions. It was realized in the 1980s that accumulation of phase unwrapping errors along a 1D integration path is a generic problem in the presence of under-sampling [88] but can be avoided in alternative algorithms for 2D phase unwrapping [89,90]. The basic source of greater robustness in higher dimensions is that in 1D only one integration path 6 can be used to connect two points t and t , while in two-and higher-dimensions multiple paths can be used to connect the same points. Assuming that Θ(x, y) = arg[G(x, y)] where G(x, y) is an analytic function, 1D unwrapping will provide identical results for phase integration from y to x that do not depend on the choice of 1D integration contour, see e.g. Ref [52]. Under this analyticity assumption, path-dependence that arises in numerical data must be the result of numerical noise and sampling several 1D unwrapping paths adds error correction through redundancy. Applications in 3D have been found to be even more robust to noise than applications in 2D, suggesting that phase unwrapping generically becomes more robust as the number of dimensions is increased [92,94,95]. A simple argument supporting this idea is that phase unwrapping makes smoothness assumptions informed by nearest neighbor phases, and as the number of dimensions increases the number of nearest neighbors that can be used to inform a phase unwrapping algorithm also increases. However, unlike in other applications, MC samples of LQFT field configurations cannot be assumed to be discrete samples of smooth functions. It is not clear without further studies of multi-dimensional LQFTs whether appropriately smeared configurations calculated on finely discretized lattices will be smooth enough for phase unwrapping algorithms to determine unwrapped phases without ambiguities arising from large phase jumps.
The action for a (0 + 1)D complex scalar field with an arbitrary U (1) invariant, spacetime translation invariant potential energy function V (|ϕ|) can be decomposed into magnitude and phase contributions as where translation invariance has been used to shift field arguments. The partition function for the interacting theory can be similarly decomposed as, The phase integral can be evaluated analytically be introducing dual lattice variables representing the differences between phases at adjacent lattice sites, This transformation is related to dual lattice variable methods that have a long history in lattice gauge theory [65] and can be viewed as a 1D analog of the O(N ) model dual lattice variable transformation introduced in Ref. [43]. To simplify the change of variables, we first use 2π-periodicity to rotate the integration domains for the sequence of θ(t) integrals,ˆπ −π dθ(0) We then change variables from θ(t) to ∆(t) for all t ≥ 1, with trivial Jacobian, Having first rotated the integration ranges, we are left with simple integration bounds for the newly-introduced dual variables,ˆπ −π dθ(0) To completely decouple the integrals, we would like to introduce the final variable ∆(0) = θ(0) − θ(L − 1). The presence of PBCs slightly complicates this transformation by introducing the constraint which can be implemented by a δ-function The compact nature of the phase variables requires this final integral to run over all reals to satisfy PBCs in all winding number sectors. To treat all dual variables on equal footing, we instead handle this sum over winding number sectors directly,ˆ∞ This path integral change of variables from θ to ∆ turns the phase integrals turns into a product of decoupled integrals. This representation allows the integral over phases to be explicitly evaluated as where we have used an integral representation to factorize the δ-function and explicitly integrated the ∆(t) to produce modified Bessel functions of the first kind, I |q| (z). The remaining integrals over |ϕ(t)| cannot be evaluated in closed form for arbitrary V (|ϕ|); however, since I |q| (z) ≥ 0 for z ≥ 0, the form of the partition function given in the final line of Eq. (A10) defines a positive-definite probability density for q, ϕ: The probability distribution P(q, |ϕ|) is a positive-definite, normalizable function that can be used for MC sampling of |ϕ| and q as described below. The integrals over phase variables can similarly be performed analytically for scalar field correlation functions. The general correlation function can first be written in terms of the new dual variables, Inserting this observable into the path integration and explicitly evaluating gives The integrand is again positive-definite and can be interpreted as an integration measure without a sign problem, in contrast to Eq. (20). It is also possible to calculate correlation functions by MC sampling field configurations according to P(q, |ϕ|) and then including the ensemble average of the product of Bessel functions in Eq. (A13) as a reweighting factor. Care must be taken in defining MC updates of q. For instance, a Metropolis scheme in which updates q → q are proposed and then accepted with probability min 1, e −S eff (q,|ϕ|)+S eff (q ,|ϕ|) with will experience "topological freezing". The minimum action q = 0 sector is sampled effectively but q = 0 sectors make O(e −L ) suppressed contributions to the partition function, because they involve products of L small factors L−1 t=0 I |q| I0 , so are scarcely or never present in a finite-N MC ensemble. This is problematic for MC calculations of correlation functions, because q = 0 contributions can provide significant contributions to correlation functions with non-zero U (1) charge. Considering Eq. (A13) for the case of the scalar field propagator G = G 1,0 , the q = −1 sector makes a contribution at large t ∼ L involving the exponentially large product I1 ∼ e L . This situation of exponentially rare MC configurations making exponentially large contributions to observables suggests this MC scheme has an overlap problem where the distribution being importance sampled has poor overlap with the region of configuration space making dominant contributions to observables of interest. Figure 15 plots the integrated autocorrelation time for O = |ϕ| 2 for all standard ensembles used in this work. Figure 16 analogously plots the integrated autocorrelation time for the dual variable ensembles. Autocorrelation times between corresponding standard and dual ensembles are similar. There is significant autocorrelation only on the finest lattice (M 2 = 0.00625), and as such all analyses of methods we introduce applied binning with bin size ≥ 10 on the finest lattice to more accurately estimate errors in the presence of this autocorrelation. Larger bin sizes were used in some cases to further improve the χ 2 /DoF estimates. Tables II-VI indicate the bin and window sizes used for all fits to energy levels required to produce the various spectrum plots in this work.
Instead, MC sampling over q can be replaced by explicit summation over a finite subset of winding numbers that make dominant contributions to particular observables. MC sampling can be performed using the modified probability distribution where Z 0 represents the q = 0 contribution to the partition function and is defined to ensure that´D|ϕ|P 0 (|ϕ|) = 1. Importance sampling with respect to P 0 (|ϕ|) can be performed with local Metropolis update steps of |ϕ(t)| and an accept-reject probability determined by changes in the action S eff (q = 0, |ϕ|). Correlation functions G Q,2P can be calculated from field configurations importance sampled according to P 0 (|ϕ|) by explicit summation over all q, Given a finite MC ensemble of scalar field magnitude |ϕ i |, i = 1, . . . , N sampled from Eq. (A15), correlation functions can be estimated from the corresponding ensemble averages where G dual Q,2P denotes ensemble estimates of G Q,2P in this dual variables approach. Significant contributions to Eq. (A17) arise for q = −Q, . . . , +Q but topological charge sectors with |q| > |Q| make sub-dominant contributions that rapidly converge to zero and allow the sum over topological charge sector to be truncated in practical calculations.

Rel. Error
FIG. 14: In the left plots, we compare free-field Z0;1,0 = |ϕ| 2 estimated using both the standard and phase-integrated ensembles versus the analytical value given in Eq. (6). The three main ensembles A0, B0, and C0, agree with the analytical prediction to the percent level. The auxiliary D0 and E0 ensembles agree to the few percent level. In the right plot, GEVP methods are used to determine the lowest six energy levels in the spectrum of the free theory, rescaled into physical units. The coarsest ensemble A0 exhibits large statistical and systematic uncertainties in fitting, and for Q = 3, 4 no plateau could be reliably fit (indicated by vertical gray lines). Where reliable estimates are possible, the data agree with analytical predictions.

Appendix C: Monte Carlo Ensembles
Ensembles are generated via Metropolis sweeps over the sites in a red-black alternating pattern for efficient execution, with N skip /2 odd and N skip /2 even updates between each measurement. Both the standard 1D complex scalar field action defined in Eq. (25) with the potential Eq. (80) and the analytically phase-integrated dual form of the theory given in Eq. (29)-(30) are used to perform MC calculations using identical values of the parameters M 2 , L, and λ given in Table. I. The phase unwrapping techniques based on smoothed numerical integration of the wrapped phase described above are applied to all correlation functions generated using the standard complex scalar action. The cumulant expansion is then used to estimate correlation functions from sample moments of the corresponding unwrapped phases and log-magnitudes, and a generalized eigenvalue problem (GEVP) is solved to numerically extract the low-lying spectrum of the theory from the resulting correlation function estimates [109].
We perform some checks for ensemble consistency. Eq. (6) gives the non-interacting expectation for Z 1;0,1 = |ϕ| 2 . We can reliably estimate this overlap on our non-interacting lattices and compare against the theoretically predicted value. The left plot of Figure 14 compares the standard ensemble and dual variable ensemble estimates versus the theoretical prediction, finding agreement to the percent level for the three main ensembles A 0 , B 0 , and C 0 , while the auxiliary ensembles (used only for investigation of lattice spacing effects) match the prediction at the few percent level. Eqs. (10) and (11) describe the non-interacting spectrum in terms of Q = 0 ground state energy E = E 0 = 2arcsinh (M/2). There are six low-lying states (energy E i ≤ 4E 0 ), with two states in each of the Q = 0 and Q = 1 channels, and one state in each of the Q = 2 and Q = 3 channels. Figure 14 further demonstrates that our free-field ensembles correctly reproduce this low-lying spectrum to within statistical and systematic fitting errors. Figure 15 plots the integrated autocorrelation time for O = |ϕ| 2 for all standard ensembles used in this work. Figure 16 analogously plots the integrated autocorrelation time for the dual variable ensembles. Autocorrelation times between corresponding standard and dual ensembles are similar. There is significant autocorrelation only on the finest lattice (M 2 = 0.00625), and as such all analyses applied binning with bin size ≥ 10 on the finest lattice to more accurately estimate errors in the presence of this autocorrelation. Larger bin sizes were used in some cases to further improve the χ 2 /DoF estimates. Tables II-VI indicate the bin and window sizes used for all fits to energy levels required to produce the various spectrum plots in this work.
Systematic errors presented in the fit tables indicate the variation in central value as the fit window is offset by up to two lattice points. Systematic errors for the phase unwrapping technique additionally include variation in the central value as higher order cumulants are included. In this work, these errors include variation up to cumulant order 6 (the second subleading order in phase variations), as higher-order cumulants include too much noise to be reliably estimated.