Testing the mechanism of lepton compositness

Strict gauge invariance requires that physical left-handed leptons are actually bound states of the elementary left-handed lepton doublet and the Higgs field within the standard model. That they nonetheless behave almost like pure elementary particles is explained by the Fr\"ohlich-Morchio-Strocchi mechanism. Using lattice gauge theory, we test and confirm this mechanism for fermions. Though, due to the current inaccessibility of non-Abelian gauged Weyl fermions on the lattice, a model which contains vectorial leptons but which obeys all other relevant symmetries has been simulated.

calculable, and have been confirmed again in lattice simulations [5,6,16]. It is therefore, at least in principle, possible to observe them in experiment.
Unfortunately, experiments can neither directly employ electroweak bosons as initial states nor detect them as final states, and thus these slight differences are obscured through production and decay processes. It is thus more natural to investigate the impact on the initial states in collider experiments, i.e., leptons [4,7] and protons [4,17]. The lepton case [2][3][4]15] will be discussed in more detail in section II.
While a direct analytical approach may be possible for leptons [4], protons can likely only be reverse-engineered [17]. In either cases a confirmation of the FMS mechanism in the fermion sector using lattice methods would be a valuable first step, and could provide at least a qualitative insight into expected effects. Unfortunately, this is currently not possible for two reasons. On the one hand, the scales are too separated to be accessible with currently available computer resources. While this could still be controlled using extrapolations, the second problem cannot be circumvented. The lattice regularization introduces a gauge anomaly in the weak interactions, as it is incompatible with parity violation [18]. This problem is unresolved despite many efforts [19][20][21][22][23].
Though a quantitative test is therefore not possible, a qualitative test of the FMS mechanism for fermions is: Parity violation is not an important part of it, and it works in the same way for vectorial fermions. It is thus possible to test the very same mechanism which requires the physical electron to be actually a bound state of an elementary electron and a Higgs field using a vectorial electron. Here, we will perform a first such test.
We do so by investigating a system containing an SU(2) Yang-Mills theory coupled to a gauged scalar doublet mimicking the weak-Higgs subsector of the standard model as well as one generation of vectorial leptons. The latter includes one flavor which is gauged under the weak interaction and two fermion flavors that are ungauged. This allows us to construct a gauge-invariant Yukawa sector which obeys the same pattern as its standard model counterpart. This model, and how the FMS mechanism operates in it, will be introduced in section III. Its lattice version, the so called Wilson-Yukawa model [24][25][26][27] follows in section IV. In this first test, the simulations will be quenched. This is justified as the dynamics of the fermions will not alter the basic principles of the FMS mechanism.
Our primary aims are twofold. On the one hand we want to determine the physical, i. e. gauge-invariant, spectrum of the theory. Thereby, we want to show that there is a genuine composite bound state of the elementary, gauged fermion and the Higgs field. On the other hand we want to demonstrate that the FMS mechanism works and predicts correctly the mass of this state. The corresponding lattice observables will be given in section V. The results for both items will then be presented in section VI.
We indeed find hints for a positive confirmation of both items and summarize our results in section VII. There, we will also discuss potential next steps as well as implications for experiments.

II. LEPTONS IN THE STANDARD MODEL
In the following, we rehearse the construction of observables for leptons in the standard model [2][3][4]15]. This is particularly useful, as this provides a perspective, especially on the symmetries and degrees of freedom, which is non-standard compared to usual treatments [28]. For simplicity and to be as close as possible to our toy model in section III, we will treat only the weak interaction, the Higgs doublet, and a single generation of leptons. A generalization to the full standard model can be found in [2,3,15]. The other parts of the standard model do not have any relevant influence on the mechanism in the following.
Thus for our purposes, the relevant part of the standard-model Lagrangian is given by Herein W a µν is the usual field-strength tensor of the weak gauge bosons W and Z, and D µ is the covariant derivative in the fundamental representation. The matrixvalued field X contains the components of the usual scalar doublet φ as and thus the standard Higgs fluctuation mode as well as the three would-be Goldstone bosons. The single lefthanded Weyl spinor ψ L is gauged under the weak interaction in the fundamental representation, ψ L = (ν L e L ) T . The two flavors of right-handed Weyl spinors e R and 1 1 We do not consider Majorana neutrinos for simplicity. ν R are not gauged and combined into a flavor doublet If the Yukawa couplings y f vanish, the theory obeys three important symmetries. 2 First, we have the local weak gauge symmetry SU(2) w . 3 Second, we have a global SU(2) Rf flavor symmetry of the right-handed Weyl fermions for our particular case. At this point we would like to emphasize that the components of the gauged lefthanded spinor cannot be identified with any flavor structure as they merely distinguish different gauge charges similar to the color charge in QCD. Finally, we have a less obvious global SU(2) c symmetry which acts only on the scalar doublet as a right-multiplication on X and leaves all other fields unchanged. Basically, this symmetry relates the scalar doublet and its charge conjugated counterpart ǫ ij φ * j (first column of X) in a nonlinear way. The advantage of the X notation is the linear realization of SU(2) c within the standard-model Higgs sector.
If degenerate Yukawa couplings are switched on (within one generation), the global SU(2) c symmetry of the scalar field and the SU(2) Rf flavor symmetry of the ungauged fermions are broken to a diagonal flavor subgroup SU(2) df , which elements d act as X → Xd and Ignoring for a moment the BEH effect, there are four physical fermionic states in the theory which are grouped into two chiral doublets. The first two states are the flavor doublet of right-handed Weyl fermions χ R . One of them is the right-handed charged lepton, χ R 2 = e R and the other the right-handed neutrino χ R 1 = ν R . The other physical doublet is a gauge-invariant, left-handed Weyl bound state, Ψ L = X † ψ L , which is a singlet with respect to the non-Abelian gauge group but carries a global SU(2) c charge. The two components of this doublet will be identified with the left-handed electron and the left-handed neutrino below. In case non-zero Yukawa couplings break the global SU(2) c symmetry and the flavor symmetry SU(2) Rf to the diagonal subgroup SU(2) df , this bound state transforms in the same way as the righthanded fermions. In this way it appears as if in the physical spectrum the diagonal subgroup acts as an effective flavor symmetry for both the left-handed and righthanded sector. Note that the two gauge-dependent components of the elementary left-handed Weyl fermion ψ L do not transform under SU(2) df , but can be transformed into each other via a gauge transformation and can therefore not be associated with physically observable particles.
Likewise, in the bosonic sector the gauge-invariant bound states trX † X and trτ a X † D µ X form the physical Higgs and the physical W and Z bosons, a singlet scalar and a triplet vector with respect to SU(2) c , see [2,3,15] for details of this sector.
With the BEH effect switched on, and fixing to 't Hooft gauge, the Higgs field can be split into its vacuum fluctuations η and its vacuum expectation value. Using the gauge freedom, we conventionally chose the vacuum expectation value to be in the real 2 direction. Thus, we have φ i = v √ 2 δ i2 . At tree-level, this yields the customary result that the gauged and ungauged Weyl spinors can be combined into two Dirac spinors, each with a mass given by m f = y f v/ √ 2, forming the usual leptons [28]. However, these objects are thus gauge-dependent.
Nonetheless, any physical state has to be unaltered by gauge-fixing. Thus only the left-handed bound state Ψ L and the right-handed χ R remain in the fermionic sector. The decisive step is now to realize the FMS mechanism [2,3] to make contact with the conventional and successful perturbative treatment: Expand any gaugeinvariant composite operator in the Higgs vacuum expectation value of the scalar field. This yields to leading order the results shown in table I. As an example, consider the physical left-handed fermion [2,3], where the matrix-valued η contains the usual fluctuation field identified with the elementary Higgs boson and the Goldstone fields in the same manner as X contains φ in (2). To leading order Ψ L thus reduces to the elementary left-handed fermions. The other physical states in table I follow in the same way. Of course, only the total sum in (3) is gauge-invariant, and the leading order alone is not. When now forming a propagator it follows Thus, to all orders in perturbation theory and to leading order in η the propagator of the physical, composite fermion state is given by the gauge-dependent elementary ones, i.e., the propagators of the left-handed charged lepton and neutrino. Especially, the poles and thus the masses coincide. This was shown to all orders in perturbation theory for the Higgs bound-state-elementarystate duality and generalizes straightforwardly to all other standard model particles [5]. In this way, the gauge-invariant Dirac spinor (Ψ L f χ R f ) T describes the physical neutrinos (f = 1) and charged leptons (f = 2) with the same properties as the usual gauge-dependent ones of perturbation theory [2,3,15]. This can be extended to include the hypercharge sector as well as to quarks [4,15].
As long as the correction O(η) is small, this is an excellent approximation. Indeed, lattice results show this to be the case in the bosonic sector [7,8], though the subleading part is not zero [5,16], and could in principle be accessible even in experiments [4,16,17]. That this is the case is a peculiarity of some theories like the standard model [15,29]. In generic theories, qualitative differences may already appear at leading order [10,11], as confirmed by lattice simulations [12][13][14].
The important bottom line is that the physical lefthanded leptons in the standard model are actually bound states of the elementary ones and the Higgs field [2,3]. That they seem to have the properties of the elementary ones is due to the FMS mechanism, as described in equations (3)(4). This could have far-reaching consequences for experiments, if confirmed [4]. The aim here is therefore to test the underlying mechanism. As noted in the introduction, this is not (yet) possible for the standard model (1). Thus, the next step is to create a theory which works in the same way regarding the FMS mechanism, but is accessible to lattice simulations.

A. The theory
The chiral nature of the weak gauge theory is the main problem to directly test the FMS mechanism for leptons via lattice simulations. In order to circumvent this technical problem, we investigate a toy model that replaces the Weyl fermions by Dirac spinors. At the same time, our model should be as close as possible to the gauged Higgs-Yukawa structure of the standard model, e.g., via imposing similar internal symmetries. Thus, we investigate a standard-model-like theory with vectorial leptons instead of chiral ones. This may be viewed as an extension where a further generation with opposite helicities is effectively added to one of the standard-model generations. See [26,27] for earlier discussions of this approach.
More precisely, we study a theory that comprises an SU(2) w gauge theory coupled to a scalar field and a vectorial fermion ψ in the fundamental representation. Further, we have two flavors of vectorial fermions χ which are singlets with respect to the gauge group. The latter are coupled to the gauged scalar and fermion field via Yukawa interaction terms. This particular setup can be used straightforwardly in lattice simulations. In particular, the gauged fermion ψ mimics the left-handed leptons of the standard model, while the two components of χ can be associated with the right-handed electron and neutrino. We will therefore refer to them as vectorial leptons in the following, or just leptons in case it is clear whether the vectorial or standard-model leptons are meant.

Name
Spin SU(2)c SU(2) Rf Operator LO FMS expansion The Lagrangian of this theory is given by It is thus structurally similar to the reduced standard model (1), except that now tree-level masses for the fermions are allowed for the different fermion species. This theory obeys the same symmetries which we discussed for the standard model in Sec. II, i.e., an SU(2) w gauge symmetry, a global SU(2) c symmetry of the scalar field if y f = 0, as well as an SU(2) f flavor symmetry if m χ1 = m χ2 and y f = 0. Except for the additional possibility to break the flavor symmetry by setting m χ1 = m χ2 , the explicit symmetry breaking patterns are the same as in the standard model if we allow for nonvanishing Yukawa couplings or a BEH mechanism via gauge fixing. In the limit of m ψ = m χ f = 0 an additional discrete chiral symmetry emerges, with γ 5 = iγ 0 γ 1 γ 2 γ 3 . As the tree-level masses do not interfere with the FMS mechanism and substantially reduce computing efforts in the lattice simulations, we will consider only finite tree-level masses, and thus this symmetry will be explicitly broken although such terms are forbidden within the standard model. The physical spectrum contains once more the bound state rather than the gauge-dependent ψ. It is again a doublet under the global SU(2) c , and what has been discussed in the previous section for the standard model on diagonal flavor symmetry for y f = 0 applies here as well. Thus, the theory effectively contains the two vectorial neutrino-like operators Ψ 1 and χ 1 as well as the two vectorial electronlike operators Ψ 2 and χ 2 in the fermionic sector.

Tree-level
Switching now the BEH effect on leads to a somewhat different behavior as in the standard model. Rather than to directly obtain two Dirac particles with separate masses, a non-diagonal mass matrix arises for four Dirac fermions due to the aforementioned doubling of degrees of freedom. Introducing a vector (ψ χ 1 χ 2 ) T , the mass matrix reads Solving the associated eigenvalue problem leads to four different mass values With the five available parameters of the fermion sector, it is possible to form all desired mass values. That these masses do not correspond to the flavors of χ and ψ is due to the mixing of ψ 1 and χ 1 as well as ψ 2 and χ 2 caused by the BEH mechanism and the Yukawa couplings. As usual, the corresponding mass eigenstates can be obtained by suitable rotations in field space where we have two different mixing angels θ f for the two flavors/components of χ f and ψ f .
where ζ ± f have mass M ± f . For the present purpose, we select a particular case to reduce the dimensionality of the phase diagram. We set m χ f = m ψ = m and y f = y, and thus reduce the parameter space to two. Then the effective flavor symmetry SU(2) df is unbroken, and the fermion fields arrange in two doublets with masses In this particular case, the mass eigenstates are obtained from the charge eigenstates via 2. Leading-order quantum corrections However, beyond tree-level, the situation changes slightly. While the relations m χ1 = m χ2 ≡ m χ and y 1 = y 2 are protected by the diagonal flavor symmetry, the relation m ψ = m χ is not. As the gauged and ungauged fermion flavors couple in different ways to the weak gauge bosons, this leads to a splitting of both mass terms once quantum fluctuations are considered. At the one-loop level, we have where c y and c W are dimensionless constants resulting from one-loop integrals with an internal fermion line as well as an internal scalar or gauge boson line, respectively. Further, α W = g 2 4π with g the weak gauge coupling.
Including these one-loop corrections for the mass terms, we obtain for the eigenvalues of the fermion mass matrix Here, we have neglected one-loop corrections to the Yukawa coupling. These are stronger suppressed in the weak coupling regime as they are ∼yα W and ∼y 3 . As a consequence, the mixing at leading-order (12) will also change. The mixing angle θ that translates the doublets in the weak charge eigenstates into the mass eigenstates, ζ + = ψ cos θ + χ sin θ and ζ − = χ cos θ − ψ sin θ reads at one-loop order, Thus, the maximal mixing of ψ and χ in the degenerated case m ψ = m χ will be altered. In particular, θ becomes small if either y is small or α W is large causing a larger split m ψ − m χ .

C. FMS prediction
The FMS mechanism can be applied in this theory in the same way as in the standard model discussed in section II. In the bosonic sector this leads to the same results. It gets more interesting for the hybrid Ψ. In analogy to (3), the FMS mechanism yields Therefore, the Ψ bound-state can be mapped on the gauge-dependent elementary fermion ψ. This implies that Ψ is not a mass eigenstate of the gauge-invariant spectrum as ψ is a linear superposition of ζ ± and we expect two poles at M ± for ΨΨ . Of course, gaugeinvariant mass eigenstates can be constructed from Ψ and χ via a suitable rotation in field space as in the elementary case. For this purpose, a full analysis of the correlation matrix is required which has to include crosscorrelators of the form Ψχ and χΨ and their corresponding FMS expansions. Treating all FMS expanded terms perturbatively, the FMS mechanism predicts that the mixing of χ and Ψ is given by the mixing of χ and ψ. This can be most easily seen by choosing on-shell renormalization conditions for the n-point functions containing elementary fields and composite operator insertions following the lines of Ref [5]. In this case it can be shown that the poles and their residues of the gauge-invariant correlators coincide with their gauge-dependent counter parts. Of course, this is a perturbative statement and nonperturbative bound state effects might alter this behavior. However, we do not expect strong deviations in the weak coupling regime.

A. Dirac operator
In order to discuss the discretization 4 of the theory which we described in section III A we rewrite the fermionic part of the action as [30] with and thus in a block-diagonal form in which the interaction with the Higgs field through the Yukawa interaction explicitly appears in the off-diagonal parts.
To obtain a lattice version, the standard discretization of the bosonic sector is used [30,33]. For the fermionic action, we give for future reference the most general version, and only specialize afterwards to our case of degenerate parameters. The first diagonal block has then been implemented as a standard SU(2) Wilson-Dirac operator [34] With this notation, we have also the relation between the parameters κ ψ = 1 2(m ψ +4) for the gauged fermion. For the γ matrices we used standard chiral Euclidean ones [34].
The second diagonal block has been implemented as a free Wilson-Dirac operator where κ χ f = 1 2(mχ f +4) are the two hopping parameters for the ungauged fermions.
The third and the fourth block are the Yukawa couplings between the fermions and the Higgs. Due to the Euclidean spacetime they obtain an overall sign factor, but remain otherwise close to the continuum form The combined lattice operator is called Wilson-Yukawa operator [24][25][26][27].
In the following, we will set κ F = κ ψ = κ χ1 = κ χ2 and Y = Y 1 = Y 2 . We will furthermore quench the fermionic sector, as will be discussed in more detail in section IV C. Note that because of the rescaling of the Higgs field and the fermion fields with their own hopping parameters the lattice Yukawa couplings obey The inverter used in our work is based on the BiCGstab method as discussed in [35]. In particular, our implementation is based on the application of the full operator on a vector v(x) α,i which acts as a fermionic source. The index α is a Dirac index, while the index i indicates the fermionic species, and runs over both the gauge components of ψ and the flavor components of χ, and thus over the four-dimensional Dirac operator (16). The implementation has been checked by calculating the trace of the full operator on a small volume, in the free case with a static Higgs field, and by comparing it with the result from an algebraic computation software. Parallelization has also been enabled for the algorithm, with the openMP API, and tested with respect to the serial results.
For our spectroscopical purposes, we used a single point source located in the origin. This strategy required 16 inversions, given the 4 Dirac indices and the 4 different fermionic species included in the multifield. This point source gave sufficient statistical accuracy for our purposes.

C. Phase diagram and simulation points
The computational costs for the full theory are, as generically for theories with dynamical fermions, very high. This can already be gathered from simulations without the ungauged fermions in the QCD-like domain [36]. However, as in the standard model, we do not expect a substantial influence of the fermions on the FMS mechanism regarding the mass spectrum of the theory. Thus, we will investigate a quenched scenario in the following for a first qualitative check.
We thus calculate the fermionic spectrum on configurations with dynamical gauge and Higgs fields created using the methods of [8,33]. This includes a subset of configurations fixed to minimal Landau-'t Hooft gauge, for which the algorithms of [8,12] were used. These were necessary to determine the propagator of the gaugedependent field ψ. Note that this propagator has much less statistical fluctuations, and we therefore needed only O(50) configurations for it, while using O(1000) configurations for the gauge-invariant quantities. As gaugefixing is very expensive, this leads to roughly the same total computing costs for both types of configurations. However, ultimately the computing time was dominated by the calculation of the fermionic observables below. We list the parameter sets in table II. They were selected for being suitably similar to the standard-model case in the bosonic sector at either weak coupling or somewhat stronger coupling at different discretizations. However, as will be seen, all parameter sets show essentially the same behavior.
The fermionic sector of the theory described in section III has still 5 parameters, three hopping parameters κ ψ , is the mass of the scalar bound state tr(X † X) corresponding to he Higgs boson of our toy model. The running coupling, and thus the vacuum expectation value, are in the miniMOM scheme [8,37] evaluated at 200 GeV. The results are from the largest volumes employed here, 24 4 . κ χ f , and the two Yukawa couplings Y 1 and Y 2 . As noted these are reduced to two by κ F = κ ψ = κ χ f and Y = Y 1 = Y 2 , leaving only κ F and Y . It remains to find values for these parameters such that the physics is the one expected for (heavy vectorial) leptons, while at the same time being accessible with available computational resources.
For this purpose we investigated a wide set of κ F and Y values. We found that the time needed for inversion increased substantially for κ F 1/8, and thus for negative tree-level masses. At the same time, for 0 < κ F ≪ 1/8 the observed masses in the fermionic sector became very heavy, above one in lattice units. Thus, as a compromise we selected κ F = 0.11 and κ F = 0.12 as simulation points. For the Yukawa couplings we choose Y = 0.01, Y = 0.05, and Y = 0.1. Again, for larger Yukawa couplings the inversion time became very long. This is to be expected, as the lower mass, according to either (11) or (14), then approaches zero, which yields high inversion times of the Dirac operator. At smaller Yukawa couplings their impact on the masses became too weak to be detectable within our precision. Note that the theory is symmetric under a change of sign of the Yukawa couplings, and thus we can keep with positive values.
To keep finite-volume effects under control, we did simulations for 5 different lattice volumes, 8 4 , 12 4 , 16 4 , 20 4 , and 24 4 . However, we find that even for the finest lattices the infinite-volume behavior has been reached, within available statistical uncertainty, essentially at 20 4 . This is discussed in appendix A in detail, alongside other lattice artifacts. Thus, these choices are sufficient for our purposes. Hence, in total 150 different sets of lattice parameters have been investigated, all in all a little more than 50.000 configurations. The dominant statistical uncertainty stems from the hybrid Ψ bound state, probably due to the strongly fluctuating [8,33] Higgs field component.

V. SPECTROSCOPIC OBSERVABLES
Regarding spectroscopy we are interested in the two point functions, which can be accessed from the inverted Wilson-Yukawa operator. We are interested in the bound state Ψ, the gauge-invariant fermion χ, and the gaugedependent fermion ψ. Strictly speaking Ψ and χ have the same quantum numbers, and thus we are looking in principle for the ground state and the first non-trivial excited state. As the excited state could potentially decay, this would require in principle a Lüscher-type analysis. However, as we will see, our results show agreement with the analytical investigation of section III, implying that the corresponding states are almost stable. This is confirmed in a few cases, where indeed decays and scattering states are kinematically forbidden, and thus we have two stable fermionic states in the physical spectrum. Therefore, a straightforward discretization of the desired propagators of ψ, χ, and Ψ turns out to be sufficient for the purposes at hand.
In the gauge-invariant sector, the full cross-correlation matrix for the two states Ψ and χ for the propagator is obtained from the Dirac operator (16) and the bound state structure (7) in a straightforward way by Wick contraction [34], Here we used the fermionic species as subscripts to distinguish the various sections of the inverted operator. They should not be confused with the elements of the Dirac operator in (16). All of the elements of the propagator matrix (20) have a full Dirac matrix structure. For spectroscopy we use the trace in the Dirac structure. Thus, we get in the off-diagonal blocks the possibility for mixing between the bound states and the elementary fermions, just as with left-handed and right-handed particles in the standard model [4]. Albeit a full variational analysis of (20) reveals that substantially more statistics is required for quantitative precision, we will already get important insights from it. Thus we use the diagonal two propagators individually as well as the eigenvalues from the variational analysis of the full matrix below.
In addition to the gauge-invariant observables we also consider the gauge-dependent ψ. In order to obtain a non-zero result for it, it is necessary to invert the operator only on gauge and scalar configurations in a fixed gauge. Conceptually, we have to investigate the cross-correlation matrix for the elementary fermion fields, i.e., and compare the results to Eq. (20) to test the FMS mechanism. Of course, matrix (21) is only meaningful within a gauge-fixed setting as otherwise all elements except (D −1 )χ χ vanish. Once more, a full variational anal-ysis is expensive and requires more statistics as we have currently available. Thus, we restrict ourselves again on the diagonal elements as they will contain for our purpose all relevant information about the spectrum. As the pure χ-field propagator χχ is gauge-invariant, it does not matter if we calculate it on a gauge-fixed or nongauge-fixed configuration. Therefore, we use the results of the unfixed configurations for this element and focus in the gauge fixed set up on the Dirac trace of the (D −1 )ψ ψ component.
Note that on our gauge-fixed configurations, where the gauge symmetry and the global symmetry are broken to a common subgroup, it is possible to also define an extended propagator matrix The additional off-diagonal elements describe the gaugedependent overlap of the elementary ψ field with its gauge-invariant bound state counterpart Ψ. In principle, one could perform a variational analysis on the matrix defined in Eq. (22) to obtain a comprehensive understanding of the fermion sector of this model. However, we are here interested in a first qualitative comparison of the FMS mechanism for fermions. Hence, we will separately analyze (20) and the (D −1 )ψ ψ component of (21) which is sufficient for our task. As usual, we will perform a zero-momentum projection, which is executed on the inverted matrix elements We define the effective masses using the quantity In this notation, we indicate both the diagonal elements of M , so that we have but also the masses obtained from the eigenvalues of M , which are also called the principal correlators The correlators M (t) and λ i (t) show a sum-of-coshs behavior in our finite volumes. As will be seen the correlators contain multiple levels, and thus fits require multiple cosh terms.

A. Impact of quenching on the FMS predictions
For the spectroscopic results, we have in principle clear predictions from the analytical investigations in sections III B 2 and III C. Nonetheless, a few more comments are in order. On the one hand, lattice simulations will definitely capture more information about the system than the basic one-loop approximations done in Sec. III. On the other hand, the analytical predictions are derived for fully dynamical fermions while we use a quenched scenario for the simulations in the following to gain a first qualitative investigation of the spectrum. For an actual comparison, we have to either do the FMS analysis within a quenched setting or to perform unquenched simulations. The latter is clearly beyond the scope of this work due to the high computational costs. The former option is challenging as well because a direct translation of the quenched approximation into a continuum formulation is involved.
As a first heuristic step into this direction, we redo our analytic calculation by assuming that the mass of the fermions is much larger than any momentum scale for internal fermion lines. This causes an effective suppression of fermion fluctuations that will properly reproduce the quenching effects in the bosonic sector of the model. However, the situation is less clear for fermionic observables. We find that the mixing effect which exists for dynamical fermions is not captured by the quenched approximation. This manifests in the pole structure of the propagators χχ and ψψ . For the dynamical case, we have two distinct poles which are present in both propagators. Within our handwaving modeling of the quenched approximation we find that each propagator contains only one pole given by the (quenched) quantum corrected versions of m ψ and m χ for ψψ and χχ , respectively. In the following, we will denote these infrared mass terms by M ψ and M χ to avoid confusion with the bare mass parameters. We will obtain similar results from the lattice investigations in a moment. Thus, our analysis seems to be consistent from that perspective.
However, at this point we would like to mention that our modeling of the quenching has some ambiguities when we resum the propagators. For instance, we obtain not only simple pole terms ∼1/(p 2 − m 2 ) but also terms of the form ∼1/(p 2 − m 2 ) 2 as we treat internal and external fermion lines differently. Of course, more sophisticated methods are available to model the quenching, e.g., via introducing additional ghost fields that precisely cancel the fermion determinant which was used in the context of quenched chiral perturbation theory [38][39][40][41][42]. However, the influence of these ghost fields will manifest only within higher loop terms in the fermion propagators when applied to our model and also suffer from unitarity issues.

B. Lattice results
In the following, we will go carefully through each of our lattice findings. We will exclusively use the infinitevolume extrapolated results, which we obtain along the lines described in appendix A. However, the results on the 20 4 lattices are already compatible with the infinitevolume results within statistical errors, so this is of little actual concern.
The first interesting question is the gauge-invariant ground state. We find the following pattern. Performing a variational analysis of the effective masses from the gauge-invariant sector, i. e. (20), we find that the lowest level is the same as would be obtained directly from investigating the lower diagonal element, i. e. the correlator of the ungauged fermion χ alone. This ungauged fermion correlator shows very little noise, and can be very well captured by a double-cosh fit, as shown for an example in figure 1. We therefore conclude that this operator has (essentially) perfect overlap with the ground-state and is the physical lightest particle in this quantum number channel. Because it is a fermionic state, it cannot decay into any of the bosonic particles, and is absolutely stable. Note that in this, and all other channels, the correlator at t/a 2 is dominated by a mode which has a mass of order π, and thus is a lattice artifact, which will not be considered further.
The next object is the gauge-fixed propagator of ψ, (21). Due to the gauge-dependency [8] we find a (very) slight non-monotonous behavior in the correlator at short times. At later times the effective mass exhibits a good plateau, which is very well described by a single cosh term. This is illustrated in figure 2. We use this plateau  to determine the mass of the gauged fermion. Thus, we conclude that there is only a single state in this channel as well.
The fact that the elementary fields are only dominated by a single mass state is in contradiction to the results of Sec. III but as outlined in Sec. VI A, we trace back this circumstance to the quenching. We find further evidence for this conclusion by the following points. First, we checked the mixed correlators on the diagonal terms of Eq. (21) which describe the overlap of ψ and χ on gauge-fixed configurations. We indeed find results that are compatible with zero at long times, see Fig. 3. 6 Second, we find that the mass M ψ extracted from ψψ is substantially larger than the mass M χ . This also fits into the picture as the mass term for ψ gets additional corrections from gauge boson fluctuations which are not present at one-loop order for M χ , cf. Eq. (13).
The results for the two masses M ψ and M χ are listed in table III. We find that M χ is effectively independent of the gauge coupling and varies only with the Yukawa coupling and indirectly with the parameters of the scalar sector. By contrast, M ψ depends on the gauge coupling and the ratio M ψ /M χ tends to be smaller with smaller gauge coupling which is expected. The masses are described almost always within 1σ statistical error by a form motivated by the Taylor expansion at small y of (13). We provide the fit parameters in table IV. Though in principle this can be translated into the form (13), the values should not be identified with the parameters there, as the present ones are the quenched lattice ones. Further, the fits have been performed in the lattice Yukawa coupling. The actual Yukawa couplings need to be rescaled by √ κκ F ≈ 0.06 on average, giving again reasonable numbers.
This leaves the bound state Ψ. This turns out to be quite complicated, especially as it is substantially more 6 Such a demixing is expected if a random ungauged flavor transformation is applied. Because in the quenched case the two global symmetries are independent at the level of the path integral, and not locked due to the interaction, this is a possible explanation. Conversely, the masses are not affected, and are thus faithfully reproduced.
# κF Y aMχ aM ψ r 1 0.11 0.01 0.421 +0.001 −0.008 0.817(3) 5.7 1 0.11 0.05 0.407 (6) 0.77(3) 0.5 1 0.11 0.1 0.353 (9) 0.54 (1) (27)(28). Note that the fit is done in the lattice Yukawa coupling Y and not the continuum y. noisy than the other two states. However, for the smallest and largest tree-level Yukawa couplings the correlator shows a mass, which is essentially the one of the ψ channel and the χ channel, respectively. This is also confirmed by the variational analysis of (20), though the errors are considerably larger. This is illustrated in figure 4. The situation is different for the intermediate Yukawa coupling. A naive analysis yields a mass in between the two other channels, with relatively larger errors. However, close scrutiny actually yields that there is a systematic trend of the correlator. The puzzling situation can be resolved by assuming that the bound state correlator is determined by a combination of the two elementary propagators. Thus, a single-parameter fit of type with r to be fitted, is performed for late times. This yields a much better agreement. This is also supported by the variational analysis, which indeed finds as second eigenvalue only the mass value M ψ in the same channel, albeit at substantially larger errors. Reiterating the same procedure also for the other two Yukawa couplings yields that also in that case the correlator is well-described by (29), though extremely strongly dominated by either of the two terms. 7 How this is realized is shown for an example setup in figure 5, and the values for r are also listed in table III. Note that for many systems the mass difference of both states are less than the mass of the scalar singlet or twice the mass of the vector triplet, and thus the heavier state cannot decay into the ground state, and it is thus stable.
There is a subtlety with determining r. The effective mass for the bound state is substantially more noisy than for the elementary fermions. However, (29) limits the upward fluctuations of the effective mass to the error of the larger mass M ψ and the downward fluctuations likewise to the error of the smaller mass M χ . Thus, any attempt to find an error band for r would require to allow for relative errors of the input masses in (29) that are larger, and of the same relative size as the one from the bound state correlator itself, and thus than actually are observed and listed in table II. This would be artificial. The situation is further complicated because of the closeness of the masses M ψ/χ . We therefore quote only an optimal value for r by varying it such that the fit (29) goes best through the effective mass error range of the bound state, as shown in figure 5.
That our results are indeed compatible with the fact that the bound state operator Ψ is a mixture of two mass states that are described by ψ and χ is consistent with the FMS picture although the quenched analysis makes the interpretation less obvious. As discussed in Sec. III C, the FMS mechanism projects the on-shell properties of Ψ onto the the on-shell properties of ψ for a model with dynamical fermions. Thus, Ψ has overlap with two mass eigenstates and the mixing angle of Ψ and χ is the same as the mixing angle between ψ and χ at least in a perturbative set up. As we have accumulated evidence that the quenched calculation does not resolve the mixing between the elementary states, one would naively conjecture that the bound state is only described by M ψ . However, a quenched (continuum) FMS analysis would actually be necessary to make a decisive statement. As a quenched continuum analysis is already involved for the elementary fields, a detailed analysis for the FMS mechanism is beyond the scope of this work. Here, we conjecture that the FMS mapping of the mixing of states gets altered as we expect a nontrivial modification of the higherorder FMS terms due to the quenching. In particular the higher-order terms allow for intermediate states that are described by M χ and M ψ causing the observed overlap with both states. The relative ratio between both states is given by the strength of the Yukawa coupling which is confirmed by our data. 8 We can therefore conclude, that our results support that the physical spectrum is given in terms of the χ fermion and a bound state of the ψ fermion and the scalar field X. The masses of these states are correctly predicted by the FMS mechanism to be the ones of the elementary fermions, both gauge-invariant and gaugedependent. 8 The lattice also contains nonperturbative information which might not be present within a perturbative treatment of the FMS mechanism. Only an unquenched analysis could reveal as to whether such effects exists and might further modify the mixing.
Eventually, the physical spectra are shown as a function of the fermionic parameters in figure 6. It is consistent with the FMS prediction in all channels, both bosonic [8], and fermionic, for all parameters. It is also visible that a suitable choice of parameters allows to have both very heavy and very light fermions in the spectrum. In fact, by suitably tuning κ F and Y , it appears there would be no obstacle in tuning through the full range of masses from the lightest neutrino to the top quark, and even beyond.

VII. CONCLUSIONS
We have collected evidence that the FMS mechanism is also working in the fermionic sector of gauge-Higgs theories, as anticipated already 40 years ago in [2,3]. Especially, we find that the physical spectrum is indeed a mix of ungauged (would-be right-handed) fermions and (would-be left-handed) fermion-Higgs bound states.
While we have yet investigated a quenched, vectorial system, there is no conceptual difference to the one in the standard model [2][3][4]15], or even beyond where also gauge-invariant fermion-Higgs bound states are expected [15,43,44]. Furthermore, even though this covered only lepton-like states, the qualitative mechanism is the same also for hadrons [4,15,17]. Of course, this is no guarantee that the mechanism works indeed in all these cases. Ultimate proof will require a detailed analysis of these systems. However, there is no obvious reason that it should not, especially given that the FMS mechanism has passed so far all (lattice) tests [7,8,[12][13][14] in various theories.
This results should serve as a field theoretical foundation for a treatment of fermions in theories with BEH effect, which take fully into account the invariance of the observables under the full gauge group. This includes the standard model. This has far-reaching consequences for phenomenology, as this substructure could be accessible at future colliders [4], like the proposed CLIC [45], FCC-ee [46], or ILC [47]. Identifying and measuring this substructure easily forms an experimental program in itself. Exploiting the Higgs component also allows to increase the reach for new physics searches which couple to the Higgs component directly, like dark matter through Higgs portals. The results here motivate strongly an experimental program aimed at these effects. Moreover, as it is uniquely tied to the field-theoretical structure of the standard model, it is a guaranteed discovery: Either this effect is found, or there is something very different at work, probably new physics.
A natural next possible step, aside from the obvious but expensive unquenching and reduction of lattice artifacts, is the study of the structure of the Higgs-fermion bound state. Especially form factors, which already illuminated the substructure and size of the vector bosons [16], are an obvious next goal. This should give a first idea of anomalous couplings to the Z, as well as the weak FIG. 6. Spectra for all setups. The bosonic sector is fixed, and plotted to the left, while the two states identified in the fermion channel, each being a doublet, are plotted for the different parameters to the right.
radius of the fermion-Higgs bound state 9 , both of which are highly interesting questions at future lepton colliders [49]. servables [8,13,16,33], and appears to be rather generic for this type of theory. Though for a detailed quantitative understanding an extended investigation of lines of constant physics will be necessary, at the qualitative level of this work this is sufficient. Of course, due to the fact that Wilson fermions require mass renormalization [34], it would be necessary to change the values of the fermion parameters to keep the physical masses fixed when moving along such lines of constant physics. Even in the present case, where the unbroken global symmetry requires some of the parameters to always coincide, this would be a further additional logistical and computational complication. Fortunately, as for the present purpose it is sufficient to have some values of the masses, this is not necessary. When in a next step a continuum extrapolation will be attempted, this will change.
In this context the use of Wilson fermions, which break chiral symmetry explicitly [34], could be problematic. In particular as the mass generation by the BEH effect is, as in the standard model, dynamic. Hence, similar interference as in QCD may be possible [34]. However, in the present model theory we use an additional explicit breaking to avoid negative tree-level masses according to (11). Since the two, quite different, tree-level masses used in the main text do not show substantial differences in behavior, we conclude that, as in heavy-quark QCD, we are still in a region where we this effect is negligible.
The strongest dependence we see are the dependencies on the volumes. Depending on the lattice parameters, these can be relevant, especially when the comparison of masses of different states is done as needed here. We therefore use a fitting ansatz [50] m N = m ∞ + a N e −bN (A1) with the lattice extension N to determine the infinitevolume mass m ∞ . This is done using the volumes of 8 4 to 20 4 , on which the change is largest. The results from the 24 4 lattices are then used to confirm the fit result. An example for the finest lattice, and thus most extreme volume effects, is shown in figure 7. The results confirm that we see already on the 20 4 lattice infinitevolume behavior for the masses. In the main text only these infinite-volume masses have been used 10 . The value of r seems to be more dependent on the volume. However, it is visible in figure 7 that the value is substantially dependent on whether the finite-volume masses or the ones extrapolated to infinite volume are used in the fit (29). There is unfortunately no expected behavior for it. However, the results slow down quicker than linear in 1/N , and a quadratic extrapolation suggests a finite value above zero. Nonetheless, even if the value would be close to zero, or essentially zero, this would only mean that the oscillation effect is strongly volume-dependent, and eventually wins for the intermediate values of the Yukawa couplings.