{'fraction_non_alphanumeric': 0.034327063959573, 'fraction_numerical': 0.01205035433424035, 'mean_word_length': 4.992474466941409, 'pattern_counts': {'":': 0, '<': 7, '<?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': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Renormalization group analysis of weakly interacting van der Waals Fermi system

Weak-coupling phenomena of the two-dimensional Hubbard model is gaining momentum as a new interesting research field due to its extraordinarily rich behavior as a function of the carrier density and model parameters. Salmhofer [{\it Commun. Math. Phys}. \textbf{194}, 249 (1998);{\it Phys. Rev. Lett}. {\bf 87}, 187004 (2001)] developed a new renormalization-group method for interacting Fermi systems and Metzner [{\it Phys. Rev. B} {\bf 61}, 7364 (2000);{\it Phys. Rev. Lett}. {\bf 85}, 5162 (2000)] implemented this renormalization group analysis of the two-dimensional Hubbard model. In this work, we demonstrate the spin-wave dependent susceptibility behavior of model graphene-phosphorene van der Waals heterostructure in the framework of renormalization group approach. We implement signlet vertex response function for the weakly interacting van der Waals Fermi system with nearest-neighbor hopping amplitudes. This analytical approach is further correlated with {\it ab initio} simulation results and extended for spin-wave dependent susceptibility behavior with possible experimental protocols. We present the resulting compressibility and phase diagram in the vicinity of half-filling, and also results for the density dependence of the critical energy scale.


I. INTRODUCTION
Hubbard model 1-7 becomes an effective tool to describe the many-body correlation effects in condensed matter physics [8][9][10][11] . It exhibits various different channels for the transition of normal metallic behavior by an independent tunability of the electronic band structure via a single tight-binding parameter and the half-filling nature on the hexagonal lattice bilayer system 8,10,12 . In this aspect, weak-coupling phenomena of the two-dimensional (2D) Hubbard model is gaining momentum as a new interesting research field because of its extraordinarily rich behavior as a function of the carrier density and model parameters 10 . Low-dimensional interacting lattice electron systems are generally featured with Fermi liquid instabilities driven by long-lived electron-electron or electron-hole excitations near the Fermi surface. These excitations lead to various types of long range correlations and destruction of Fermi liquid behavior.
However, for low dimensional Fermi systems, we recall the Renormalization Group (RG) is the most promising and best controlled analytical approach [13][14][15][16] . In 2D or bulk systems, RG methods are probably the least biased approach to study such systems at weak coupling and low energy scales in a transparent way [17][18][19] . This has led to the development of alternative flow schemes 20,21 . They provide the possibility to tackle the competing instabilities from different angles. Fortunately, it turns out that besides well-understood conceptual differences the flow of these schemes agrees quite favorably over a wider parameter range of models like the 2D Hubbard model. This connection of RG method and 2D Hubbard model encourages to further proceed on a firm approach to understand the instabilities in 2D system. In this regard, van der Waals (vdW) heterostructure systems 22 are coming under direct focus, where RG Theory and 2D Hubbard model can be implemented for practical realization to explore the concept of spin-wave assisted device application 23 .
To realize and analyze prototypic 2D model for the electronic degrees of freedom, Hubbard model including RG theory is promising in case of 2D planes 24,25 . In particular, 2D Hubbard model 26 is appropriate to explain interesting phenomena of weakly interacting systems which can be tuned by their carrier density. Besides, large scale numerical algorithms are worked out to solve such 2D interacting Hubbard model including related model parameters 27 . In such systems, weak interaction reveals antiferromagnetic ground state close to half-filling and symmetric s-wave phase away from halffilling region. In such cases, conventional perturbation theory fails at a certain region close to the half-filling region, where challenging and frequent infrared deviations raise due to weakly interacting Fermion system 17,25,28,29 and van Hove singularities 30 , non-smooth critical points in the density of states (DOS) of solids in their Brillouin zone, 31 including the 2D monolayer sheet system 32 . Moreover, appropriate model for 2D weakly interacting vdW Fermi systems with singularities within weakly interacting region is still lacking in the current analytical prospects. To overcome this fact, RG theory is presently the most efficient and promising analytical technique for such 2D weakly interacting Fermi systems 23,25,[33][34][35] . 2D Hubbard model is the simplest of a class of models describing electrons moving on a lattice having weakly-interacting electron-electron interactions. However, electron-electron interactions can break one or more symmetries of the Hamiltonian, and this leads to phases such as anti-ferromagnets, charge or density of spin-waves, etc. Besides, RG theory is presently the most promising and best controlled approach to low dimensional systems with competing singularities at weak coupling region. In this regard, direct device application through RG theory is unrealistic in pristine form, until and unless corroborated with vdWs heterolayer system using ab initio pseudopotential method. Thus, it is preferable to implement RG theory approach, unlike Dirac-like field theory to study the stability of the system supported by pseudopotential based density functional theory (DFT) simulation, unlike full potential basis of tight-binding method for practical realization supporting experimental scope.
In this work, we implement RG theory to the 2D Hubbard model of Graphene-Phosphorene (Grp-P) heterostructure considering nearest neighbor hopping amplitudes of dense electron distribution profile close to halffilling region. Here, we particularly calculate the trend of collective two particle couplings in a single loop stage, completely ignoring the inappropriate energy and linear momentum dependency, but considering vital tangential momentum dependency being a dynamic parameter. We further extend the RG theory to obtain the susceptibility values to analyze the physical response of uncertainties gestured by diverging couplings in weakly interacting vdW heterostructure. Charge-density screening, spindensity and singlet function based susceptibilities have been considered during analytical calculations at various pairing symmetries. Here, spin-density phase diagram has been presented close to half-filling region of this 2D Hubbard model. The obtained results show the critical energy scale behavior away from the half-filling region with perfect nesting and room temperature stability of the designed system for further functionality.

II. RENORMALIZATION GROUP EQUATIONS
In this section we revisit Salmhofer's renormalization group approach 13 for general weakly interacting Fermi systems with the explicit flow equations for effective two particle interactions on a one-loop level, and finally derive one-loop flow equations for several physical observables, like compressibility, spontaneous diamagnetic free susceptibility, etc., which will later be used for our stability analysis of the 2D Hubbard model.

A. Presentation of Functional Integral form
We have taken weakly interacting spin-half Fermions of single-particle basis state with crystal momentum (p c ), up or down spin-wave projection, σ [↑, ↓] and kinetic energy (E q ). Thus, the behavioral aspects of the system can be determined via following Eq. 1, (1) where, Q=(f 0 ,p,σ) presents a multi index parameter connecting the frequency (f 0 ) with each particle quantum numbers; Grassmann variables (χ Q and χ Q ) are allied to the formation and extinction operators, χ Q = Q -µ is the energy related to a particle depending on the chemical potential (µ) and ζ[χ, χ ] is a random many body coupling factor 36 .
In the calculation steps, all the connected and dependent transient flow equations have been acheived from the functional via implimenting Green's functions with the normalized Gaussian measure (given below in Eq. 2).
Green's function is an efficient and yet accurate tool to predict the variables in flow equations and maintain proper correlation among the RG theory equations. The infinite hierarchy in flow equations is essential to check at a single-loop level to predict dominant structural uncertainties in the vdW system within weakly interacting limit, neglecting other irrelevant components of the effective interactions (details shown in Supplimentary Material 37 ). and kF3 for the Grp-P heterostructure system with a comparison to the pristine graphene and phosphorene sheet. The model parameters are U=t and t=0, and the chemical potential |µ|=0.005.

B. Vertex Response Function
Renormalization Group method is presently the most promising and best controlled analytical approach to low dimensional Fermi systems with competing singularities at weak coupling region. Renormalization groups are commonly used as the existence of power-law singularities near critical points and half-filling regions. Here, the classic predictions are sometimes disturbed and logarithmic and exponential corrections are required in step-bystep process to address van der Waals low dimensional systems. To formulate the mechanism of Renormalization group theory is an interesting aspect to systematically explain the families of weakly-interacting Fermi systems united by common scaling variables. In the weak coupling case, the low-energy physics can be obtained via naive Hartree-Fock perturbation theory with weak coupling RG analysis.
Moreover, energy corrections of each iteration step are confined to one-particle vertex function (Γ) and used for further calculations. The spin-wave structure of Γ can be used for a spin rotation invariant system with response flow equation, which is given as follow in Eq. 3, where, K 1 and K 2 presents the multi indexed particle quantum numbers for both spin arrangements, singlet S and triplet T projection operators are used for two particle spin space, Λ is the energy scale varying at a functional trend for singlet and triplet operators.

C. Free susceptibility
We consider the spontaneous response to a field perpendicular to the plane. Here, the susceptibility follows like Eq. 4, Meanwhile, translation and spin rotation invariant systems are considered with conserved charge in the absence of an external field (F E ). As a consequence, in such normal and non-symmetrical broken phase, Φ(q) , the expectation value nullifies at F E (q) → 0.

D. Compressibility
In case of the density D = [n], the compressibility C = ∂[n]/∂µ D . The scaling relation then implies Eq. 5 where, x stands for a constant quantity, which leads to Eq. 6 with Eq. 7 This is used for the weakly interacting systems to evaluate the compressibility factor.

A. Linear Response Behaviour
We extend the RG equations to obtain linear response behaviour (i.e. spin-wave susceptibility value) of the system towards finite applied external electric field. Finite external field (F E ) leads to an additional contribution to the action of the flow equation, which is given as follow in Eq. 8 Thus, the Grassmann variables in bilinear form (Φ(q)), Hermitian conjugate (Φ (q)) and F E (q) is the complex conjugate of F E (q). One can calculate the response of the field interaction towards electronic charge density, ρ(q)= k,σ (ξ k−q,σ ξ k,σ ) and z-direction dependent spin density, s z (q)= k [ξ k−q,↑ ξ k,↑ -ξ k−q,↓ ξ k,↓ ] and response to coupling fields which interact with the change in singlet operator, where ∂(k) is an even parity function satisfying the orbital symmetry of the coupling operator (i.e. s wave). The linear response of the expectation value with the flow equation of the partition function.

B. Weakly interacting Fermi system
We study the effects of weak interaction in this 2D vdW system by RG theory analysis. It is obvious that over a wide range of energies the screened interaction is a relevant perturbation, and the system fails in the strong coupling regime displaying non-Fermi liquid behavior with power law energy dependence of the quasiparticle residue (discussed in the Supplimentary Material 37 in details). We describe a semimetal-semiconductor transition and associated quantum criticality in a 2D vdW system with particular emphasis on graphene monolayer sheet [38][39][40] . Moreover, the system enters into a weak coupling regime below a certain energy range. In this regime, only logarithmic renormalizations are possible 25,41-45 . As a consequence, all the possible calculations have been performed for physical observables, like compressibility, spontaneous free susceptibility, etc. 42,46,47 . Other physical observables, like free susceptibility, phase transition, etc. are discussed in the Supplimentary Material 37 in details.

IV. VERTEX FUNCTION FINDINGS
A. Vertex response function corrections One can verify that higher-order corrections contain higher powers of the upper limit 48 . The presence of anomalous dimensions in the RG theory proves that physical observables consdiered here int he work, such as spontaneous susceptibility, compressibility, etc., show indirect dependencies on response function. This strong coupling results differ quantitatively, but not qualitatively from the case of pristine graphene.
The flow of the vertex function (Γ) (Eq. 3) and the free susceptibility (χ F ) (Eq. 4) values have been calculated analytically using RG theory for several choices of the bare interaction, termed as field potential (U>0), the next nearest neighbor hopping amplitude (t<0) and the chemical potential (µ), where t and µ have been fixed so that the Fermi surface and particle density are approaching to k and half-filling (i.e. n≈1, |µ| →0), respectively. Momentum dependent flow of the vertex functions are obtained for the three considered systems at low dimensional scale. It is observed that the flow of the vertex function achieves a density value (n≈0.992), i.e. close to half-filling region at t=0 and |µ|=0.005 (plotted in Fig. 1). Here, the singlet part of the vertex function is taken at specified values of momenta on the Fermi surface including the values of momenta, where strong renormalization of Γ occurs. Most probably, the singlet vertex function possesses its leading values supporting umklapp scattering along the diagonal direction of the kspace Brillouin zone. It is observed from the RG theory calculation that the density decreases away from halffilling (i.e. n≈1, |µ| →0) region in a step-like manner and the calculations fail into a region, where s-wave symmetry dominates with pairing correlations at necessarily low scale energy levels indicating quantum confinement (i.e. nanostructuring) in the Grp-P heterostructure system. As a result, differential step like behaviour is prominent from the plots in Fig. 1.

B. Spin-wave dependent susceptibility behaviour
Two types of antiferromagnetic spin-wave (i.e. incommensurate and commensurate) susceptibilities have been obtained from the numerical calculations to validate the physical instability associated with the diverging vertex function as a consequence of response functions 49 . It is observed that the incommensurate spin susceptibilities can be neglected compared to their respective commensurate counterparts (shown in Fig. 2) due to the incommensurability parameter (δ →0) nearing to half-filling (n≈1, |µ| →0) as a function of energy scale. Moreover, the ratios between the two susceptibility values of monolayer system in case of weakly interacting s-wave region are in the same scale, whereas the ratios of free susceptibility are much reduced indicating the validity of response functions for 2D weakly coupled Fermi systems. In addition, the ratio in case of hetero bilayer is much greater compared to its monolayer system. Thus, antiferromagnetic wave vector based spin-wave susceptibility dominates over coupled susceptibilities at low scale energy level (i.e. ground state), confirming the existence of antiferromagnetic ground state with rotational spin invariance.

C. Spin-wave dependent compressibility behaviour
The values of compressibility (Eqn. 5) and homogeneous spin susceptibility have been calculated depending on the energy scale (shown in Fig. 3(a) and Fig. 3(b)), respectively for the three systems. These two parameters can be realized directly from the forward scattering of vertex function. The free compressibility, κ 0 , and spin-wave susceptibility, χ(s,0), are defined ignoring the probable infrared cut-off energy values. As a result, the flow equations of the vertex response function is entirely checked by the flow of the Landau functions at the initial stage of the Hubbard model at critical energy scale value, Λ=Λ 0 .
The compressibility is inhibited at low scale energy levels close to half-filling region, which is expected for systems having a finite bandgap near to µ, with dominant uncertainty in spin-wave density. Besides, the values of the compressibility diverge with respect to the reduc- tion in homogeneous spin-wave susceptibility values away from half-filling region, where s-wave uncertainties play major role indicating a finite spin-wave bandgap opening trend in any singlet or triplet spin-wave 2D Fermi systems. It is found that the spin-wave susceptibility flows from zero to negative values which implies that our single loop calculation fails, when it reaches strongly interacting regime. This trend is validated only in weakly interacting region applied for 2D Fermi system. Thus, one diverging compressibility is found mentioning a strong tendency toward phase transition and separation. As a result, the increase in κ is very close to the uncertainty point, where the renormalized interactions achieve large value that the single loop findings are not further validated. The stability is prominent at higher energy scales (i.e.

0.04≤
Λ t ≥0.03) from the plot in Fig. 3(a) and Fig. 3(b). Moreover, we have shown the trend of free susceptibilities (χ 0 ) at specific choice of t and µ values (shown in Supplimentary Fig. S1 37 ), strongly supporting the trend of homogeneous spin-wave susceptibility.

D. Spin-wave density dependent phase diagram
We show the phase diagram of µ and U at half-filling measured by the prevailing uncertainty from the flow equation (shown in Fig. 4). The leading uncertainty region in case of commensurate spin-wave density is detached from the s-wave coupling region via a narrow strip where incommensurate density fluctuations show dominant effect with q=(π, π-δ). For small value of U, the region surrounding half-filling is exponentially small, where uncertainties in spin-wave density play a dominant effect. The value of (t, U), in-plane phase diagram is shown in Fig. 4 with |µ|=4t and µ<0 (below half-filling). It is observed that the chemical potential is fixed at the van Hove singularity 50 irrespective of the value of spin-wave density, which is decreasing away from half-filling region with growing value of |µ| at 0.048, one critical point in the flow equation shown in FIG. 4. Spin-wave density dependent phase diagram for (µ, U) phase diagram for t=0 close to half-filling (n≈1, |µ| →0) region; the solid line represents the spin-density-wave regime of the Grp-P heterostructure system with a comparison to the pristine graphene and phosphorene sheet. Inset is showing the heterostructure system. Supplimentary Fig. 2 37 ). We observe from the phase diagram of µ and U at half-filling measured by the prevailing uncertainty from the flow equation (shown in Fig. 4). The leading uncertainty region in case of commensurate spin-wave density is detached from the s-wave coupling region via a narrow strip where incommensurate spin-wave density fluctuations show dominant effect with q=(π, π-δ). For small value of U, the region surrounding half-filling is exponentially small, where uncertainties in spin-wave density play a dominant effect. The value of (t, U), in-plane phase diagram is shown in Fig. 4 with |µ|=4t and µ<0 (below half-filling).

E. Behavioral trend of critical energy scale
The behavioral trend of critical energy scale, Ξ c , depending on the function at t<0 where, |µ| →0 and U=t is plotted in Fig. 5(a). The shrinkage in the value of Ξ c with growing trend of |µ| confirms that the E F residues on the van Hove singularity within the weakly interacting regime having only Fermi surface nesting dominance. Momentum variables of Γ have been projected upon the Fermi surfaces at specific points, where the strong renormalization of the vertex function occurs. Here, total 8 fixed points considered for the normalization process. Besides, additional 8 points have been considered in a more refined projection scheme upon the van Hove surface to check the collective effect.
As a result, sum total 16 points (8 for Fermi surface and 8 for van Hove surface) have been taken during the flow equation calculations. These refinements, in the response function of the spin-wave density scheme, enhance the analytical calculation efficiency significantly with a reasonable reduction in critical energy scale, without altering Γ and χ values, qualitatively. The variation trend of Λ c as a function of the inverse number of discretization points, N 0 , upon the Fermi surface at fixed choice of model parameters is shown in Fig. 5(b). It is noticed that the critical energy scale is varying in a discrete manner with 8 discretization points, and resumes stability for next 8 points with right order of magnitude indicating the controlled trend of flow equation for weakly interacting systems. This last result shows that the weak coupling behavior of Grp-P system is qualitatively different from that in pristine graphene. In graphene, the quasiparticle Y factor tends to a finite value at zero energy, i.e. the system retains Fermi-liquid behavior with well defined quasiparticles, on the contrary, do not become sharp quasiparticles, even at the lowest energies.

V. EXPERIMENTAL SCOPE OF THIS STUDY
In this presnt work, we call the Renormalization Group (RG) analysis as the most promising and best controlled analytical approach for low dimensional Fermi systems. In 2D or bulk systems, RG methods are probably the least biased approach to study such systems at weak coupling and low energy scales in a transparent way. This has led to the development of alternative flow scheme. It provides the possibility to tackle the competing instabilities from different angle. It agrees quite favorably over a wider parameter range of models like the 2D Hubbard model. This connection of RG method and 2D Hubbard model encourages to further proceed on a firm approach to understand the instabilities in 2D system. In this present study, we have considered van der Waals heterostructure system, where RG Theory and 2D Hubbard model are implemented for practical realization to explore the concept of spin-wave assisted device application supporting strongly towards the experimental realization.
To realize the practical application of this system, we demonstrate the hetero bilayer dual gate field effect transistors (DG-FETs) with high on-off ratio considering Grp-P bilayer as a model system (shown in Fig. 6(a)). The mechanism of Giant Stark Effect is implemented here to design the vdW heterolayer-based dual-gate FET, where unsaturated Grp-P bilayers are acting as electrodes 51,52 at room temperature. In this aspect, direct device application through RG theory is unrealistic in pristine monolayer form, until and unless corroborated with heterolayer system using ab initio pseudopotential method (refer to 37 ). It is not a secret that there is no experimental evidence in condensed matter physics for the realization of the model solenoid that could produce the field in monolayer system. Besides, the ability to realize such specified device can be proceeded via perturbed system using a pseudopotential method (i.e. DFT simulation), which is similar to the vector field potential 53 . Here, saturated Grp-P bilayers serve as electrode-cum-tunneling barriers in the crucial part of the model device having dual (i.e. top and bottom) gates to generate a transverse electric field. Besides, several devices utilize an orbital field to achieve spin splitting 54,55 via the interference of the spinor wave functions in the 2D systems. Moreover, few more reports are there showing implications of this vector field corroborated with ab initio DFT simulations both experimentally and theoretically for device applications 56,57 . It is evident that all phosphorus based FETs can neglect the metal to semiconductor contact at the interface for this transport behaviour. The transmission band spectrum of the model dual-gate FET at 0 V/nm and 5 V/nm, electric fields free from source to drain voltage is shown in Fig. 6(b) at 300 K. The semiconducting characteristic of the saturated Grp-P heterolayer blocks the transmission process near the Fermi level at 0 V/nm electric field, originating a finite transmission gap of 1.12 eV, as a result, there is no transmission states near Fermi level. This absence of states indicates the off state of the device. Similarly, a transmission peak emerges near the Fermi level with a wide range of varying dispersion (i.e. -0.07 to +0.08 eV), at 5 V/nm electric field value, which is an unusual Van Hove-like singularity in 2D materials. This supports the on state of the designed dual-gate FET, which is stable at 300 K. The transmission eigenchannels at E and at the Γ (0, 0) point of the reciprocal k-space, presented in Fig. 6(c) & (d), intensely illustrate the off and the on state of the designed FET controlled by a dual-gate induced electric field, respectively.
Similarly, an external transverse electric field of 0 V/nm, the transmission eigenvalue achieves the value . The ratio of the on-off state is found to be 10 2 which indicates the persistent of quantum confinement effect due to nanostructuring in the vdW heterostructure. The value of this model hetero bilayer based dual-gate FET is 3 order higher than pristine phosphorene FETs 58 and 1 order higher compared to graphene-MoS 2 based FETs at nanoscale 58,59 . The calculated density of states of the scattering region at 5 V/nm as an external electric field is inset in Fig. 6(b). The calculated DOS pattern indicates a distinct finite peak near the Fermi level, implying a strong correla-tion between the transmission channels. The physics behind such strong correlation shows that the transport phenomena is mainly dominated by resonant tunneling through interface states [shown in Fig. 6(c) and (d)], not barrier tunneling at the Fermi level.

VI. SUMMARY AND CONCLUSIONS
A study on the spin-wave response function behavior of vdW Grp-P heterostructure was presented via fieldtheory perspective implementing analytical RG theory. This formulation can be used systematically to identify coupling uncertainties in 2D weakly interacting Fermi systems. In this system, critical energy can be scaled at specific singular points where vertex functions and susceptibilities diverge and become finite. Moreover, the instabilities are modulated via short range antiferromagnetic correlations present in the system, which is noticed from the flow of vertex functions and spontaneous susceptibility. The analytical results show an anomalous singularity of van Hove kind in 2D heterostructure. This finding supports the qualitative understanding of vdW heterostructure systems from the anomalous critical points of the response function. This will also facilitate quantitative analysis of room temperature stability of such systems and on-state electronic device modeling in future with a step towards experimental realization.

SUPPLEMENTARY MATERIAL
See the Supplimentary Material for the physical variables and related equations relating with the phase diagram of the graphene-phosphorene heterostructure system.

SUPPLEMENTARY MATERIAL
Rational control of individual molecular spins in nanoelectronics devices is a pivotal prerequisite to fulfill its potential promised by molecular spintronics 41,60 . There is a rapidly growing interest in the transmission of information encoded by electron and atomic spins which, in principle, can be realized by the propagation of spin waves or magnetization waves through low-dimensional sheets. In this aspect, weak-coupling phenomena of the two-dimensional Hubbard model is gaining momentum as a new interesting research field recently because of its extraordinarily rich behavior as a function of the carrier density and model parameters 42,61 .
In this regard, establishing an effective ab initio Schrodinger equation for spin waves incorporating analytic response functions with suitable boundary conditions for spin wave controlled transmission is absolute requirement to validate the theoretical understanding about antiferromagnetic spin-wave field-effect transistor [43][44][45] . In this aspect, this report is the first of its kind where we propose that the spin wave assisted field-effect transmission may be selectively used for controlling individual molecules 46 embedded in the dual-gate field-effect transistor architecture.

Free susceptibility
Moreover, we have shown the trend of free susceptibilities (χ 0 ) at specific choice of t and µ values (shown in Fig.  7), strongly supporting the trend of homogeneous susceptibility.
FIG. 7. Free (non-interacting) susceptibilities of Grp-P heterostructure system with a comparison to the pristine graphene and phosphorene sheet for t=0 and |µ|=0.005.

Phase Transition
Let us consider the phase transition in the weakly coupled heterostructure due to the instabilities prevailing on the surface.
FIG. 8. Spin-wave density dependent phase diagram for (a) The (µ, U) phase diagram for t=0 close to half-filling (n≈1, |µ| →0) region; the solid line represents the spin-density-wave regime of the Grp-P heterostructure system with a comparison to the pristine graphene and phosphorene sheet. Similarly, (b) The (t, U) phase diagram for 4t<0 value below the half-filling (|µ| <0) region; the dotted lines represent the spin-density-wave regime of the same three systems.
We show the phase diagram of µ and U at half-filling measured by the prevailing uncertainty from the flow equation (shown in Fig. 8(a)). The leading uncertainty region in case of commensurate spin-wave density is detached from the s-wave coupling region via a narrow strip where incommensurate spin-wave density fluctuations show dominant effect. For small value of U, the region surrounding half-filling is exponentially small, where uncertainties in spin-wave density play a dominant effect. The value of (t, U), in-plane phase diagram is shown in Fig. 8(b) with |µ|=4t and µ<0 (below half-filling).
In this case, all the steps are required for a clear udnerstanding of the facts. As a result, rational control of individual molecular spins in nanoelectronics devices is a pivotal prerequisite to fulfill its potential applicability. Meanwhile, there is a rapidly growing interest in the transmission of information encoded by electron and atomic spins which, in principle, can be realized by the propagation of spin waves or magnetization waves through low-dimensional systems. In this aspect, weak-coupling phenomena of the two-dimensional Hubbard model is gaining momentum as a new interesting research field recently because of its extraordinarily rich behavior as a function of the carrier density and model parameters. We show that this finding supports the quantitative realization of experimental and theoretical systems farther from their critical points before the underlying field theory is well understood for potential application.

First Principle based DFT Calculations
QUANTUM Espresso codes 62 based on Perdew-Burke-Ernzerhof (PBE) with generalized gradient approximation (GGA) is used in our DFT calculations 63 in the framework of local spin density approximation (LSDA) i.e. all calculations are carried out taking into account the presence of electron spins with time-reversal (TR) symmetry breaking along z-axis. Supercells having larger than 20Å lattices are considered for first principle calculations to avoid periodic images. Brillouin zone integration has been performed with a 3×3×1 k points grid sampling at an energy cutoff of 60 Ry based on the Monkhorst and Pack scheme 64 . Densities of states have been obtained on a k-point sampling of 9×9×1 with smearing effect at 540 Ry as charge density cut-off energy. The supercell structure has been optimized till the force and energy per each atom attains the value within 0.01 eV/Å and 2.0×10 −7 eV, respectively. Three types of systems, such as graphene monolayer, phosphorene monolayer and graphene/phosphorene bilayer have been considered in our calculations. The bilayer heterostructure has been made of the host (i.e. phosphorene) on graphene substrate (i.e. bilayer of graphene/phosphorene). In this case, the supercell, 4 √ 3×4 √ 3 graphene on √ 27 × √ 27 BP with a rotation angle of 12 0 was chosen with 48 Carbon atoms in graphene layer. The lattice mismatch is less than 0.24 %. In case of monolayer systems, same number of atoms are taken in the simulation cell.