Tunable Stripe Order and Weak Superconductivity in the Moir\'e Hubbard Model

The moir\'e Hubbard model describes correlations in certain homobilayer twisted transition metal dichalcogenides. Using exact diagonalization and density matrix renormalization group methods, we find magnetic Mott insulating and metallic phases, which, upon doping exhibit intertwined charge and spin ordering and, in some regimes, pair binding of holes. The phases are highly tunable via an interlayer potential difference. Remarkably, the hole binding energy is found to be highly tunable revealing an experimentally accessible regime where holes become attractive. In this attractive regime, we study the superconducting correlation function and point out the possibility of weak superconductivity.

The moiré Hubbard model describes correlations in certain homobilayer twisted transition metal dichalcogenides. Using exact diagonalization and density matrix renormalization group methods, we find magnetic Mott insulating and metallic phases, which, upon doping exhibit intertwined charge and spin ordering and, in some regimes, pair binding of holes. The phases are highly tunable via an interlayer potential difference. Remarkably, the hole binding energy is found to be highly tunable revealing an experimentally accessible regime where holes become attractive. In this attractive regime, we study the superconducting correlation function and point out the possibility of weak superconductivity.

I. INTRODUCTION
Twisted moiré materials have attracted intense recent attention due to the range of correlated phenomena they exhibit and their versatile experimental tunability. Exotic emergent phenomena including correlated insulating states [1,2], quantum critically and tunable metal-insulator transitions (MIT) [3,4] have been recently realized in the twisted WSe 2 system, a typical class of transition metal dichalcogenides. The low energy physics of twisted WSe 2 is well captured by the so-called moiré Hubbard model, a variant of the standard triangular lattice Hubbard model in which the electron hopping amplitude acquires a spin-dependent phase [5][6][7][8][9][10]. The Hamiltonian H of the moiré Hubbard model is the sum of kinetic (H 0 ) and interaction (H I ) terms with H 0 = −t σ=↑,↓ r,j=1,2,3 e iσφ c † r+aj ,σ c r,σ + h.c. , (1) and H I = U r n r↑ n r↓ is the standard onsite repulsive Hubbard interaction. In Eq. 1, a 1,2 are primitive lattice vectors of the moiré triangular lattice with relative angle 2π/3, a 3 = −a 1 − a 2 shown in Fig. 1, and σ =↑, ↓ represents the spin of the electron.
The new feature of the model is the spin-dependent phase φ in hopping amplitude, te iσφ . Initially obtained from DFT calculations [2], it captures the effect of the displacement field (gate voltage difference) between two layers. The displacement field breaks the inversion symmetry between two layers. The introduced phase is the simplest way of breaking inversion symmetry, which breaks the C 6 rotation symmetry of the original model * awietek@flatironinstitute.org down to C 3 , and the fact that the phase is spin-dependent guarantees the time-reversal symmetry. φ produces a flux of 3φσ per triangular plaquette, alternating in sign between adjacent triangles and opposite for spin up and spin down. Since a flux of 2π per triangular plaquette is equivalent to zero flux, the model is invariant under φ → φ + 2π/3. Further, a flux of π ≡ 3π per triangle corresponds to changing the sign of the hopping along each bond, equivalent to a particle-hole transformation. Thus the spectrum of the model at density n and phase φ is the same as at density 2−n and flux φ±π/3 [5][6][7]9].
Importantly, both the carrier concentration and the phase are greatly tunable experimentally by the voltages associated with each layers. The sum of the two layer voltages determines the chemical potential while their difference, i.e. the displacement field, tunes the phase φ [2,[5][6][7]9]. In tWSe 2 , physically achievable values of the displacement field correspond to changes of φ over the range of − π 3 φ π 3 . In this work, we explore the physics of the moiré Hubbard model along the aforementioned experimental tunable degrees of freedom: doping and displacement field. Through a combined exact diagonalization (ED) and density matrix renormalization group (DMRG) study, we obtain unambiguous numerical results that provide insight into the physics of the moiré Hubbard model. Key results include an approximate phase diagram at half-filling, carrier pairing, and indications of superconductivity away from half-filling at non-zero displacement fields.

II. PHASE DIAGRAM AT HALF-FILLING
To establish the ground state physics of the model Eq. The colors denote different quantum numbers of the first excited state. Red (Blue) colors correspond to Sz = 0 (Sz = 1) states at different momenta and point group symmetry representations, see appendix A. The magnetic insulating regime at large-U extends close to U = 0 for φ = π/6, π/2, 5π/6. Two 120 • Néel ordered states in the x-y plane with opposite chirality are realized as well as a x-y ferromagnet. The flux pattern of the moiré Hubbard model and the high symmetry momenta in the Brillouin zone are shown bottom left and right.
U/t and the flux φ, we employ ED on finite periodic lattices [11]. We focus on the N s = 12 cluster, since this cluster is highly symmetric, can stabilize 120 • degree magnetic orders, and also features the M point in the Brillouin zone, which we find to be important for the low-energy physics of the system. The next larger cluster featuring both the K and M point while having a C 3 symmetry would have N s = 36, which is not within reach of ED.
The upper panel of Fig. 1 shows the phase diagram obtained from ED, along with color-coded indications of the momentum, point group representation, and the total magnetization S z . In symmetry breaking phases, the quantum numbers of the low-lying "tower-of-states", which become the degenerate ground states in the thermodynamic limit, can be computed from group representation theory, Ref. [12].
To complement this characterization of ground state phases we computed the ground state magnetic structure factor (in the x-y plane), as well as the single-particle gap ∆ p , and the spin gap ∆ s , where E 0 (m, n) denotes the ground state energy of the system with a number of m electrons with spin up, and n electrons with spin down. The results as a function of U/t and φ are shown in Fig. 2. At large U/t, we find three different insulating regimes, with three different magnetic orders in the x-y plane. For phase φ ∈ (π/3, 2π/3) the spins orient ferromagnetically. This is evident in the large value of the magnetic structure factor at momentum Γ ≡ (0, 0) in Fig. 2(b). Consistent with this interpretation, the first excited state has discrete momentum, k = Γ = (0, 0) and non-zero magnetization S z = 1, as expected for a translationally invariant state with ferromagnetic order. The insulating magnetic regimes at φ ∈ (0, π/3) and φ ∈ (2π/3, π) exhibit two peaks in the structure factor at k = K 0 ≡ (4π/3, 0) and k = K 1 ≡ (2π/3, 2π/ √ 3). This is indicative of the planar 120 • Néel order, which can have two different chiralities. While the structure factor cannot distinguish between the two chiralities, the momentum of the first excitation indicates which chirality is realized, with the first excitation for φ ∈ (0, π/3) exhibiting excitation wave vector k = K 0 while for φ ∈ (2π/3, π) yields momentum k = K 1 , which distinguishes the two chiralities.
In order to distinguish the metallic from the insulating regime, results for the single-particle gap ∆ p as defined in Eq. (3) are shown in Fig. 2(c). We find that indeed the magnetically ordered regions coincide with regions of an enhanced single-particle gap. The magnetic insulating regimes extend close to U/t = 0 for values of φ = π/6, π/2, 5π/6. These findings are consistent with previous results from a Hartree-Fock approximation [6], where these three magnetic orders at large-U as well as a metallic regime at small-U were found. However, the MIT found here is shifted towards larger values of U/t than the Hartree-Fock transition, consistent with recent dynamical meanfield results [9] and with more elaborate numerical results of the case φ = 0, where the critical interaction strength for the MIT has been estimated to be U c /t ≈ 8.7 [13].
Apart from the magnetic regimes Fig. 1 shows several other regions with different quantum numbers of the first excitations, each of which could indicate a different phase being realized. At φ = 0, we find that in the intermediate coupling regime 9 U/t 11 the lowest excited state has S z = 0, momentum q = Γ and, importantly, a point group representation E of the discrete Magnetic order and gaps of the Ns = 12 cluster from Exact Diagonalization. The magnetic structure factor Sm(q) evaluated at points q = K (a) and Brillouin zone center q = Γ (b), indicating 120 • Néel order and ferromagnetic order in the x-y plane, respectively. We observe a discontinuity as small values of U/t around special values of φ = π/6, π/2, 5π/6. (c) The single-particle gap ∆p indicates the insulating and metallic regimes. (d) An enhanced spin gap ∆s is observed close to the transition from the magnetic insulator to the metallic regime, possibly indicating non-magnetic insulating states.
dihedral point group D 3 . This point group representation has been observed for chiral spin liquids [14]. This observation would be consistent with recent proposals of a chiral spin liquid in the intermediate coupling regime at φ = 0 [15][16][17]. This putative chiral spin liquid regime is suppressed by a finite displacement field, which is consistent with the expectation that breaking the SU(2) symmetry disfavors the formation of spin singlet and hence makes the spinon condensation harder to realize [8].
At the particular values φ = π/6, π/2, 5π/6 close to U/t = 0 we observe several regimes where the first excitation has momentum q = M ≡ (π, π/ √ 3). The representation of the ground state is also different from the surrounding regime. Moreover, we observe a discontinuity in the structure factor in Fig. 2(a,b). Therefore, this region could constitute a separate, likely metallic, phase, whose precise nature is yet to be determined.
At the boundary of the magnetic insulators, we observe a state with momentum q = M and S z = 0, which could indicate a non-magnetic insulator. This would be also supported by the fact that the spin gap shown in Fig. 2(d) is rather sizable in this region and exhibits a discontinuity. However, we cannot fully exclude that this is an artifact of the finite cluster size, which could render this particular state to be the lowest excitation. Especially, we do not observe any particular feature in the magnetic structure factor in (a,b) in this regime.

III. PHYSICS AT SMALL HOLE DOPING
A. Intertwined spin and charge ordering Away from half-filling, charge and spin modulations of considerably larger wavelength than accessible by ED become important. We, therefore, apply the DMRG method to study the ground state properties of the system at small hole-doping on elongated cylindrical geometries. In this study, we focus on the YC3 geometry shown in Fig. 3. The lattice is chosen to have periodic boundary conditions along the short cylinder length and open boundary conditions in the long direction, as conventionally chosen for use in matrix product state techniques. The cylinder is well suited to study the ordered phases since it resolves the momenta K 0 and K 1 . While YC4 and YC5 geometries of the triangular lattice Hubbard model at half-filling have been studied to unravel a chiral spin liquid [15,16], these geometries do not resolve the K points, and would hence introduce unphysical frustration. Similarly, twisted boundary conditions shift the resolved momenta. The YC6 cylinder would be a suitable candidate system. However, DMRG simulations at finite hole-doping in addition to the staggered magnetic field on the YC6 cylinder are currently beyond our reach.
We focus on the case U/t = 8, which is believed to be relevant in twisted WSe 2 [2,4], where the system is metallic for φ = 0 while insulating for φ = π/6 at halffilling. We pick particular values of φ = 0, π/6, π/3, and π/2. The case φ = 0 corresponds to the pure Hubbard model with nearest-neighbor hopping, where previous DMRG studies suggested the possibilities of pairdensity wave [18] or chiral metal [19] physics in this intermediate U regime. At φ = π/6 we are doping the system which according to Fig. 1 exhibits 120 • Néel order in the x-y plane at half-filling. For φ = π/2 the system exhibits x-y ferromagnetic order at half-filling. In contrast, at φ = π/3 the system is metallic at half-filling.
We show ground state properties on the YC3 upon doping the ferromagnetic state at φ = π/2 in Fig. 3. We observe a regular charge density modulation, where two holes form one stripe. The sign of the magnetic correlation switches on the maximum of the respective holedensity. Thus, the system exhibits typical intertwined spin and charge order. Originally, such orders were proposed by Hartree-Fock studies on the square lattice [20][21][22][23] and have as of now been firmly established as the ground states in certain parameter regimes of the Hubbard model [24][25][26][27][28].
Hole-density 1 − ni and spin correlation S x i S x 0 + S y i S y 0 of the ground state on the 36 × 3 YC3 cylinder for U/t = 8 and small hole-doping p ≈ 0.074, corresponding to n h = 8 holes at φ = π/2. The reference site of the spin correlations is marked with the black cross. We observe charge modulations with two holes per stripe. The spin correlation switches the sign at the maxima of the hole density, indicating intertwined spin and charge ordering. are shifted by δ = π √ 3p from the ordering vectors K0 (a) of the 120 • Néel order for φ = π/6 and Γ = (0, 0) in (b) for φ = π/2. The charge structure factor Sc(k) in (c) is identical for φ = π/6 and φ = π/2 and is peaked at a wave vector 2δ, indicating stripe order. For φ = π/3 in (d) no peak indicating stripe order is observed. The behavior Sc(k) ≈ α|kx| at small |kx| indicates a metallic state.
To quantify these observations, we computed the magnetic structure factor S m (k) and the charge structure factor, on the 72×3 YC3 cylinder at hole-doping p = 1/18 (n h = 12) for different values of φ. Fig. 4(a) shows the magnetic structure factor S m (k) for φ = π/6. We observed that its peak is shifted from the ordering vector K 0 of the 120 • Néel order by δ = π √ 3p = π √ 3(1 − n). Similarly, we observe a peak in the charge structure factor S c (k) at a wave vector of 2δ in Fig. 4(c). Hence, stripe order where the wave length of the charge modulations is half the wave length of the spin modulation is also realized for φ = π/6 and the spin modulation is a modulation of the 120 • order. To further verify the case of stripe ordering for φ = π/2 shown in Fig. 3, S m (k) shown in Fig. 4(b) is similarly peaked at a small shifted wave vector δ instead of Γ = (0, 0), which would indicate uniform ferromagnetism. The charge structure factor Fig. 4(c) is identical for both φ = π/6 and φ = π/2, as expected from symmetry. Finally, at φ = π/3 we do not observe any charge ordering, as can be seen from the structure factor S c (k) in Fig. 4(d). Instead, we clearly observe that for small values of |k x |. This is a key feature of a metallic state [29][30][31][32], as opposed to an insulating state which would be indicated by S c (k) ≈ αk 2 x . This agrees with the fact that the parent state at half-filling from ED in Fig. 1 is already metallic.

B. Effects of the flux φ on the hole-binding
The charge structure factor S c (k) in Fig. 4(c,d) illustrates a non-trivial effect of the flux φ on the charge degrees of freedom upon hole-doping the parent states at half-filling. To further investigate the effect of φ, we investigate the ground state electron density n R from DMRG for different values of φ in Fig. 5 for U/t = 8. Results are shown for hole-doping with n h = 4, 6, 8 holes and φ = 0, π/6, π/3. The density profile of φ = π/6 is identical to φ = π/2 due to symmetry. For φ = 0 in Fig. 5(a) we find that the holes remain separated from one another. Hence, we observe one hole per stripe as previously reported in Ref. [18]. In contrast, at φ = π/6 two holes bind together to form the charge modulation, so two holes per stripe are observed. Moreover, the holedensity at φ = π/3 does not show regular charge density wave modulations, consistent with the absence of a peak in the S c (k) in Fig. 4(d). Therefore, by changing the value of φ the system can be tuned from having repulsive, isolated holes to paired holes forming stripes to a more uniform charge density in a metallic state.
To quantify this effect of the flux φ on the charge degrees of freedom, we investigate the hole-binding energy x at small hole-doping for U/t = 8.0 in the case of (a) φ = 0, (b) φ = π/6, π/2, and (c) φ = π/3. The number of holes is denoted by n h . The φ = π/6 and φ = π/2 yield identical density profiles as guaranteed by the symmetry discussed in the main text. This density profile indicate that single holes are separated at φ = 0, while pairs of holes bind together to form a charge modulation at φ = π/6, π/2. At φ = π/3 no clear charge density wave patterns are formed.
E h−h and the electron-binding energy E e−e defined by, where E(m, n) denotes the ground state energy in the sector with m up-and n down-electrons. The dependence of these energies on the N s = 12 cluster from ED for various values of U/t is shown in Fig. 6. The holebinding energy E h−h is strongly positive at φ = 0 which implies that it is energetically more favorable to introduce two separate holes than a pair of holes, consistent with the isolated holes shown in Fig. 5(a). The value of E h−h decreases as a function of φ, eventually becoming negative and attaining a minimum between φ = π/6 and φ = π/3. Negative hole-binding energies indicate that binding of two holes is energetically prefereable to having two isolated holes. The small hole-binding energy at π/6 is, thus, consistent with having bound hole pairs as shown in Fig. 5(b). Due to symmetry, the electronbinding energies E e−e are identical to the hole binding energies E h−h up to a shift of φ = π/3.

C. Weak superconductivity
To determine whether the stripe state is accompanied by superconductivity, we investigate the pairing properties of the system. We focus on a particular set of parameters upon doping the 120 • magnetically ordered phase, U/t = 8, p = 1/18 for both φ = 0 and φ = π/6. Superconductivity is diagnosed by two means. First, we demonstrate off-diagonal (quasi) long range order in the pairing correlations. We consider the singlet pairing matrix ρ S (r i α|r j β), where α, β denote the direction of the nearest-neighbor lattice site on the triangular lattice and denotes the singlet-pairing operator. On elongated quasi one-dimensional geometries, long-range order even at T = 0 is excluded by the Mermin-Wagner theorem. Quasi long-range order, i.e. algebraic decaying correlation functions, is interpreted as an indication of true long-range order in the fully two dimensional system. We performed ground state simulations for bond dimensions D = 2000, . . . , 8000 and extrapolated towards infinite bond dimension. Results on the 72 × 3 cylinder are shown in Fig. 7. The extrapolated correlation funciton for φ = 0 is well-described by an exponential decay with a correlation length of ≈ 4.47, shown in Fig. 7(a). For φ = π/6 the correlation function is well described by either algebraic decay with an exponent ≈ 3.34 or an We extrapolate data from finite bond dimension D to infinite bond dimension by fitting a second order polynomial to the truncated weight ξ in DMRG, as shown in the insets. We fit both algebraic as well as exponential decay to the extrapolated correlations. The long-distance behavior of φ = 0 is well-described by an exponential decay with correlation length ≈ 4.47 while for φ = π/6 both an algebraic decay with exponent ≈ 3.34 and an exponential decay with correlation length ≈ 9.63 can be fitted.
exponential decay with correlation length ≈ 9.63. It is, thus, difficult to discern whether algebraic or exponential decay is realized based on the present data.
As compared to similar studies on different superconducting phases, this exponent is rather large (which means the superconductivity correlation is weak). Ref. [33], for example, reported an exponent η ≈ −0.96 for the superconducting state realized in the t-J model on square lattice on a width W = 4 cylinder. Hence, we interpret an exponent of η = −3.34 as a sign of weak superconductivity.
As a second means of diagnosing superconductivity, we investigate the spectrum of the non-local singlet density matrix, Upon Cooper pair condensation, we expect one or more eigenvalues to become dominant over the residual continuum of eigenvalues [34]. The spectra ofρ S (r i α|r j β) for U/t = 8, φ = π/6, and n h = 0, 6, 12, 18 on the 72 × 3 cylinder are shown in Fig. 8. For φ = 0 in Fig. 8(a) we observe simply continuum of eigenvalues and no dominant eigenvalues, indicating absence of superconductivity. Similarly, for φ = π/3 at half-filling n h = 0 shown in Fig. 8(b), we observe a continuum of eigenvalues. However, for n h = 6 we observe a small gap of three dominant eigenvalues to the continuum of residual eigenvalues. Similarly, we observe six dominant eigenvalues for n h = 12 and nine for n h = 18. Hence, the number of dominant eigenvalues equals the number of stripes in the system. While this phenomenon is analogous to the superconducting state in the two-dimensional t-J model on the square lattice, the gaps are smaller in magnitude. Whereas, for the robust superconductivity in the t-J model on the square lattice gaps of the order of ∆ ≈ 0.1 have been observed [34], here we report a gap of order ∆ ≈ 0.01. We again interpret this small gap as a sign of weak superconductivity. However, the spectrum in the hole-doped cases n h = 6, 12, 18 is clearly gapped in Fig. 8, as compared to the half-filled case or φ = 0 without superconductivity, which indicates a condensate of Cooper pairs forming on the stripes of the system.

IV. CONCLUSION AND DISCUSSION
The moiré Hubbard model is an effective description of the low-energy physics of twisted WSe 2 , which features a spin-dependent staggered flux through the plaquettes of a triangular lattice. This model exhibits a rich phenomenology whose central aspects we have now established by our combined ED and DMRG study. We have established an approximate phase diagram at half-filling, where at larger values of U/t two magnetic regimes with 120 • Néel order and one regime with xy-ferromagnetic order is realized. In comparison to the previous Hartree-Fock study [6], the critical values of U/t of the metalinsulator transition in shifted to larger values, in agreement with more elaborate studies on the pure triangular lattice Heisenberg model at φ = 0 [13,15]. At the particular values φ = π/6, π/2, 5π/6 the magnetic insulating phases extend up to U/t = 0 (within the numerical precision), which is also apparent from the single-particle gap in Fig. 2(c). Since at φ = 0 at intermediate values of 8 U/t 11 the model has been shown to feature a non-magnetic insulating (possibly spin liquid) state, a natural question is where else to expect putative spin liquid regimes. The spin gap in Fig. 2(d) is pronounced close to the metal-insulator transitions, which could be a first indication of a non-magnetic insulating state. However, further studies will be required to establish the phase diagram close to the metal-insulator transition at nonzero values of φ.
Doping the parent magnetic insulating states leads to the formation of intertwined spin and charge ordering such that the wavelength of the charge modulations is half the wavelength of the spin modulation at φ = π/6 and φ = π/2. Interestingly, a nonzero flux φ leads to the formation of hole-pairs which we relate to a strong dependence of the hole-binding energies on φ. While for φ = 0 individual holes are strongly repulsive, at intermediate and large values of U/t the hole-binding energy is found to be negative, leading to an attractive force between the holes. Analogously, we find that a finite value of φ can enhance the superconducting pair correlation and lead to a gap in the eigenvalues of the two-body density matrix. However, the pairing correlations at U/t = 8 and φ = π/6 are still weak and can be fitted by an algebraic decay with exponent ≈ 3.34 or an exponential decay with correlation length ≈ 9.63 in units of the lattice spacing, in either case too rapidly decaying to be consistent with a physical superconducting phase. Similarly, the gap in the spectrum of the two-body density matrix remains small. It will be interesting to determine how the superconductivity can be further enhanced, for example by further frustrating the magnetic order or tuning the value of the flux φ.
Pair-density wave (PDW) orders have been previously discussed in particular parameter regimes of Eq. (1). For φ = 0 in the case of the pure triangular lattice Hubbard model we observe an alternating sign of the pairing correlation ρ S (r, x) in Fig. 7(a), similar as has been observed in Ref. [18]. In our case, however, the absolute values |ρ S (r, x)| decay exponentially fast with a rather short correlation length. From this we cannot conclude that a PDW is realized at U/t = 8 and p = 1/18. At nonzero φ, PDW order has also been suggested for small values of U/t in DMRG [35] and renormalization group [36] studies. For the parameters studied in this manuscript, U/t = 8 and φ = π/6, π/2 we observe uniform pairing correlation inconsistent with a PDW state.
We note that, perhaps consistently with our finding of only weak superconductivity, no superconductivity has yet been observed in these materials, although superconductivity has been observed in the closely related twisted bilayer graphene materials [37]. This is in interesting counterpoint to the high -T c cuprate materials, a square lattice material family in which robust superconductivity is observed. For the cuprates the accepted theoretical model is the square lattice Hubbard model. Whether the ground state in certain parameter regimes is superconducting is subject of ongoing research [25] where an absence of superconductivity in the unfrustrated case has been noted [38]. However, in the closely related squarelattice t-J model robust superconductivity was recently established [33,34,39,40].
More generally, the great tunability of the moiré materials, in particular the ability to vary both carrier concentration and the hopping phase over wide ranges in situ, offers promise of a detailed comparison to theory. Displacement field-tunable metal insulator transitions with interesting precursor phenomena have been reported [1,4] along with indications of metallic magnetic phases [1]. The plethora of interesting phases found in our calculations encourage further experimental searches for the stripe and potentially superconducting phases predicted here. In favorable cases stripe phases may be observed via anisotropies in transport measurements although multidomain structures commonly occur and complicate the observations. The spatial modulation of the charge density occurring in a stripe may also be accessible to scanning capacitance probes. Deeper theoretical and experimental understanding of the spin liquid phase that may occur near U = 9t at φ = 0 and of the anomalous 'transition' phases separating the insulating magnet and non-magnetic metal phases are also important open questions.

ACKNOWLEDGMENTS
We would like to thank Edwin M. Stoudenmire and Matthew Fishman for insightful discussions. The DMRG calculations in this manuscript have been performed using the ITensor library [41]. The exact diagonaliza-  tion calculations have been performed using the Hydra library [42]. J.C., J.Z. and A.J.M acknowledge support from the NSF MRSEC program through the Center for Precision-Assembled Quantum Materials (PAQM) -DMR-2011738. The Flatiron Institute is a division of the Simons Foundation.