Complex Langevin Simulations of Zero-dimensional Supersymmetric Quantum Field Theories

We investigate the possibility of spontaneous supersymmetry breaking in a class of zero-dimensional ${\cal N} = 2$ supersymmetric quantum field theories, with complex actions, using complex Langevin dynamics and stochastic quantization. Our simulations successfully capture the presence or absence of supersymmetry breaking in these models. The expectation value of the auxiliary field under twisted boundary conditions was used as an order parameter to capture spontaneous supersymmetry breaking in these models.

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 density or with a theta term, Chern-Simons gauge theories or chiral gauge theories, it is not straightforward to apply path integral Monte Carlo. In these cases we encounter a complex action problem or sign problem. The basic aim of complex Langevin method [1][2][3][4] is to overcome this problem by extending the idea of stochastic quantization for ordinary field theoretic systems with real actions to the cases with complex actions. This also leads to complexification of the real dynamical field variables that appear in the original path integral. We can define a stochastic process for the complexified field variables by Langevin equation with a complex action. Then the expectation values in the original path integral are calculated from an average of corresponding quantities over this stochastic process 1 . See Ref. [11] for a pedagogical review on this method and Ref. [12] for a recent review in the context of the sign problem in quantum many-body physics.
The central theme of stochastic quantization is that expectation values of observables are obtained as equilibrium values of a stochastic process. In Langevin dynamics, this is implemented by evolving the system in a fictitious time direction, τ , subject to a stochastic noise.
We could think of applying Langevin dynamics when the actions under consideration are complex. In such cases, the field variables become complexified during Langevin evolution since the gradient of the action, the drift term, is complex.
In our simulations, we use real Gaussian stochastic noise to tame excursions in the imaginary directions of the field configurations [30][31][32].
For an arbitrary operator O, we can define a noise averaged expectation value where the probability distribution P [φ(τ )] satisfies the Fokker-Planck equation When the action is real, it can be shown that in the limit τ → ∞, the stationary solution of the Fokker-Planck equation will be reached guaranteeing convergence of the Langevin dynamics to the correct equilibrium distribution. When the action is complex we will end up in a not so easy situation. The drift term will be complex and thus if we consider Langevin dynamics based on the above equation we will end up with complexified fields: φ = Reφ + iImφ. We can still consider Langevin dynamics with complex probabilities [4,[33][34][35] but proofs towards convergence to the complex weight, exp(−S), will be non-trivial. The paper is organized as follows. In Sec. II we apply complex Langevin dynamics to a class of zerodimensional bosonic field theories with complex actions, to compute expectation values of correlators and then compare them with analytical results. We discuss supersymmetry breaking in a zero-dimensional model with N = 2 supersymmetry and with a general form of the superpotential in Sec. III. In Sec. IV, using complex Langevin dynamics, we explore supersymmetry breaking in these models with real and complex actions for different forms of superpotentials. In Sec. V we conclude and provide possible future directions. In Appendix. A 1 we study a correctness criterion of our simulations using the Fokker-Planck operator. In Appendix. A 2 we study reliability of our simulations by examining the probability distributions of the magnitude of the drift terms. In Appendix. B we provide the set of simulation data tables.

II. BOSONIC MODELS WITH COMPLEX ACTIONS
Let us consider actions of zero-dimensional quantum field theories derived from a general potential of the form with φ being a real scalar field, g a coupling parameter and δ a real number. A class of (Euclidean) scalar quantum field theories, that are not symmetric under parity reflection, has been investigated in the literature using the above form of the potential [36]. We can, for example, write down a twodimensional Euclidean Lagrangian of the form for a scalar field with mass m. Such theories are very interesting from the point of view that they exhibit non-Hermitian Hamiltonians. Even more interesting is that there is numerous evidence that these theories possess energy spectra that are real and bounded below.
One can think of making the above Lagrangian supersymmetric by adding the right amount of fermions. The supersymmetric two-dimensional Lagrangian takes the form where ψ, ψ are Majorana fermions. This supersymmetric Lagrangian also breaks parity symmetry. It would be interesting to ask whether the breaking of parity symmetry induces a breaking of supersymmetry. This question was answered in Ref. [36]. There, through a perturbative expansion in δ, the authors found that supersymmetry remains unbroken in this model. We could think of performing nonperturbative investigations on SUSY breaking in this model using complex Langevin method. We leave this investigation for future work [37]. (Clearly, a nonperturbative investigation based on path integral Monte Carlo fails since the action of this model can be complex, in general.) Let us consider the 0-dimensional version of the bosonic Lagrangian with m = 0. The Euclidean action is the same as the one given in Eq. (6) where N = 2 + δ.
The partition function of this model is We can look at the k-point correlation functions, G k of this model. We have The one-point correlation function, G 1 can be evaluated as [38] and the two-point correlation function, G 2 as Similarly we can compute higher moments of φ. In Table I we compare our results from complex Langevin simulations for G 1 and G 2 with their corresponding analytical results.
In Fig. 1 we show the complexified φ field configurations on the complex φ R − φ I plane as it evolves in Langevin time. The Langevin time history of G 1 for the case N = 3 is shown in Fig. 2. In Fig. 3 we show the Langevin time history of G 1 and G 2 for the case N = 4.

FIG. 2:
Langevin time history of the field variable (one-point correlation function G 1 ) for the i g 3 φ 3 theory at coupling parameter g = 0.5. Simulations were performed with adaptive Langevin step size ∆τ ≤ 0.02, generation steps N gen = 10 6 and measurements taken every 100 steps. Simulated field configurations are an average of 10 2 such simulation chains with random initialization. Solid and dashed lines represent the exact values.

III. SUPERSYMMETRY BREAKING IN ZERO-DIMENSIONAL FIELD THEORIES
Let us consider a 0-dimensional supersymmetric model. For a general supersymmetric potential, W (φ), the action is given by where φ is a bosonic field, ψ andψ are fermionic fields, and B is an auxiliary field. The prime denotes derivative of the superpotential with respect to φ. There is a symmetry in the above action that exchanges fermionic fields with bosonic fields and this symmetry is known as supersymmetry. We define two independent supersymmetry charges Q and Q corresponding to an N = 2 supersymmetry. This action can be derived from dimensional reduction of a one-dimensional theory, that is, a supersymmetric quantum mechanics with two supercharges. We can see that the above action is invariant under the following supersymmetry transformations and The supercharges Q and Q satisfy the algebra We also note that the action can be expressed in Q-or QQ-exact forms. That is, The auxiliary field B has been introduced for off-shell completion of the supersymmetry algebra. It is possible to integrate out this field using its equation of motion The simulations were performed with coupling parameter g = 0.5, adaptive Langevin step size ∆τ ≤ 0.02, thermalization steps N therm = 10 4 , generation steps N gen = 10 6 and measurements taken every 100 steps. We have used an average of 10 2 such simulation chains with random initial configurations. The table compares these numerically simulated values with the exact results.
It is easy to show that the action is invariant under the two supersymmetry charges The partition function of the model is Completing the square and integrating over the auxiliary field it becomes Integrating over the fermions it takes the form When SUSY is broken, the supersymmetric partition function vanishes. In that case, the expectation values of observables normalized by the partition function could be ill-defined.
The expectation value of the auxiliary field B is crucial in investigating SUSY breaking. It can be evaluated as Thus, in this model, the normalized expectation value of B is indefinite (it is of the form 0/0) when SUSY is broken.
In order to overcome this difficulty we can introduce an external field and then eventually take a limit where it goes to zero. We usually introduce some external field to detect spontaneous breaking of ordinary symmetry so that the ground state degeneracy is lifted to specify a single broken ground state. We take the thermodynamic limit of the theory, and after that, the external field is turned off. The value of the corresponding order parameter then would tell us if spontaneous symmetry breaking happens in the model or not. (Note that to detect the spontaneous magnetization in the Ising model, we use the external field as a magnetic field, and the corresponding order parameter then would be the expectation value of the spin operator.) We will also perform an analogues method to detect SUSY breaking in the system. Introduction of an external field can be achieved by changing the boundary conditions for the fermions to twisted boundary conditions.
A. Theory on a one-site lattice Let us consider the above 0-dimensional theory as a dimensional reduction of a one-dimensional theory, which is a supersymmetric quantum mechanics. The action of the one-dimensional theory is an integral over a compactified time circle of circumference β in Euclidean space. We have the action Here the dot denotes derivative with respect to Euclidean time τ ∈ [0, β]. Note that the Q supersymmetry will not be preserved in the quantum mechanics theory.
Let us discretize the theory on a one-dimensional lattice with T sites, using finite differences for derivatives.
We have the lattice action with n denoting the lattice site. We have rescaled the fields and coupling parameters such that the lattice action is expressed in terms of dimensionless variables. The lattice action preserves one of the supercharges, Q. The Q supersymmetry will not be a symmetry on the lattice when T ≥ 2.
Let us consider the simplest case of one lattice point, that is, when T = 1. The action becomes where φ(1) and ψ(1) are dependent on the boundary conditions. In the case of periodic boundary conditions, the action reduces to Thus the action for the 0-dimensional supersymmetric model with N = 2 supersymmetry is equivalent to the dimensional reduction of a one-dimensional theory (a supersymmetric quantum mechanics) with periodic boundary conditions.

B. Twisted boundary conditions
Now, instead of periodic boundary conditions, let us introduce twisted boundary conditions for fermions (analogues to turning on an external field), with the motivation to regularize the indefinite form of the expectation values we encountered earlier 2 . We have The action in this case has the form We see that supersymmetry is softly broken by the introduction of the twist α In the limit α → 0 supersymmetry is recovered. The partition function is The expectation of auxiliary field B is given by It is important to note that the quantity B α is now well defined. Here, the external field α plays the role of a regularization parameter and it regularizes the indefinite form, B = 0/0, of the expectation value under periodic 2 Twisted boundary conditions were considered in the context of supersymmetric models by Kuroki and Sugino in Refs. [39,40].
boundary conditions and leads to the non-trivial result. Vanishing expectation value of auxiliary field, B α in the limit α → 0 indicates that SUSY is not broken, while a non-zero value indicates SUSY breaking. We can write down the effective action of the model with twisted boundary conditions as The drift term needed for the application of complex Langevin method in Sec. IV has the form

IV. MODELS WITH VARIOUS SUPERPOTENTIALS
In this section, we investigate spontaneous supersymmetry breaking in various zero-dimensional models using complex Langevin method. Wherever possible, we also compare our numerical results with corresponding analytical results.

A. Double-well potential
Let us begin with a case where the action is real. We consider the case when the derivative of the superpotential is a double-well potential where g and µ are two parameters in the theory.
When µ 2 > 0, the classical minimum is given by the field configuration φ = 0 with energy implying spontaneous SUSY breaking. The ground state energy can be computed as the expectation value of the bosonic action at the classical minimum We can also see from SUSY transformations that SUSY is broken in the model.

The twisted partition function is
When α → 0 we have Hence, SUSY is broken for W = g (φ 2 + µ 2 ).
Let us consider the observable The above expression, once evaluated, becomes In Fig. 4 we show our results from Langevin simulations of this model. We show linear and quadratic extrapolations to α → 0 limit in Figs. 5 and 6. The results are tabulated in Table II. The simulation results are in good agreement with the analytical predictions, and strongly suggest that SUSY is broken for this model.
We also consider the case when the derivative of the superpotential is complex, where g and µ are again two parameters in the theory. We show Langevin time history of the auxiliary B field, and linear and quadratic extrapolations to α → 0 limit in Figs. 7, 8 and 9, respectively. The results are tabulated in Table III. We have successfully simulated the complex double-well superpotential using complex Langevin and our results strongly suggest that SUSY is preserved for this model.
The results mentioned above can be partly motivated by classical dynamics, that is, in the absence of stochastic noise. In Fig. 10, we show the classical flow diagrams on the φ R − φ I plane for the above discussed double-well models. The arrows indicate normalized drift term evaluated at the particular field point. In the same figure, we have also shown the scatter plot of complexified field configurations. These plots demonstrate how equilibrium configurations are attained during complex Langevin dynamics.

B. General polynomial potential
Let us extend our analyses to the case where the derivative of superpotential, W , is a general polynomial of degree k, and imaginary (Right) parts of B α against the regularization parameter, α for supersymmetric potential W = g (φ 2 + µ 2 ). Simulations were performed with g = 1 and µ = 2. We have used adaptive Langevin step size ∆τ ≤ 10 −4 , thermalization steps N therm = 10 4 , generation steps N gen = 10 6 and measurements were taken every 100 steps. The dashed red lines are the linear fits to B α in α, and filled red squares are the linear extrapolation values at α = 0. The solid black lines represent the quadratic fits to B α in α, and filled black diamonds are the quadratic extrapolation values at α = 0. The α → 0 limit values obtained from these plots are given in Table II. The twisted partition function is written as  6: Plot of real (Left) and imaginary (Right) parts of B α against the regularization parameter, α for supersymmetric potential W = g (φ 2 + µ 2 ). Simulations were performed with g = 3 and µ = 2. We have used adaptive Langevin step size ∆τ ≤ 10 −4 , thermalization steps N therm = 10 4 , generation steps N gen = 10 6 and measurements were taken every 100 steps. The dashed red lines are the linear fits to B α in α, and filled red squares are the linear extrapolation values at α = 0. The solid black lines represent the quadratic fits to B α in α, and filled black diamonds are the quadratic extrapolation values at α = 0 . The α → 0 limit values obtained from these plots are given in Table II. For the second term in the above equation, assuming the coefficients of the polynomial potential to be real, we have Upon turning off the external field, the first term of Eq. (51) vanishes, hence Thus, for a general polynomial superpotential, W of the degree even (odd), the SUSY is broken (preserved).
The expectation value of the auxiliary B field is  Table III. and imaginary (Right) parts of B α against the regularization parameter α for supersymmetric potential W = ig (φ 2 + µ 2 ). The simulations were performed with g = 3 and µ = 2. We have used adaptive Langevin step size ∆τ ≤ 10 −4 , thermalization steps N therm = 10 4 , generation steps N gen = 10 6 and measurements were taken every 100 steps. The dashed red lines are the linear fits to B α in α, and filled red squares are the linear extrapolation values at α = 0. The solid black lines represent the quadratic fits to B α in α, and filled black diamonds are the quadratic extrapolation values at α = 0. The α → 0 limit values obtained from these plots are given in Table III.
The second term of Eq. (54) vanishes for a polynomial superpotential. ( Since we have twisted partition func- Stable fixed point Unstable fixed point FIG. 10: The scatter plot of field configurations (red dots) and classical flow diagram (arrows) on the φ R − φ I plane. The red dots represent trajectories of the fields during Langevin evolution for superpotential (Left) W = g(φ 2 + µ 2 ) and (Right) W = ig(φ 2 + µ 2 ). In these simulations, we have used g = 1.0, µ = 2.0 and α = 0.4. The first 10 5 points are plotted with measurements taken every 10 2 steps. In both cases, the field start at point (1.5, −1.0) and with the aid of a stochastic noise it drifts towards equilibrium configuration. Filled circles and squares, represent the stable and unstable fixed points, respectively. tion in denominator, this term is not indefinite.) Hence, we have k : even (55) Now, turning external field off, α → 0, W 2 ] = 0 k : even (56) The above expression confirms that SUSY is preserved (broken) for odd (even) degree of derivative of a real general polynomial superpotential.
Let us consider polynomial superpotential with real coefficients. In this case the above argument for SUSY breaking is valid. Later, we will also discuss a specific case of complex polynomial potential. For simplicity we assume that g k = g k−1 = · · · = g 0 = 1, then for k = 3, 4 we have and We have learned from Eq. (53) and (56) that SUSY is broken (preserved) for k = 4 (k = 3). In Fig. 11 we show Langevin time history of B α for the above two polynomial models. We show linear and quadratic extrapolations to α → 0 limit in Fig. 12. The results are tabulated in Table IV. The simulation results are in good agreement with the corresponding analytical predictions. Now, let us consider the case with complex polynomial superpotential. We modify the real double-well potential discussed in the previous section as follows, In this complex potential case, the argument given in Eq. (53) and (56) are not valid. We investigate SUSY breaking using complex Langevin dynamics. In Fig. 13, we show Langevin time history of the auxiliary B field for regularization parameter, α = 0.4. We show linear and quadratic extrapolations to α → 0 limit in Figs. 14 and 15. The results are tabulated in Table V. Our simulation results imply that expectation value of auxiliary field, B α does not vanish in the limit, α → 0. Hence SUSY is broken in this model.

C. PT -symmetric models inspired δ-potentials
Let us consider the superpotential which is the same as the one we considered earlier for the case of the bosonic models. Simulations were performed for superpotential W (φ) = g k φ k + g k−1 φ k−1 + · · · + g 0 with g k = g k−1 = · · · = g 0 = 1. In these simulations, we have used adaptive Langevin step size ∆τ ≤ 5 × 10 −5 , generation steps N gen = 10 7 and measurements were taken every 500 steps.
The twisted partition function takes the form The expectation value of the auxiliary field is Let us consider various integer cases of δ and check whether SUSY is broken or preserved in these cases.
For the case, δ = 0 one can easily perform analytical evaluations. We have the twisted partition function Turning the external field off, α → 0, we get a non-zero value for the partition function implying that SUSY is preserved in the system. Also we have Since B α [δ = 0] = 0, we infer that SUSY is preserved in the theory when δ = 0.
For the case δ = 2, we have the twisted partition function Turning the external field off, α → 0, we get a non-zero partition function FIG. 12: The expectation value, B α against the regularization parameter, α for superpotential W (φ) = g k φ k + g k−1 φ k−1 + · · · + g 0 with g k = g k−1 = · · · = g 0 = 1. The simulations were performed with adaptive Langevin step size ∆τ ≤ 5 × 10 −5 , thermalization steps N therm = 5 × 10 4 , generation steps N gen = 10 7 and measurements taken every 500 steps. The dashed red lines are the linear fits to B α in α, and red dots are the linear extrapolation value at α = 0. The solid black lines represent the quadratic fits to B α in α, and black dots are the quadratic extrapolation value at α = 0. The α → 0 limit values obtained from these plots are given in Table IV. indicating that SUSY is preserved in the system.

The expectation value of the B field is
confirming that SUSY is preserved for the case δ = 2.
One can perform similar calculations for the case δ = 4 and show that SUSY is preserved in the theory.
We simulate the δ-potential using complex Langevin dynamics for δ = 1, 2, 3 and 4. The drift term coming from the δ-potential is In Fig. 16 we show the Langevin time history of the auxiliary B field for δ = 1, 2, 3 and 4. We show linear and quadratic extrapolations to α → 0 limit in Fig. 17 for δ = 1, 3 and Fig. 18 for δ = 2, 4, respectively. The results are tabulated in Table VI and VII. It is clear from  FIG. 14: Real (Left) and imaginary (Right) parts of B α against the regularization parameter, α for supersymmetric potential W = igφ (φ 2 + µ 2 ). Simulations were performed with g = 1 and µ = 2. We have used adaptive Langevin step size ∆τ ≤ 5 × 10 −5 , thermalization steps N therm = 5 × 10 4 , generation steps N gen = 10 7 and measurements were taken every 500 steps. The dashed red lines are the linear fits to B α in α, and red dots are the linear extrapolation value at α = 0. The solid black lines represent the quadratic fits to B α in α, and black dots are the quadratic extrapolation value at α = 0. The α → 0 limit values obtained from these plots are given in Table V. our simulation results that the expectation value of auxiliary field, B α , vanishes in the limit α → 0. Hence we conclude that SUSY is not broken in the model with δ-potential for values of δ = 1, 2, 3, 4.

V. CONCLUSIONS AND FUTURE DIRECTIONS
We have successfully used complex Langevin dynamics with stochastic quantization to investigate supersymmetry breaking in a class of zero-dimensional N = 2 supersymmetric models with real and complex actions.
We looked at double-well superpotential, general polynomial superpotential and also PT -symmetric models inspired δ-potentials. In some cases we were able to cross check the presence or absence of supersymmetry breaking wherever analytical results were available. Our simulations strongly suggest that SUSY is preserved for PTsymmetric models inspired δ-potentials. We have also investigated the reliability of complex Langevin simulations by monitoring Fokker-Planck equation as correctness criterion (in Appendix. A 1) and also by looking at the probability distributions of the magnitude of the drift terms (in Appendix. A 2).
It would be interesting to study complex Langevin dy- and imaginary (Right) parts of B α against the regularization parameter, α for supersymmetric potential W = igφ (φ 2 + µ 2 ). Simulations were performed with g = 3 and µ = 2. We have used adaptive Langevin step size ∆τ ≤ 5 × 10 −5 , thermalization steps N therm = 5 × 10 4 , generation steps N gen = 10 7 and measurements were taken every 500 steps. The dashed red lines are the linear fits to B α in α, and red dots are the linear extrapolation value at α = 0. The solid black lines represent the quadratic fits to B α in α, and black dots are the quadratic extrapolation value at α = 0. The α → 0 limit values obtained from these plots are given in Table V. namics in the above models, generalized to non-Abelian cases, for example with SU (N ) symmetry. Supersymmetry may be restored in the large-N limit of these models. It would also be interesting to explore spontaneous SUSY breaking when δ in the superpotential is a continuous parameter. Other possibilities include extending our investigations to 1-and 2-dimensional models with same superpotentials. These results will appear in an upcoming work [37].

ACKNOWLEDGMENTS
Work of AJ was partially supported by the Seed Grant from IISER Mohali. AK was partially supported by IISER Mohali and a CSIR Research Fellowship (Fellowship No. 517019).

Appendix A: Reliability of complex Langevin simulations
In this section we would like to justify the simulations used in this work. We look at two of the methods proposed in the recent literature. One is based on the Fokker-Planck equation as a correctness criterion and the other is based on the probability distribution of the magnitude of the drift term.

Fokker-Planck equation as correctness criterion
The holomorphic observables of the theory O[φ, τ ] evolve according to [30,31,41] where L is the Langevin operator Once the equilibrium distribution is reached, assuming that it exists, we can remove the τ dependence from the observables. Then we have and this can be used as a criterion for correctness of the complex Langevin method. This criterion has been investigated in various models in Refs. [30,31,41]. The criterion for correctness, in principle, needs to be satisfied for a complete set of observables O[φ], in a suitably chosen basis [31]. It leads to an infinite tower of identities, which as a collection, resembles to the Schwinger-Dyson equations.
For the observable O, as the auxiliary B field, we have We show the Langevin history of the above mentioned correctness criterion, LB, with regularization parameter  Table VIII we provide the simulated values of LB α for superpotential W = g(φ 2 + µ 2 ) with coupling parameter g = 1, 3 and various values of regularization parameter, α. In Tables. IX and X, we tabulate the simulated values of LB α for superpotential W (φ) = −ig(iφ) (1+δ) with coupling parameter g = 0.5 and various values of regularization parameter, α.

Decay of the drift terms
Another method to check the correctness of the complex Langevin dynamics, as proposed in Refs. [42,43], is to look at the probability distribution P (u) of the magnitude of the drift term u at large values of the drift. We have the magnitude of the drift term In Refs. [42,43] the authors demonstrated, in a few simple models, that the probability of the drift term should be suppressed exponentially at larger magnitudes in order to guarantee the correctness of complex Langevin method. However, in the models we investigated in this work we see that the probability distribution falls off like a power law with u, even though we have excellent agreements with corresponding analytical results, wherever applicable. In Fig. 21 we show the probability distribution P (u) against u for the superpotential W (φ) = g(φ 2 + µ 2 ) on a log-log plot. In Fig. 22 we show the probability distribution P (u) of the magnitude of drift term u for superpotential W (φ) = −ig(iφ) (1+δ) on a log-log plot. In both cases we see that the distribution falls off like a power law for large u values. This  Table VI. needs further investigations, and we save it for future work. The simulations were performed with adaptive Langevin step size ∆τ ≤ 5 × 10 −5 , thermalization steps N therm = 5 × 10 4 , generation steps N gen = 10 7 and measurements were taken every 500 steps. The dashed red lines are the linear fits to B α in α, and red dots are the linear extrapolation value at α = 0. The solid black lines represent the quadratic fits to B α in α, and black dots are the quadratic extrapolation value at α = 0. The α → 0 limit values obtained from these plots are given in Table VII.  The expectation values B α obtained using complex Langevin simulations for the model with superpotential W = g φ 2 + µ 2 . In the limit α → 0, B α = 0. Thus SUSY is broken in this model.