{'fraction_non_alphanumeric': 0.029799597705430978, 'fraction_numerical': 0.009933199235143658, 'mean_word_length': 4.365041300293099, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Twisted mass quarks and the phase structure of lattice QCD

The phase structure of zero temperature twisted mass lattice QCD is investigated. We find strong metastabilities in the plaquette observable when the untwisted quark mass assumes positive or negative values. We provide interpretations of this phenomenon in terms of chiral symmetry breaking and the effective potential model of Sharpe and Singleton.


Introduction
As a consequence of (soft) chiral symmetry breaking, nature has arranged itself such that three of the pseudo-scalar mesons are light, with masses around 140 MeV. This lightness of the pion mass becomes important also when we think of numerical simulations in lattice QCD. Approaching the "physical point", at which the pion mass assumes its value as measured in experiment, the algorithms used in lattice simulations suffer from a substantial slowing down [1,2] which restricts present simulations to rather high and unphysical values of the quark mass.
In addition to this slowing-down of the algorithms for Wilson fermions, the quark mass does not act as an infrared regulator allowing thus for the appearance of very small unphysical eigenvalues of the lattice Wilson-Dirac operator. These eigenvalues render the simulations more difficult and sometimes even impossible.
Staggered fermions solve this problem but it is not clear how to use this approach to simulate N f = 2 or odd number of flavours [3]. Overlap fermions [4] also solve the problem but they are computationally very demanding and, unless new algorithms are invented, they are very difficult to use for dynamical simulations.
An elegant way out may be the use of so-called twisted mass fermions [5,6]. This formulation of lattice QCD is obtained when the Wilson term and the physical quark mass term are taken not parallel in flavour chiral space, but rotated by a relative twist angle ω. If the Wilson term is given the usual form, such a chiral rotation leads to a twisted mass parameter µ, in addition to the standard Wilson quark mass m 0 ("untwisted" quark mass). Lattice QCD with a twisted mass was first employed for O(a)-improved Wilson fermions with the nice feature that the improvement coefficients and the renormalisation constants are the same as for O(a)-improved Wilson fermions without twisted mass and hence they did not need to be recalculated [7]. The main advantage of the twisted mass fermions is that the twisted quark mass provides a natural infrared cut-off and avoids problems with accidental small eigenvalues rendering therefore the simulations safe. Of course, the slowing down of the algorithms when approaching small quark masses will remain, although it is expected to be less severe.
Later on it was realized that a full O(a)-improvement of correlation functions can be obtained by using the twisted mass alone without additional improvement terms when, as a special case, m 0 is set to the critical value m crit and the above mentioned twist angle is equal to ω = π/2 [8]. In this way the demanding computation of many improvement coefficients can be avoided rendering the simulations much easier both from a conceptual as well as from a practical point of view.
The Wilson twisted mass formulation has been tested numerically in the quenched ap-proximation already [9]. The results are very encouraging. The O(a) corrections appear indeed to be cancelled and even higher order effects seem to be small, at least for the quantities and the value of the quark mass considered in ref. [9]. A word of caution has to be added at this point. Although, as mentioned above, the twisted quark mass can be decreased towards zero without simulations breaking down due to exceptional configurations, there is an important interplay between the lattice cut-off, Λ = a −1 with a the lattice spacing, and the quark mass m q (see equation (10) below). In the continuum in the presence of spontaneous chiral symmetry breaking the chiral symmetry is not realizedà la Wigner and, as the quark mass goes to zero, the chiral phase of the vacuum is driven by the phase of the quark mass term. The same must be true on the lattice, thus the scaling limit a → 0 should be taken before letting m q → 0. As a result, taking the chiral limit is a numerically delicate matter.
In order to ensure in practice that on the lattice the chiral phase of the vacuum is determined by the quark mass term, proportional to m q , and not by the Wilson term, the lattice parameters should satisfy the order of magnitude inequality [8] This same condition emerges from many different corners of the lattice theory when the physical world is approached. A very simple argument leading to the bound (1) is obtained by comparing the magnitude of the critical Wilson term to that of the quark mass term and requiring the first to be negligibly small compared to the second one, in order to be sure that lattice physics matches the requirements of the continuum theory. From the order of magnitude inequality a(Λ QCD ) 5 ≪ m q (Λ QCD ) 3 , one immediately gets the condition (1). It is important to observe, however, that the less restrictive condition may be sufficient if one is dealing with O(a)-improved quantities. It should be remarked that, since aΛ QCD can be (non-perturbatively) expressed in terms of g 2 0 , eqs. (1) and (2) are actually (order of magnitude) conditions for the values of the dimensionless bare lattice parameters am q and g 2 0 . Contact with dimensionful quantities can be made by comparing simulation data with physical inputs.
What is in practice important is to know for which range of the bare lattice parameters one can avoid troubles from chiral breaking cutoff effects, even if parametrically of order a 2 or higher. This issue has to be settled by numerical investigations aimed at establishing both the structure of the phase diagram of the lattice model in study and the size of residual scaling violations on the physical observables.
In this perspective, twisted mass fermions offer a unique opportunity to explore the phase diagram of Wilson fermions. By fixing the twisted mass parameter µ, one may vary (m 0 − m crit ) from positive to negative values. In this way, the phase diagram of zero temperature lattice QCD can be explored. It should be emphasized that, on large lattices, such an investigation would be very difficult without having µ = 0, since else the algorithms would slow down dramatically approaching the critical quark mass.
In this work we have performed simulations to explore the phase diagram of zero temperature QCD. As it will be shown in the following, we find strong metastabilities in the plaquette expectation value. We determined in both metastable branches a number of quantities such as the (untwisted) PCAC quark mass and pseudo-scalar meson masses. The results presented in this paper are obtained at only one value of β = 5.2, with β related to the bare gauge coupling g 0 by β = 6/g 2 0 . Since the value of β = 5.2 corresponds to a rather coarse value of the lattice spacing (a ≈ 0.16 fm) our work can only be considered as a starting point for a more detailed investigation of the phase diagram of lattice QCD. In particular, the β-dependence of the strength of the observed metastabilities has to be determined. We believe that a qualitative and even quantitative understanding of the phase diagram is a necessary prerequisite for phenomenologically relevant numerical simulations.
The paper is organized as follows. In section 2 we introduce Wilson twisted mass fermions and give our notations. This is followed by a short discussion of the algorithms used. In section 3 we provide our evidence for metastabilities by hysteresis effects and long living metastable states. There, we also show results for a selected set of physical quantities. In section 4, we give a possible interpretation of these results in terms of chiral symmetry breaking and the Sharpe-Singleton effective potential model. We conclude finally in section 5. In the appendix some details of the applied update algorithms are explained.
2 Lattice action and basic variables 2.1 Lattice action Let us start by writing the Wilson tmQCD action as In eq. (3) m 0 is the quark mass parameter and µ is the twisted quark mass parameter. The operator D[U] is given bȳ with r the Wilson parameter which will be set to r = 1 in our simulations. The action as it stands in eq. (3) can, of course, be studied in the full parameter space (m 0 , µ). A special case arises, however, when m 0 is tuned towards a critical bare quark mass m crit . In such, and only in such a situation all physical quantities are, or can easily be, O(a) improved. It is hence natural to rewrite m 0 = m crit +m (5) withm an offset quark mass. The values of m crit need only to be known with O(a) accuracy [8] and can be, for instance, taken from the pure Wilson theory at µ = 0. For standard Wilson fermions usually the hopping parameter representation is taken in the numerical simulations. This representation is easily obtained from eq. (3) by a rescaling of the fields We then obtain the form of the action that is actually used in our simulations Although not needed for the discussion of the numerical data presented below, we give for completeness here the action in the so-called physical basis. This action is obtained by introducing new fields ψ(x) andψ(x) which are related to the fields in eq. (3) by a chiral transformation The action then reads where we have identified

Simulation algorithms
In our numerical simulations we used two different optimized updating algorithms for producing samples of gauge configurations: the Hybrid Monte Carlo (HMC) algorithm with up to three pseudo-fermion fields as suggested in [10,11] and the two-step multi-boson (TSMB) algorithm [12].
In the standard HMC algorithm we used even-odd preconditioning, which in presence of a twisted mass is only a slight modification of the standard preconditioning technique [13]. We give the relevant equations in appendix A.1 of this paper. As a subsequent improvement of the algorithm, we implemented the idea of ref. [10] and used shifted fermion matrices to "precondition" the original fermion matrix. These shifted matrices are treated by introducing additional pseudo-fermion fields. In the shifted fermion matrix we simply used larger values of the twisted mass parameter than the value of µ that is to be simulated. Using two pseudo-fermion fields we experienced a substantial improvement of the HMC algorithm by at least a factor of two. The addition of a third pseudo-fermion field gave only another 10-20% improvement. Again we list the relevant equations, how the shifted matrices are implemented, in appendix A.1. As a further algorithmic trick we used the Sexton-Weingarten leap-frog integrator as proposed in ref. [14].
Our alternative algorithm, the TSMB algorithm [12], is based on the multi-boson representation of the fermion determinant [15]. Optimized polynomial approximations are used, both in the first update step and in the second global accept-reject correction step, for reproducing the dynamical effect of fermions on the gauge field. We apply high order least-square optimization and obtain the necessary polynomials using high precision arithmetics [16]. Concerning the optimization of TSMB for QCD see, for instance, [17].
A useful improvement of the TSMB update algorithms can again be achieved by even-odd preconditioning. This can be implemented in TSMB for twisted mass quarks along the lines of ref. [18]. For the even-odd preconditioning of the TSMB update the flavour indices of the quark fields have to be kept. This means that the multi-boson fields have 24 components per lattice site (2 for flavour, 3 for colour and 4 for Dirac spinor indices). Correspondingly, the polynomials are approximating the function x − 1 2 as in the case of a single Dirac flavour with untwisted quark mass. We give some more details of our even-odd implementation of the TSMB algorithm in appendix A.2.
In the region of light quarks an important part of the numerical effort has to be spent on equilibrating the gauge configuration in a new simulation point. This is particularly relevant in studies of the phase structure where many different points in the parameter space have to be investigated. In case of TSMB the equilibration time is substantially longer than the autocorrelation of relevant physical quantities in equilibrium: on our lattices equilibration can take ten or more times the autocorrelation time of the plaquette observable. The autocorrelation times in equilibrium themselves are similar but most of the time by factors of 2-3 shorter with our twisted quark masses than with untwisted quark masses of similar magnitude. For an approximate formula of the computational cost see ref. [19].
The use of two different optimized update algorithms was very helpful in checking our results. We did not try to obtain a precise performance comparison. Qualitatively, we did not see a noticeable difference in the speed once equilibrium was reached, but the HMC algorithm with multiple pseudofermion fields (MPHMC) turned out to be faster in the equilibration process. In particular, crossing the transition region below and above the metastability region is faster with MPHMC. Nevertheless, the extension in κ of the metastability region is the same with both algorithms.
The data used for preparing the figures in this publication were obtained with MPHMC, except for the upper four panels in fig. 2, which were obtained with TSMB. The thermal cycles in fig. 1 were only run with MPHMC. In the other figures the results of the TSMB runs, whenever performed, were always consistent within errors with the shown MPHMC results.

Numerical results
In this section we give our numerical evidence for the phenomenon of metastability mentioned in the introduction. As a first step and for an orientation we have investigated thermal cycles in the hopping parameter κ. We then discuss metastable states in the plaquette expectation value. Finally we determine quantities such as the pion mass and the untwisted PCAC quark mass in the metastable branches in order to obtain a picture of the physical properties in the different states. In most cases we perform the simulations at a twisted mass aµ = 0.01 but in a few cases we also put aµ = 0, which is possible on the lattice sizes we consider.

Thermal cycles
We started our investigation of the phase diagram of zero temperature lattice QCD by performing thermal cycles in κ while keeping fixed β = 5.2 and the value of the twisted mass parameter aµ. These cycles are performed such that a starting value of κ start is chosen and then κ is incremented, without performing further intermediate thermalization sweeps, until a final value of κ final is reached. At this point the procedure is reversed and κ is decremented until the starting value κ start is obtained back. At each value of κ 150 configurations are produced and averaged over.
In fig. 1 we show three such thermal cycles, performed at aµ = 0, aµ = 0.01 and aµ = 0.1 from bottom to top. In the cycles signs of hysteresis effects can be seen for aµ = 0 and aµ = 0.01 while for the largest value of aµ = 0.1 such effects are hardly visible. Hysteresis effects in thermal cycles may be signs of the existence of a first order phase transition. However, they should only be taken as first indications. Nevertheless, they provide most useful hints for further studies to search for metastable states.

Metastability
Guided by the results from the thermal cycles, we next performed simulations at fixed values of aµ and κ, starting with ordered and disordered configurations, staying again at β = 5.2. In fig. 2 we show the Monte Carlo time evolution of the plaquette expectation value, in most cases on a 12 3 × 24 lattice. For several values of κ we find coexisting branches with different average values of the plaquette. The gap (the "latent heat") appears to be rather large. At κ = 0.1717 we show the history of the plaquette expectation value also on a larger (16 3 × 32) lattice. It seems that the gap in the plaquette expectation value does not depend much on the lattice size, suggesting that the metastability we observe here is not a finite volume effect. In most cases the twisted mass is aµ = 0.01, except for the picture left in the bottom line where it is aµ = 0.
The lifetime of a metastable state, i.e. the time before a tunneling to the stable branch occurs, depends on the algorithm used. In fact, one may wonder, whether the appearance of the metastable states seen in fig. 2 may not be purely an artefact of our algorithms. We cannot completely exclude this possibility but we believe it is very unlikely: we employed two very different kinds of algorithms in our simulations as explained in section 2.2. We observe the metastable states with both of them. We also inter-changed configurations between the two algorithms: a configuration generated with the algorithm A was iterated further with   algorithm B and vice versa. We find that in such situations the plaquette expectation value remains in the state where it has been before the interchange of configurations took place. In addition, as we shall see below, the two states can be characterized by well defined and markedly different values of basic quantities. We therefore conclude that the metastable states are a generic phenomenon of lattice QCD in our formulation.

Pion and quark masses
By selecting separately configurations with high and with low plaquette expectation value, we measured the pion mass and the untwisted PCAC quark mass to study the physical properties in the two metastable states.
We obtained the pseudo-scalar ("pion") mass from suitable correlation functions. These are constructed from the standard composite fields defined in terms of the fieldsψ and ψ in eq. (9): Here τ α , α = 1, 2, 3 are the usual Pauli-matrices in isospin space. The corresponding composite fields in terms of the quark fields χ andχ of eq. (3) are then given by the transformation in eq. (8). For instance, for α = 1, 2 ("charged pions") the pseudo-scalar density has the same form in the χ-basis as in the ψ-basis. Therefore, the mass of the charged pions can be extracted from correlators in the χ-basis in the usual way. The charged axial vector and vector currents are rotated into each other by the angle ω in such a way that at ω = π/2 they are interchanged. (For more details see the literature, e.g. [5,8].) Besides the pion mass, we measured the PCAC quark mass from the axial vector current in the χ-basis: HereÔ ∓ is a suitable operator that we have chosen to be the pseudo-scalar densityÔ ∓ = χ τ ∓ 2 γ 5 χ(x), ∂ * µ is the lattice backward derivative defined as usual and τ ± = τ 1 ± iτ 2 . One can show that in the limit a → 0 the quantity m PCAC χ is asymptotically proportional, through finite renormalization constants, tom.
In fig. 3 we show the pion mass squared in lattice units as function of (2κ) −1 . We observe that the pion mass is rather large and the most striking effect in the graph is that it can have two different values at the same κ. If we consider the quark mass m PCAC χ in fig. 4, we see that in the states with a low plaquette expectation value the mass is positive while for high values of the plaquette expectation it is negative. These quark masses with opposite sign coexist for some values of κ. Plotting the pion mass versus m PCAC χ one obtains fig. 5. Figs. 2-4 clearly reveal that for small enough µ metastabilities show up in the quantities we have investigated, such as m π , m PCAC χ and the average plaquette, if m 0 is close to its critical value. What "small enough µ" means is likely to change with β. Simulations at larger values of β are in progress. As a matter of fact, when m 0 is significantly larger (smaller) than m crit we find m PCAC χ to be positive (negative) and no signal of metastabilities. The remark that metastabilities take place for m 0 close to its critical value will be important both in sect. 4.1 to understand why they affect also a purely gluonic observable such as the plaquette and in sect. 4.3, where it leads to a plausible explanation of the observed metastability phenomena in terms of spectral properties of the lattice tmQCD Dirac matrix (suppression of the "eigenvalue cloud crossing" phenomenon by the fermionic determinant).
The remarks in sect. 4.1 may provide further insight also on the similar metastability phenomena reported in [20] for the N f = 3 untwisted Wilson theory and on the reason why they "disappear" when changing the gluonic action or details, e.g. c SW -value, of the fermionic action. A possible reason is that these changes might shift the range of m 0 where metastabilities appear to values where no data are yet available.

Physical interpretation
The observed strong metastabilities discussed in the previous section clearly suggest that we are working either directly at a first order phase transition or at least very close to it such that we see the remnants of a close-by first order phase transition. With the present data we cannot really differentiate between these two scenarios and in the following we will therefore discuss both of them.

Jump in the plaquette and chiral symmetry breaking
Generally speaking, a jump in the plaquette as seen in our data can arise owing to the lack of chiral symmetry for chirally non-invariant formulations of lattice QCD. The argument relies on the key observation that, when working with chirally twisted Wilson fermions, there are two distinct sources of chirality breaking. The first source is the combination of the untwisted Wilson and mass terms    The second source of chirality breaking is the twisted mass term µχiγ 5 τ 3 χ. As pointed out in section 2, one may trade the bare parameters m 0 and µ in eq. (3) for the equivalent bare parameters m q and ω of eq. (9). The latter are best suited to discuss the connection with continuum QCD physics, as ω is an unphysical parameter, while m q represents the bare quark mass. Assuming spontaneous chiral symmetry breaking in infinite volume, the pion mass squared is expected to vanish linearly in m q (up to lattice artifacts) as m q → 0. Moreover in the continuum limit the physical scalar condensate is expected to show a discontinuity and changes sign as m q passes through zero: where by [ψψ] R we mean the appropriately subtracted and renormalized scalar density. We recall that for ω = 0 this is a non-trivial linear combination ofχχ,χiγ 5 τ 3 χ and the constant field (see below for details).
In order to make contact with the observed metastability phenomena in the regime of spontaneous chiral symmetry breaking, two further remarks are important: 1. at non-zero lattice spacing the twisted mass term µχiγ 5 τ 3 χ induces the twisted condensate [χiγ 5 τ 3 χ] R , while the untwisted mass termsχM[U]χ of eq. (13) determines the untwisted condensate [χχ] R .

the local plaquette field
admits on the basis of lattice symmetries an operator expansion of the form with [...] sub denoting a subtracted, multiplicatively renormalizable, operator and F the continuum gauge field strength tensor. The plaquette expectation value P (r, am 0 , aµ) can be correspondingly written in the form The important point about the representation (17) is that it shows that P is actually sensitive to the value of the subtracted condensates [χiγ 5 τ 3 χ] sub and [χχ] sub .
Before continuing it is useful to pause a moment and discuss the structure of eqs. (16) and (17) and nature of the various terms appearing in it.
• We first notice that the contributions from the identity and the F · F operator are put together within a square parenthesis in eqs. (16)- (17) to remind us that there is no unambiguous way to subtract from the latter its power divergent mixing with the identity. Ultimately this is due to the fact that, unlike the chiral condensates, the vacuum expectation value of F · F is not an order parameter of any symmetry.
• For the reason we have just recalled, it is instead perfectly possible to unambigously define, in the massless limit, multiplicatively renormalizable operators [χχ] sub and [χiγ 5 τ 3 χ] sub , by following the procedure outlined in refs. [21]. More generally, such quark bilinears can be defined as finite operators even at non-vanishing masses, though not uniquely. This can be done by setting, for instance with the dimensionless coefficient functions C S 0 and C P determined at some finite spacetime volume V = V 0 by the conditions Both the coefficients C S 0 and C P admit a finite polynomial expansion in am and aµ (actually in (aµ) 2 for parity reasons).
• In terms of [χχ] sub and [χiγ 5 τ 3 χ] sub , the renormalized scalar density in the physical basis, [ψψ] R , reads where z m = Z P /Z S 0 , Z M = [z 2 m cos 2 ω + sin 2 ω] 1/2 and Z Γ denotes the renormalization constant ofχΓχ in the standard Wilson regularization computed in a mass independent renormalization scheme 1 . Consistently with the general arguments given above, we remark that only the leading a −3 divergent subtraction is uniquely fixed by the symmetries of the theory (WTI's and spurionic transformations). Consequently these properties can be used to make the chiral scalar condensate,ψψ, multiplicative renormalizable in the massless limit, by defining it, e.g. as the Wilson average over the expectation values computed with opposite values of the coefficient of the Wilson term [22,8].
After this little digression let us go back and discuss the implications of eq. (17). If we are on the lattice and take the action of eq. (3) for values of µ orm much larger than O(aΛ 2 QCD ), the condensates [χiγ 5 τ 3 χ] sub or [χχ] sub are expected to show no metastability and thus the same should be true for the plaquette expectation value. However, if µ is smaller than O(aΛ 2 QCD ) the physical scalar condensate signaling spontaneous chiral symmetry breaking is not simply given by [χiγ 5 τ 3 χ] sub , but has in general also an untwisted component, [χχ] sub . Both components have an impact on the value of the plaquette (see eq. (17)). Whenm passes from positive to negative values the expectation value of the untwisted operator [χχ] sub should also change sign and, at non-vanishingly small values of µ, eventually become very small for almost critical values of m 0 . In this situation, owing to the presence of the chiral symmetry breaking term (13)

Effective potential model
The scenario of a jump in the scalar condensate for Wilson fermions on the lattice has actually been given already some time ago by Sharpe and Singleton [23]. As it has been shown in that work, the phase structure of lattice QCD for µ = 0 with Wilson-type quarks can be understood in the low energy chiral theory of pseudo-Goldstone bosons if the influence of leading lattice artifacts of O(a) and O(a 2 ) is taken into account. There are two alternatives: either there exists an Aoki phase [22] or there is a first order phase transition between the phases with positive and negative quark mass and the Aoki phase does not exist.
The relevant part of the effective potential is written in [23] as Here A denotes the flavour singlet component of the SU(2) matrix valued field Σ in the low energy effective chiral Lagrangian: Because of the relation 1 = A 2 + 3 r=1 B r B r the variable A is constrained to lie between -1 and +1 inclusive. In the vicinity of the critical quark mass the constant c 2 = O(a 2 ) and the other parameter c 1 is proportional to the bare quark mass (in our notations c 1 ∝m).
In order to find the ground state ("vacuum") the effective potential has to be minimized. Without repeating the details of the discussion in ref. [23] let us just summarize the result.
In case of positive c 2 there exists an Aoki phase in the region of bare quark masses defined by −2c 2 ≤ c 1 ≤ 2c 2 . At the boundaries c 1 = ±2c 2 all three pion masses vanish. Inside the Aoki phase the charged pions are massless because they are the Goldstone bosons of spontaneous flavour symmetry breaking but the neutral pion is massive. Outside the Aoki phase (|c 1 | > 2c 2 ) the flavour symmetry is preserved by the ground state and the three degenerate pions are massive (see figure 6).
The other alternative is that c 2 is negative. In this case the flavour symmetry is preserved everywhere but there exists a minimal pion mass because the pion mass is given by At c 1 = 0 the vacuum expectation value jumps from Σ = A = +1 to Σ = A = −1. Since the jump of this "order parameter" happens at non-zero pion mass (i.e. finite correlation length) the thermodynamical description of the behaviour near c 1 = 0 corresponds to a first order phase transition.
An interesting intermediate situation is defined by c 2 = 0. In this case the vacuum expectation value jumps between Σ = A = +1 and Σ = A = −1 at a single first order phase transition point. This limiting case is the ideal situation, when the phase structure in the Sharpe-Singleton model is identical to the expected one in the continuum. It can be characterized either by saying that the Aoki phase has zero extension or that the minimal pion mass is zero (see figure 6). Of course, this behaviour is valid only up to O(a 3 ) effects, neglecting higher orders in the chiral expansion.

Scenarios
Our numerical results reveal that we clearly observe metastabilities in various quantities. Thus our conclusion is that at least for vanishing twisted mass parameter, i.e. for the standard Wilson lattice theory, there is a first order phase transition. For non-vanishing values of µ we can have two scenarios.
The first is that the first order phase transition persists for µ = 0 but sufficiently small in absolute value. For large µ the theory approaches the quenched limit with a constant quark determinant and therefore it is plausible that no phase transition is expected. This scenario suggests that the first order phase transition line in the (m 0 , µ) plane has an end point: the two phases with positive and negative quark masses are analytically connected (see figure  6). The situation is in this sense analogous to the phase structure of the SU(2) fundamental Higgs model (see chapter 6 of [24] and references therein).
The second scenario is that for any non-vanishing value of µ the first order phase transition disappears. In this scenario, when varying m 0 , one passes at some small distance from the first order phase transition at µ = 0 and just feels this close-by phase transition.
We can at present not differentiate between these two scenarios. From the numerical side we would need to know better the µ and β dependence of the metastability phenomena. From the analytical side an analysisà la Sharpe and Singleton including the twisted mass parameter µ is helpful. 3 The first order phase transition between the phases with positive and negative quark masses observed in the previous section is consistent with the no-Aoki-phase alternative (c 2 < 0) of Sharpe and Singleton.
Our exclusion of the Aoki phase is in agreement with the results of a recent paper [25] which suggests that in case of the unimproved Wilson action the Aoki phase is restricted to the region of strong gauge couplings (β ≤ 4.6). Note that in an early paper on QCD thermodynamics with Wilson quarks [26] a first order "bulk" phase transition has also been observed at β = 4.8 which is consistent both with ref. [25] and with our observations. For further numerical work on the Aoki phase, see refs. [27] The rather strong metastability of the two phases with positive and negative quark mass can be understood on the basis of the properties of the eigenvalue spectrum of the (nonhermitean) Wilson-fermion matrix in the twisted mass basis corresponding to eq. (3). For zero twisted mass (µ = 0) at small positive quark masses there is a "cloud" of eigenvalues close to the origin near the real axis. (For a numerical study see section 4 of [17].) In order to reach negative quark masses this "cloud" has to cross near the origin to the other side with negative real parts. This eigenvalue cloud crossing is strongly suppressed by the zero of the determinant. This, we believe, is the reason at the microscopical level for the observed strong metastability. For non-zero twisted mass there is a strip of width 2|µ| around the real axis where there are no eigenvalues. If this strip is wide enough the eigenvalues are sufficiently far away from the origin and the first order phase transition disappears.
As it was already emphasized in [23], the sign of the coefficient c 2 in the low energy pion effective potential is not universal, it depends on the way the action is discretized. Therefore a clever choice of the lattice action may weaken the first order phase transition and, for instance, decrease the minimal pion mass at it. Previous results of the JLQCD Collaboration [20] support the conjecture that changing the gauge action alone has an important effect. If, indeed, one could find some parameter in the lattice gauge action which at some value would change the sign of c 2 an appealing possibility would be to tune the lattice action to this value. The features of a discretization with c 2 = 0 seem to be quite favourable from the point of view of light quark simulations when, up to O(a 2 ), there would be just a single point in the (m 0 , µ) plane with vanishing pion mass -an ideal situation corresponding to the expected phase structure in the continuum.

Conclusion
In this paper we have explored Wilson twisted mass fermions restricting ourselves to simulations at only one value of β = 5.2. By fixing the twisted mass parameter µ and changing the untwisted Wilson quark mass m 0 , or equivalently the hopping parameter κ, we encountered strong metastabilities in the plaquette expectation value, visible both in thermal cycles as well as in long-living metastable states. At the same time, the pion mass does not vanish but has a minimum at a rather large value. The PCAC quark mass m PCAC χ in the different metastable branches is positive for the branch with low plaquette expectation value and it is negative for the branch with high plaquette expectation value.
The detection of these metastabilities became possible by employing a twisted mass term. Only a non-vanishing value of µ allowed us to cross the critical quark mass. We showed that for lattice theories that break chiral symmetry explicitly the jump of the scalar condensate, when changing the sign of the quark mass, induces a jump of the plaquette expectation value with associated signs of metastability. For µ = 0 these metastabilities find a natural interpretation in the effective potential model of Sharpe and Singleton, arising from spontaneous symmetry breaking and using a low energy effective Lagrangian which also describes lattice artifacts. The agreement with the Sharpe-Singleton model is remarkable because in the continuum limit in this model the phase structure of lattice QCD with Wilson quarks approaches fast -at a rate O(a 2 ) -the expected phase structure of QCD near zero quark mass. This is an important property which has to be required from any lattice regularization of QCD.
It should be clear that our work can only represent a first step in a detailed understanding of the QCD phase diagram at zero temperature near vanishing quark masses. Clearly, substantially more work has to be done to resolve this phase structure and its behaviour in the continuum limit. For instance, at present for µ = 0 we are unable to differentiate between a scenario where the first order phase transition persists and another one where at µ = 0 only a remnant of the phase transition at µ = 0 is seen. In this respect an analysis like in ref. [23] for µ = 0 is very helpful [28].
Among the many open questions there are: How fast does the gap vanish when the continuum limit at higher values of β is approached? How are the signs of metastability related to the ones observed using the Wilson plaquette action and clover-improved Wilson fermions? How precisely do the eigenvalues re-arrange when the critical quark mass is crossed? Do different gauge actions change the couplings of the effective potential and may hence lead to avoid the phenomena of metastability and reproduce the ideal phase structure at vanishing quark mass already for non-zero lattice spacing?
The most important question is, of course, how phenomenology can be done given the metastability phenomenon seen in our present results; i.e. what is the lowest value of the quark mass that can be reached before one enters the regime of metastabilities and how does this change with decreasing value of the lattice spacing.

A Appendix
A.1 Even-odd preconditioning for the HMC algorithm where the matrix M can be easily read from eq. (7). Using M one can define the hermitian operator where the submatrices Q ± can be factorised as follows: and we have definedμ ≡ 2κµ. Note that (M ± ee ) −1 can be easily computed to Using det(Q) = det(Q + ) det(Q − ) one can now derive the following relation (an equation apart from an irrelevant factor): whereQ ± is only defined on the odd sites of the lattice. In the HMC algorithm the determinant is stochastically estimated using pseudo-fermion fields φ o : where the fields φ o are defined only on the odd sites of the lattice. In order to compute the force corresponding to the effective action S b we need the variation of S b with respect to the gauge fields (using δ(A −1 ) = −A −1 δAA −1 ): with X o and Y o defined on the odd sites as whereQ † ± =Q ∓ has been used. The variation ofQ ± reads and one finds where X, Y is now defined over the full lattice as In addition, δQ + = δQ − , M † eo = γ 5 M oe γ 5 and M † oe = γ 5 M eo γ 5 has been used. Since the bosonic part is quadratic in the φ o fields, the φ o are generated at the beginning of each molecular dynamics trajectory with where R is a random spinor field taken from a Gaussian distribution with norm one.

A.1.1 Hasenbusch trick
The trick first presented in [10] is based on the observation that writing is advantageous for the HMC, if the condition number of W + W − and of (Q + Q − )/(W + W − ) is significantly reduced compared to the condition number of only (Q + Q − ). In order to achieve this we define Withμ 2 =μ + ∆μ it follows immediately that the condition number of W + W − is lower than the one of Q + Q − if for λ min and λ max the lowest and the largest eigenvalue of Q + Q − , respectively, |λ min | ≪μ 2 2 ≪ |λ max | holds: the condition number of W + W − is |λ max |/μ 2 2 while the one of (W + W − ) −1 (Q + Q − ) 2 contrariwise isμ 2 2 /|λ min |. We can takeμ which is a lower bound for |λ min | to write down the following estimates for the condition numbers k: which leads to an optimal choice forμ 2 2 = |λ max | ·μ 2 . As has been shown in [11] also the force contribution coming from (Q + Q − )/(W + W − ) is reduced. This is true also for tmQCD and can bee seen in the following way: noticing that Q + Q − = Q 2 +μ 2 and it follows that Since the corresponding effective action reads one can see that one gets an explicit factor (μ 2 2 −μ 2 ) ≪ 1 multiplying the force contribution compared to the original effective action which will reduce the force and therefore lead to a smoother evolution of the algorithm.
Let us remark that the procedure explained above can be immediately applied to the even-odd preconditioned system. Furthermore the trick can be iterated to two or even more additional operators. In fig. 7 the cost C in units of CG iterations and the acceptance rate P A is plotted versus µ 2 = 2κµ 2 at fixed HMC stepsize and trajectory length. One can see that as expected the acceptance rate increases by introducing an additional operator and reaches a maximum aroundμ = 0.2. Of course also the costs increase when compared to the HMC without additional operators. But the costs are still much less than what is needed to reach an acceptance rate of about 90% without the additional operator (see the dashed line in fig. 7). One can see that the gain in the costs is about a factor of two.
Denoting the projector on the odd subspace by P o we finally obtain the multi-boson representation