Complex Langevin Dynamics and Supersymmetric Quantum Mechanics

Using complex Langevin method we probe the possibility of dynamical supersymmetry breaking in supersymmetric quantum mechanics models with complex actions. The models we consider are invariant under the combined operation of parity and time reversal, in addition to supersymmetry. When actions are complex traditional Monte Carlo methods based on importance sampling fail. Models with dynamically broken supersymmetry can exhibit sign problem due to the vanishing of the partition function. Complex Langevin method can successfully evade the sign problem. Our simulations suggest that complex Langevin method can reliably predict the absence or presence of dynamical supersymmetry breaking in these one-dimensional models with complex actions.


Introduction
The phenomenon of spontaneous breaking of supersymmetry (SUSY) requires an understanding of non-perturbative aspects of quantum field theory [1,2]. Lattice regularized form of the field theory path integral provides a systematic tool to investigate numerous non-perturbative features of quantum field theories. Monte Carlo methods have been used to reliably extract the physics of such systems. The basic idea behind path integral Monte Carlo is to generate field configurations with a positive definite probability weight, given by the exponential of the negative of the action (in Euclidean spacetime), and then compute the path integral by statistically averaging these importance sampled ensemble of field configurations. However, when the action is complex, for example, when studying QCD at finite temperature and baryon/quark chemical potential, QCD with a theta term, Chern-Simons gauge theories, or chiral gauge theories, the fermion determinant of the theory can be complex, and this feature results in the so-called sign problem. This makes the simulation algorithms based on path integral Monte Carlo unreliable.
There exist very few methods in the literature that can effectively handle or evade this infamous sign problem. Some of these include methods such as analytical continuation, Taylor series expansion [3], methods based on the complexification of the integration variables such as the Lefschetz thimble method [4] and the complex Langevin (CL) method [5][6][7][8][9]. (See Refs. [10][11][12][13][14][15] for simulations based on the Lefschetz thimble method.) The CL method is based on stochastic quantization, and it is a straightforward generalization of the real Langevin method into complexified field configurations.
The central theme of stochastic quantization is that the expectation values of observables are obtained as the equilibrium values of a stochastic process. In Langevin dynamics, this is achieved by evolving the system in a fictitious time direction θ (Langevin time) subject to Gaussian stochastic noise. The idea of stochastic quantization for ordinary field theory systems with real actions is extended to the cases with complex actions, and this overcomes the sign problem by defining a stochastic process, with complexified field variables, using Langevin equations for the complex action [16][17][18]. The use of Langevin equations, rather than importance sampling, waives the restriction of real and positive semi-definite measures. Such a strategy makes the real field variables φ complex, during the Langevin evolution, since the gradient of the action, the drift term, becomes complex.
The complex Langevin equation in Euler discretized form for the field variable φ reads is the action of the theory, the Langevin time is discretized as an integer multiple of the Langevin step size , and η(θ) is a Gaussian noise satisfying the constraints η(θ) = 0 and η(θ)η(θ ) = 2δ θθ . In our simulations, we use real Gaussian stochastic noise to control excursions in the imaginary directions of the field configurations [19]. In the case of real action it can be easily shown that in the large Langevin time limit, θ → ∞, the stationary solution of the Fokker-Planck equation will be reached, and thus ensuring the convergence of the Langevin dynamics to the correct equilibrium distribution. When the action is complex, the situation is challenging. The field configurations are complexified due to the drift term being a complex quantity and thus, we end up with complex probabilities. The CL method became very popular when it was first proposed in the 1980s, but certain problems were encountered shortly after. First was the problem of runaways, where the field configurations would not converge, and the second was the convergence to a wrong limit. In recent years, the CL method has been successfully revived, producing correct results even when the sign problem is severe. When applying the CL method, it is essential that the expectation values computed with a real probability P [φ] along the complex trajectory φ = x + iy are in agreement with the original expectation values with complex weight ρ(φ) = e −S(φ) . That is, . (1.5) Recently, the proof towards convergence to the complex weight in the above equivalence has been an object of thorough study, and correctness criteria have been proposed for the validity of the method [20,21]. In Appendix. A, we discuss these criteria and investigate the reliability of the simulations carried out in this work. The revived interest in this field has led to promising developments in the optimization of the method. It has been noted that the complexification of field variables can drastically influence the Langevin trajectory, making large excursions in the imaginary direction. For instance, it might get closer to unstable directions. An efficient algorithm to evade such a situation involves considering adaptive Langevin step size in the numerical integration of Langevin equations [22]. Ref. [9] provides a pedagogical review on this method, and Ref. [23] is a recent review in the context of the sign problem in quantum many-body physics.
The CL method has been used successfully in the study of various models in the recent past [24][25][26][27][28][29][30][31][32][33]. There have also been studies of spontaneous symmetry breaking (SSB) induced by complex fermion determinant [34][35][36] and a recent study addressing the SSB of SO(D) rotational symmetry based on CL method [37]. In Ref. [38] the authors used CL simulations to observe the Gross-Witten-Wadia (GWW) [39][40][41] phase transition in certain large-N matrix models. In Ref. [42] the authors looked at certain classes of zero-dimensional supersymmetric quantum field theories with complex actions using the CL method. In this work we will investigate dynamical SUSY breaking in an interesting class of supersymmetric actions that are complex and exhibiting invariance under the combine operation of parity and time reversal (PT ) symmetry. This class of actions cannot be studied using traditional Monte Carlo methods.
The paper is organized as follows. In Sec. 2 we briefly introduce supersymmetric quantum mechanics with two supercharges. In Sec. 3 we discuss the lattice regularization of the model. There, we also provide a set of useful observables that can probe SUSY breaking in the system. They include correlation functions and Ward identities. In Sec. 4 we present the simulation results of various models including the PT symmetric model. Conclusions are provided in Sec. 5. In Appendix A.1 we study the reliability of simulations using the Langevin operator and in Appendix A.2 we examine the correctness of simulations by observing the fall-off behavior of the probability distributions of the drift term magnitudes.

Supersymmetric quantum mechanics
Let us consider the action S[φ, ψ, ψ] of a supersymmetric quantum mechanics with a general superpotential W (φ). The model is invariant under two supercharges. The degrees of freedom are a scalar field φ and two fermions ψ and ψ. We take the action to be an integral over a compactified time circle of circumference β in Euclidean time. It has the form where B is an auxiliary field and the derivatives with respect to τ and φ are denoted by a dot and a prime, respectively. Denoting the two supercharges by Q and Q, the SUSY transformations have the form and The supercharges satisfy the algebra We also note that the action can be expressed in Q-and QQ-exact forms. That is, The partition function in path integral formalism is defined by with periodic temporal boundary conditions for all the fields.
In a system with unbroken SUSY, we can consider a Hamiltonian H corresponding to the Lagrangian in Eq. (2.1) having energy levels E n where n = 0, 1, 2, . . . , such that the ground state energy E 0 = 0. Then, the bosonic and fermionic excited states form a SUSY multiplet satisfying the algebra given in Eq. (2.4), with |b 0 being the ground state of the system. Assuming that the states |b n and |f n have the fermion number charges F = 0 and F = 1, respectively, when periodic temporal boundary conditions are imposed for both the bosonic and fermionic fields, Z is equivalent to the Witten index ∆ W (β) [2]. It is easy to see that the partition function defined as does not vanish due to the existence of a normalizable ground state. As a result, the normalized expectation values of observables are well-defined.
The normalized expectation value of the auxiliary field defined as can be used as an order parameter to probe SUSY breaking [43]. (We note that the unpaired state appearing in Eqs. (2.8) and (2.9) need not have to be a bosonic state or a unique state.) The auxiliary field was introduced for the off-shell completion of the SUSY algebra.
The Q transformation of B in Eq. (2.2) and the fact that the ground state is annihilated by the supercharges together imply that in the absence of SUSY breaking the normalized expectation value of the auxiliary field vanishes. However, in the SUSY broken case, we end up in a not-so-trivial situation.
In a system with SUSY spontaneously broken, the Hamiltonian H corresponding to the Lagrangian in Eq. (2.1) has a positive ground state energy (0 < E 0 < E 1 < E 2 . . . ), and the SUSY multiplet is defined as satisfying the algebra given in Eq. (2.4). Differently from the unbroken SUSY case, when SUSY is broken, the supersymmetric partition function vanishes due to the cancellation between bosonic and fermionic states. As a consequence, the normalized expectation values of observables will be ill-defined. We consider the auxiliary field as an observable, and the normalized expectation value can be computed as where the numerator vanishes from Q-supersymmetry ( b n |B|b n = f n |B|f n ). Thus, in a system with broken SUSY, the normalized expectation of auxiliary field admit a 0/0 indefinite form. In Ref. [43] Kuroki and Sugino introduced a regulator that explicitly breaks SUSY and resolves the degeneracy by fixing a single vacuum state in which SUSY is broken. This regulator α (the twist parameter) can be implemented by imposing twisted boundary conditions (TBC) for fermions. That is, ψ(τ + β) = e iα ψ(τ ) and ψ(τ + β) = e −iα ψ(τ ). (2.13) It was shown in Ref. [43] that for a non-zero α, the partition function does not vanish, and the normalized expectation value of the auxiliary field is well-defined. In the limit α → 0, PBCs are recovered and supersymmetry is restored. Thus α regularizes the indefinite form in Eq. (2.12). Now, vanishing expectation value of the auxiliary field in the limit α → 0 suggests that SUSY is not broken, while a non-zero value suggests that SUSY is broken. We incorporate the twist α when we introduce the lattice regularized theory in Sec. 3.1. It is possible to integrate out the auxiliary field using its equation of motion to get the on-shell form of the action Upon using the Leibniz integral rule and discarding the resultant total derivative term, the action takes the form In the above expression, the total derivative term we omitted was (∂φ/∂τ )W (φ). Note that such an omission is only possible in the continuum theory. When we discretize the theory on a lattice, this term does not vanish, and its presence is crucial to ensure the Q-exact lattice supersymmetry. Thus, in our lattice analysis, we will use Eq. (2.15) as the continuum target theory.

Lattice regularized models
We discretize the action given in Eq. (2.15) on a one-dimensional lattice. Let us take the lattice to be Λ, having T number of equally spaced sites with lattice spacing a. The integral and continuum derivatives are replaced by a Riemann sum aΣ and a lattice difference operator ∇, respectively. The physical extent of the lattice is defined as β ≡ T a.
There are several ways to regularize a given theory on a lattice. We will choose the prescription in which the derivatives appearing in the action take the form of a symmetric difference operator where are the forward and backward difference operators, respectively, and i, j represent lattice sites. However, it is known that the symmetric derivative leads to the so-called fermion doubling problem and this, in turn, leads to a non-supersymmetric lattice theory. We can use the Wilson discretization prescription to decouple these extra fermionic modes from the system. The difference operator is modified as where ij = ∇ + ik ∇ − kj is the usual lattice Laplacian and the Wilson parameter r ∈ [−1, 1] / {0} [44]. For one-dimensional derivatives it turns out that the standard choice of r = ±1 yields ∇ W ij (±1) = ∇ ∓ ij , thereby suggesting that the doubling problem can be resolved by simply using forward or backward difference operator. The reason being that for any choice of the lattice difference operator, the theories can be made manifestly supersymmetric upon the addition of appropriate improvement terms corresponding to the discretization of continuum surface integrals [45].
In our analysis, for the standard choice of the Wilson parameter, we follow the symmetric derivative with a Wilson mass matrix suggested in Ref. [46]. The lattice regularized action then takes the form where the quantity Ω i is defined as We can make the variables dimensionless by performing appropriate rescaling. Let us consider the following set of redefinitions for the variables Under these rescalings the action becomes

Theory on a lattice
For convenience, we will not be using the tilde sign on the dimensionless variables; all variables and fields mentioned from now on are understood to be dimensionless. Physical quantities will be labeled differently. The supersymmetry transformations are modified to contain the Wilson mass terms. For a given lattice site k they are given by and The main obstacle that prevents the preservation of exact lattice SUSY is the failure of the Leibniz rule for lattice derivatives. Unlike the continuum action given in Eq. (2.15), the lattice regularized action preserves only the Q supercharge. The Q supersymmetry is broken for T ≥ 2. It can also be shown that the action is only Q invariant. That is, The Q supersymmetry is broken for finite lattice size T because it is not possible to define a corresponding Q invariant transformation on lattice variables such that the algebra {Q, Q} = 2∇ S still holds [47]. However, the Q-exactness is essential and sufficient to kill any SUSY breaking counter-terms and thereby suppress lattice artifacts. See Refs. [43,46,48,49] for more discussions on this.
As mentioned in Sec. 2 for the continuum theory, when SUSY is broken, the partition function vanishes. In that case, the expectation values of the observables normalized by the partition function could be ill-defined. To overcome this difficulty in our lattice regularized theory, we will apply periodic boundary conditions for bosons and twisted boundary conditions for fermions [43,50].
Introducing the twist, we have φ T = φ 0 , ψ T = e iα ψ 0 , ψ T = e −iα ψ 0 . For the case of dynamically broken SUSY, when α = 0, the Witten index vanishes, which dictates that the fermion determinant changes its sign depending on the boson field configurations. This is the sign problem in models with dynamical SUSY breaking.
The partition function given in Eq. (2.6) takes the following form where S α is the lattice regularized action that respects the twisted boundary conditions. We have Let us absorb the mass term into the potential W , and define a new potential Ξ as The action with twisted boundary conditions now takes the form Also the expressions for N i and N i become After integrating out fermions, the fermionic contribution to the partition function given in Eq. (3.14) has the form This is nothing but the determinant of the twisted Wilson fermion matrix W F For periodic boundary conditions (α = 0 case) this is in agreement with the expression obtained in Ref. [46]. The full partition function takes the form Given an observable O, we can compute its expectation value as Note that the gradient of the action has to be computed to update the field configurations in the CL method. The drift term, given by the negative of the gradient of the action, contains the fermion determinant in the denominator, whose zeroes in the complexified space cause the subtlety for the conditions required for the justification of CL method. The dynamical variables may come close to the singularity of the drift term (the singular drift problem). A recent study highlighted that such a problem is not restricted to logarithmic singularities but is rather generic and may arise where the stochastic process involves a singular drift term [51].

Correlation functions
Using the expression given in Eq. (3.24) we can compute the correlation functions. The bosonic and fermionic correlation functions are defined as and 26) respectively, at the site k.
The fermionic correlation function can be shown to be Upon comparison with Eq. (3.24) we define ψ 0 ψ k L α as our Langevin observable corresponding to the fermionic correlator ψ 0 ψ k α . That is, (3.28) Now, for the bosonic correlation function, the computation is rather straightforward. The Langevin observable is the bosonic correlation function itself. For the k-th lattice site,

Ward identities
Another set of observables that would help us in the investigations on SUSY breaking is the Ward identities. For the supersymmetric variation of the fields, Eqs. (3.9) and (3.10), the invariance of the lattice action guides us to a set of Ward identities that connect the bosonic and fermionic correlators. The partition function in Eq. (3.14), upon addition of the source terms (J, θ, θ), becomes It is easy to see that the variation of the partition function under the Q-transformations vanishes upon turning off the external sources. That is, QZ α J, θ, θ = 0. In fact, the variation of any derivative of the partition function with respect to these external source terms also vanishes (upon turning off the sources). Taking the derivative of the partition function with respect to the source terms J j and θ i we get following set of non-trivial supersymmetric Ward identities, (3.31) We will consider to probe the SUSY breaking.

Complex Langevin simulations
Let us look at the relevant observables we used in the simulations. One crucial observable is the expectation value of the auxiliary field Studies have shown that the auxiliary field expectation value can be used as an order parameter to reliably predict dynamical SUSY breaking [42,43,50]. The non-vanishing (vanishing) nature of the auxiliary field indicates that SUSY is broken (preserved) in the system. That is, However, this vanishing nature could be accidental for some models. In [43], higher powers of B were considered to confirm SUSY breaking for zero-dimensional models. Thus, we also analyze other significant observables to confirm SUSY breaking predictions. We next consider the bosonic action S B α . It has been studied that for exact lattice SUSY, the expectation value of the bosonic action is independent of the interaction couplings [52]. Thus, the bosonic action expectation value simply counts the number of degrees of freedom on the lattice [52,53]. That is, S B = 1 2 N d.o.f . In supersymmetric quantum mechanics, it is expected that, S = T , and S B = T /2, where T is the number of sites. Thus we have The third indicator is the equality of the fermionic and bosonic mass gaps. The mass gaps can be extracted either by a cosh ma(t − T 2 ) fit for the t-th lattice site, or a simple exponential fit over say, the first or last T /4 time slices of the respective correlation functions [49,52].
The last set of observables involve the Ward identity. They can be used to confirm exact lattice supersymmetry successfully. See Refs. [46,48,52,54]. We expect that W 1 , given in Eq. (3.32), to hold (not to hold) for theories with SUSY preserved (broken). That is, (4.4) Only in the limit α → 0 we can comment, using the above set of observables, if the system possesses exact lattice SUSY. Since the partition function is a well-defined quantity for models with SUSY preserved, as expected, we were able to compute the normalized expectation values of observables, and hence perform numerical investigations for α = 0 (PBC) case. The issue in working without the twist field arises only in models where SUSY is spontaneously broken, since the partition function vanishes, and normalized expectation values of the observables are ill-defined. Hence in Sec. 4.2 we perform CL simulations for various values of the twist parameter to verify the consistency of our results for the α = 0 case.

Supersymmetric anharmonic oscillator
The model we consider in this section, the supersymmetric anharmonic oscillator, has a real action. Our goal is to put the simulation code to test by comparing our results with those given in Ref. [46].
The model has the potential This model has been investigated in great detail with the help of Hybrid Monte Carlo (HMC) algorithm in the lattice and non-lattice formalisms, respectively, in Ref. [46] and Ref. [55]. It was concluded that SUSY is preserved in this model for a finite value of the coupling. In the case with α = 0, the action is real and the evolution of field configurations in the system is governed by real Langevin dynamics. First, we simulate SUSY harmonic oscillator for physical parameters m phys = 10 and g phys = 0, and α = 0. Simulations were performed for different lattice spacings keeping the (physical) circle size β = 1. Figure 1 shows bosonic (blue triangle) and fermionic (red square) physical mass gaps versus lattice spacing (a) and lattice size (T ). Black dashed line shows the continuum value of SUSY harmonic oscillator mass gaps for the physical parameters m phys = 10, g phys = 0, that is, m exact = 10. We see that boson and fermion masses are degenerate within statistical errors, and furthermore, as lattice spacing a → 0, the common mass gap approaches the correct continuum value. The simulations confirm that the free action has an exact SUSY at finite lattice spacing, which is responsible for the degenerate mass gaps. Now we simulate the SUSY anharmonic oscillator for physical parameters m phys = 10 and g phys = 100, and α = 0. Simulations were performed for different lattice spacings keeping the (physical) circle size β = 1. In table 1 we provide the bosonic and fermionic mass gaps. Here also we have m B phys ≈ m F phys indicating that SUSY is preserved in the model. Table 2 contains the expectation values of the auxiliary field B α and bosonic action S B α . The mean expectation value B vanishes in the simulations and thus indicate exact lattice SUSY. This table also contains the expectation value of the bosonic action S B α . For this model, we observe S B = 1 2 T , and it was independent of physical parameters g phys and m phys , which again suggests SUSY is preserved in this model. Figure 2 shows the bosonic (left) and fermionic (right) correlation functions (used to compute the respective mass gaps) versus lattice site (t) for various lattice size (T ) for the  SUSY anharmonic oscillator. In Fig. 3 we show the bosonic and fermionic physical mass gaps versus lattice spacing (a) and lattice size (T ). Here we also compare our results with those obtained by Catterall and Gregory [46], (m B CG , m F CG ), and find excellent agreement. Black dashed line shows the continuum value of SUSY anharmonic oscillator mass gaps for the physical parameters m phys = 10, g phys = 100 that is m exact = 16.865 [45]. We see that boson and fermion masses are degenerate within statistical errors, and furthermore, as lattice spacing a → 0, the common mass gap approaches the correct continuum value. In Fig. 4 we plot the real part of Ward identity W 1 (left), and its bosonic and fermionic contributions (right), given in Eq. (3.32), versus the lattice site t for lattices with T values. We observe that the respective bosonic and fermionic contributions cancel each other out within statistical errors, and hence W 1 is satisfied. Our results confirm that the SUSY anharmonic oscillator has an exact SUSY, which is responsible for the degenerate mass gaps.

General polynomial potential
As a check of our code we consider the model with a degree-k polynomial potential with real coefficients Ξ (k) = g k φ k + g k−1 φ k−1 + · · · + g 0 .      In these systems, it is well known that SUSY is spontaneously broken if the count of zeroes of the potential is even, that is when Ξ (k) (−∞) and Ξ (k) (+∞) have opposite signs [1]. For simplicity, we assume the form g k = g, g k−1 = · · · = g 2 = 0, g 1 = m and g 0 = gµ 2 . Then, for degree k = 4 and 5, we have Using complex CL method, we confirm the analytical prediction that SUSY is broken for even-(k = 4) and preserved for odd-(k = 5) degree real-polynomial potentials.
For even-degree potential we encounter singular drift problem when α = 0. This shows up in the simulations as the complex Langevin correctness criteria involving the probability of absolute drift not being satisfied. For α = 0, we obtained a power law decay instead of an exponential decay. Therefore we have performed simulations for various non-zero α. Our simulations suggest that above a particular α value the correctness criteria is satisfied and the probability of absolute drift falls off exponentially.
In Fig. 5 we show the expectation value of the auxiliary field B α (left) and the bosonic action S B α (right) on a T = 8 lattice for various α values. We consider the parameter space where CL simulations are justified and take the limit α → 0. The filled data points (red triangles for imaginary part and blue circles for real part) represent expectation values of observables for parameter space where CL can be trusted, while for unfilled data points CL correctness criteria is not satisfied. The lines represent a linear fit of observables for parameter space where CL simulations are justified and solid squares represent values of respective observables in the limit α → 0. These results indicate that B does not vanish. The expectation value of the bosonic action S B = 1 2 T , and it is found to be dependent on the physical parameters. These results indicate that SUSY is broken for the even-degree real-polynomial superpotential.  Table 3. Expectation value of the auxiliary field B α and the bosonic action S B α for the odd-degree superpotential Ξ (5) . The parameters used are g phys = 3, m phys = 1, and µ phys = 4.0.
For the model with odd-degree potential, we observe that for α = 0 the CL correctness criteria are satisfied. In table 3 we provide the simulation results for physical parameter g phys = 3, m phys = 1 and µ phys = 4.0, on lattices with T = 8, 12, 16, and α = 0. Our simulations gave vanishing value for the auxiliary field expectation value B . We also observed that the expectation value of bosonic action is S B = 1 2 T within errors. It is also independent of the coupling g. These results indicate that SUSY is preserved for these models. In Fig. 6 we plot the Ward identities for this model. In these plots, on the left, we show the complete Ward identity (real part), and on the right, we present the respective bosonic and fermionic contributions to the Ward identity. For this model, the bosonic and fermionic contributions cancel each other out within statistical errors, and the Ward identity is thus satisfied. These results indicate unbroken SUSY in this model.

PT -symmetric models
In this section, we simulate an interesting class of complex actions that exhibit, in addition to supersymmetry, PT -symmetry.
Quantum mechanics and quantum field theory are conventionally formulated using Hermitian Hamiltonians and Lagrangians, respectively. In recent years there has been an increasing interest in extensions to non-Hermitian quantum theories [56], particularly those with PT -symmetry [57,58], which have real spectra. Such theories have found applications in many areas such as optonics [59,60] and phase transitions [61,62] Re [<ψ - t k=5: g phys =3.0, m phys = 1.0, µ phys =4.0, α=0.0 been shown that it is possible to carry over the familiar concepts from Hermitian quantum field theory, such as the spontaneous symmetry breaking and the Higgs mechanism in gauge theories, to PT -symmetric non-Hermitian theories [63][64][65]. In Ref. [66], the authors constructed PT -symmetric N = 1 supersymmetric quantum field theories in 3 + 1 dimensions. There, they found that even though the construction of the models are explicitly supersymmetric, they offer a novel non-Hermitian mechanism for soft SUSY breaking. Imposing PT -symmetric boundary conditions on the functional-integral representation of the four-dimensional −λφ 4 theory can give a spectrum that is bounded below [67]. Such an interaction leads to a quantum field theory that is perturbatively renormalizable and asymptotically free, with a real and bounded spectrum. These properties suggest that a −λφ 4 quantum field theory might be useful in describing the Higgs sector of the Standard Model.
We hope that our investigations would serve as a starting point for exploring the nonperturbative structure of these types of theories in higher dimensions.
The models we consider here have the following potential with Ξ (φ) = −ig (iφ) (1+δ) and δ is a continuous parameter. The supersymmetric Lagrangian for this PT -symmetric theory breaks parity symmetry, and it would be interesting to ask whether the breaking of parity symmetry induces a breaking of supersymmetry. This question was answered for the case of a two-dimensional model in Ref. [68]. There, through a perturbative expansion in δ, the authors found that supersymmetry remains unbroken in this model. We investigate the absence or presence of non-perturbative SUSY breaking in the one-dimensional cousins of these models using CL method. Clearly, such an investigation based on path integral Monte Carlo fails since the action of this model can be complex, in general 1 .  Table 4. Expectation values of the auxiliary field B α and the bosonic action S B α for the PTsymmetric models with δ = 2, 4.  Re [<ψ - Im [<ψ - Even δ case We note that when δ = 0 the model becomes the supersymmetric harmonic oscillator discussed in Sec. 4.1.
In table 4 we provide the simulation results for δ = 2, 4. Simulations were performed for physical parameter g phys = 0.5, on lattices with T = 4, 8, 12, and α = 0. We noticed that the auxiliary field expectation value B vanishes. We also observe that the expectation value of bosonic action is S B = 1 2 T within errors. It is also independent of the coupling g. These results indicate that SUSY is preserved in these models.
In Figs. 7 and 8 we show the Ward identities for δ = 2, 4, respectively, on a T = 8 lattice. The full Ward identity W 1 is shown on the left panel, and the real and imaginary parts of the bosonic and fermionic contributions to the Ward identity are shown in the middle and right panels. Our simulations show that the bosonic and fermionic contributions cancel each other out within statistical errors, and hence the Ward identities are satisfied. All these results clearly suggest that SUSY is preserved in models with PT -symmetry inspired δ-even potentials.

Odd δ case
In these models we had to introduce a deformation parameter to handle the singular drift problem. The results are extracted in the vanishing limit of this parameter. models with δ = 0, 2, 4 using Monte Carlo simulations. See Refs. [54,[70][71][72] for other related work on supersymmetric quantum mechanics on the lattice.  Re [<ψ - δ=4.0, g phys =0.5, α=0.0 Im [<ψ - δ=4.0, g phys =0.5, α=0.0  The simulations were carried out for various non-zero µ, and then the δ = 1 model is recovered by taking µ phys → 0 limit. Our simulations suggest that when µ phys is above a particular value the correctness criteria is satisfied and the probability of absolute drift falls off exponentially. We take in account the µ phys parameter space where CL simulations are justified and consider the limit µ phys → 0.
In Fig. 9, we show the expectation values B (left) and S B (right) for a lattice with T = 8 and for various µ phys values. The filled data points (red triangles for imaginary part and blue circles for real part) represent expectation values of observables for the parameter space where CL can be trusted, while for unfilled data points where the CL correctness criteria is not satisfied. Lines represent linear fit of observables for parameter space where CL simulations are justified and the solid squares represent values of respective observables in the lim µ phys → 0.
These simulation results indicate that B vanishes in the limit µ phys → 0. Also, the expectation value of the bosonic action S B = 1 2 T in this limit. It is also found to be independent of the physical parameters. Thus we conclude that SUSY is preserved in the δ = 1 model. Now for the δ = 3 model, inspired by the idea mentioned in Ref. [34] (which was successfully applied in Refs. [36,37]) to handle the singular drift problem, we introduce a where d f is the deformation parameter.
The values of d f are chosen such that CL correctness criteria are satisfied. The δ = 3 model is recovered as d f → 0. Our simulations suggest that above a particular d f value the correctness criteria is satisfied and the probability of absolute drift falls off exponentially.
In Fig. 10, we show B (left) and S B (right) on a T = 8 lattice for various values of d f . The filled data points (red triangles for the imaginary part and blue circles for the real part) represent the expectation values of the observables for the parameter space where CL can be trusted, while the corresponding unfilled data points represent the simulation data that do not satisfy the correctness criteria. The dashed curves represent a linear fit for B data and quadratic fit for S B data. The solid squares represent the values of respective observables in the lim d f → 0 limit. The simulation results suggest that B vanishes in the µ phys → 0 limit. Also, S B = 1 2 T within error bars in this limit. It is also found to be independent of the physical parameters of the model. These results indicate that SUSY is preserved in the δ = 3 model.

Conclusions
In this paper, using complex Langevin method, we have investigated dynamical SUSY breaking in quantum mechanics models with real and complex actions.
When periodic boundary conditions are used in quantum mechanics with broken SUSY the expectation values of observables are ill-defined. We resolved this problem by using twisted boundary conditions [43] in our non-perturbative lattice simulations.
We find that for the case of real actions, SUSY is preserved in models with harmonic and anharmonic oscillator potentials, and with odd-powered polynomial potential. SUSY is dynamically broken for the case of even-powered polynomial potential. For the case of supersymmetric anharmonic oscillator, our simulations reproduced the earlier results obtained through Hamiltonian Monte Carlo [46].
We then moved on to simulate an interesting class of supersymmetric models with complex actions that are PT -symmetric. Our simulations suggested that SUSY is preserved in these models.
We noticed that during the complex Langevin simulations some of these models suffered from the singular drift problem. In order to overcome this difficulty we introduced appropriate deformation parameters in the models, such that the criteria of correctness of complex Langevin simulations are respected. The target theories are recovered by taking the limits in which the deformation parameters go to zero.
The reliability of the simulations were checked by studying the Langevin operator on observables and the exponential fall off of the drift terms. The results of our investigations on the correctness of the simulations are provided in Appendix A. Our conclusion is that the complex Langevin method can be used reliably to probe non-perturbative SUSY breaking in various quantum mechanics models with real and complex actions.
It would be interesting to extend our investigations to models in higher dimensions, especially quantum field theory systems in four dimensions, such as QCD with finite temperature and baryon/quark chemical potentials. Another long term hope would be to apply these methods to PT -symmetric supersymmetric quantum field theories in higher dimensions. We hope that these studies may find applications in fundamental physics.

A Reliability of simulations
In order to monitor the reliability of simulations we use the two recently proposed methods: one tracks the vanishing nature of the Fokker-Planck operator [19,20,73] and the other monitors the decay of the probability distribution of the drift term magnitude [74,75].

A.1 Langevin operator on observables
The observables of the theory O i [φ, θ] at i-th site evolve in the following way where L i is the Langevin operator for the i-th site. It is defined as  Table 6. Expectation values of LB α for the model with odd-degree (k = 5) real-polynomial potential given in Eq. (4.8). The parameters used are g phys = 3, m phys = 1, and µ phys = 4.0.
Once the equilibrium distribution is reached, we can remove the θ dependence from the observables. Then we have = 0, and this can be used as a criterion for correctness of the simulations. This criterion was used to study various interesting models in Refs. [19,20,73].
If we take B i at the i-th site as the observable then we have Note that LB respects translational symmetry on the lattice, and hence we can monitor the value obtained by averaging over all lattice sites. In table 5 we provide the expectation values of LB for the simulations of supersymmetric anharmonic oscillator model discussed in Sec. 4.1. We see that this observable is zero within error bars and thus we can trust the simulations.
In Fig. 11, we show the expectation values LB for various twist parameter values for the model with case of real even-degree (k = 4) polynomial potential (discussed in Sec. 4.2). The filled data points (red triangles for imaginary part and blue circles for real part) represent expectation values of LB for the parameter space where CL can be trusted, when the second criterion, decay-of-the-drift-term, to be discussed in the next section, is applied. The unfilled data points represent the expectation values of LB in the parameter space where the decay-of-the-drift-term criterion is not satisfied.
In table 6 we provide the expectation values of LB for the simulations of real odd-degree (k = 5) polynomial potentials discussed in Sec. 4.2. We see that this observable is zero within error bars and thus we can trust the simulations.
In table 7 we show the expectation values of LB for PT symmetric models with δ = 2 and 4.  Table 7. Expectation value of LB α for the PT -symmetric potentials given in Eq. (4.9) with δ = 2 and 4. The simulation parameters used are β = 1, g phys = 0.5, and α = 0.
In Fig. 12 (left), we show the expectation values for various µ phys values of LB for the simulations of δ = 1 model. The filled (unfilled) data points represent the simulations that does (does not) respect the decay-of-the-drift-term criterion. In Fig. 12, (right) we show LB for various deformation parameter d f for the δ = 3 model. Again, the filled (unfilled) data points represent the simulations that does (does not) respect the decay-of-the-drift-term criterion.

A.2 Decay of the drift terms
The decay-of-the-drift-term criterion was proposed in Refs. [74,75]. There, the authors demonstrated, in a few simple models, that the probability of the drift term should be suppressed exponentially or faster at larger magnitudes to guarantee the correctness of the CL method.
The magnitude of the mean drift is defined as (A.4) In our work, to avoid the singular drift problem, we introduced appropriate deformation parameters in the theory. The final results are obtained after extrapolating to the vanishing limits of deformation parameters.
In Fig. 13, we show the probability distribution P (u) against u for the simulations of supersymmetric anharmonic (left) and harmonic (right) potentials discussed in Sec. 4.1. We see that the drift terms decay exponentially or faster in these models and thus the simulations can be trusted.
In Fig. 14 (left) we show the decay of drift terms for various α values, for the model with even-degree (k = 4) real polynomial potential. We see that the drift terms decay exponentially or faster when α ≥ 0.8 and the simulations can be trusted in this parameter regime. In Fig. 14 (right) we show the decay of drift terms for various α values, for the model with even-degree (k = 5) real polynomial potential. We see that the drift terms decay exponentially or faster when α = 0 and the simulations can be trusted in this parameter regime.
In Fig. 15 we show The drift term decay for the PT symmetric models with δ = 2 (left) and δ = 4 (right). The drift terms decay exponentially or faster when α = 0. Thus the simulations can be trusted in this parameter regime. In Fig. 16 (left) we show the drift term decay for the PT symmetric model with δ = 1, for various µ phys values. The drift terms decay exponentially or faster when µ phys ≥ 0.6. Thus the simulations can be trusted in this parameter regime. In Fig. 16 (right) we show the decay of the drift terms for various values of the fermionic mass deformation parameter d f , for the PT symmetric model with δ = 3. Here we see that the drift terms decay exponentially or faster when d f > 1.0 and thus the simulations can be trusted in this parameter regime.