Flavien Thierry

Giuseppe Rastiello
email: giuseppe.rastiello@cea.fr

Cédric Giry

Fabrice Gatuingt

One-dimensional Eikonal Non-Local (ENL) damage models: Influence of the integration rule for computing interaction distances and indirect loading control on damage localization

Keywords: Eikonal Non-Local (ENL) damage model, Damage dependent interactions, Integration formulas, Path-following constraint equations, Regularization

 assumed that the damage field is linear between nearby integration points. Mainly motivated by improved integration accuracy, this method leads to a more gradual damaging process and makes it easier to follow the structural response when damage localizes. In this contribution, we show that despite these advantages, such a choice induces parasite damage diffusion issues, making the use of a standard trapezoidal rule more relevant. Moreover, once this choice is made, the structural response during damage localization can be conveniently described using a path-following method based on controlling the non-local strain variation of the localizing element instead of the local strain as in the works cited above.

Introduction

Continuum damage models are often used to represent failure processes in different kinds of materials. From a mathematical viewpoint, it is well known that material softening leads to the ill-posedness of the equilibrium problem. In a numerical framework, this translates into spurious energy dissipation and non-objective results with respect to the spatial discretization of the considered problem.

Several methods were proposed in the literature as localization limiters. One can cite, among them, Integral Non-Local (INL) formulations on the internal variables of the constitutive model [START_REF] Pijaudier-Cabot | Nonlocal Damage Theory[END_REF][START_REF] Ganghoffer | A reappraisal of nonlocal mechanics[END_REF][START_REF] Borino | A symmetric nonlocal damage theory[END_REF][START_REF] Grassl | Damage-plastic model for concrete failure[END_REF], gradient enhanced models [START_REF] Frémond | Damage, gradient of damage and principle of virtual power[END_REF][START_REF] Peerlings | A critical comparison of nonlocal and gradient-enhanced softening continua[END_REF], phase-field formulations [START_REF] Francfort | Revisiting brittle fracture as an energy minimization problem[END_REF][START_REF] Bourdin | The variational approach to fracture[END_REF][START_REF] Miehe | Thermodynamically consistent phase-field models of fracture: Variational principles and multi-field FE implementations[END_REF][START_REF] Verhoosel | A phase-field model for cohesive fracture[END_REF], or Thick Level Set models [START_REF] Moës | A level set based model for damage growth: The thick level set approach[END_REF][START_REF] Bernard | Damage growth modeling using the thick level set (tls) approach: Efficient discretization for quasistatic loadings[END_REF][START_REF] Moreau | Concurrent development of local and non-local damage with the thick level set approach: Implementation aspects and application to quasi-brittle failure[END_REF][START_REF] Lé | Coupling damage and cohesive zone models with the thick level set approach to fracture[END_REF]. Despite their differences, the common feature of these theories is that a characteristic length is introduced into the formulation allowing to diffuse dissipation over a zone of finite size. From a numerical viewpoint, this enables recovering mesh insensitivity.

In particular, INL models [START_REF] Pijaudier-Cabot | Nonlocal Damage Theory[END_REF][START_REF] Bažant | Nonlocal Integral Formulations of Plasticity and Damage: Survey of Progress[END_REF] are widely used due to their simplicity of implementation and numerical robustness. According to these formulations, damage evolution is driven by a non-local field (e.g., a non-local equivalent strain field), obtained by weighted averaging of the corresponding local field over the domain. In standard theories, averaging is performed such that non-local interactions between material points are controlled by the Euclidean distance between them (i.e., the higher the distance, the lower the interaction).

Limitations and deficiencies of INL formulations have been pointed out by many authors (e.g., [START_REF] Geers | Strain-based transient-gradient damage model for failure analyses[END_REF][START_REF] Jirásek | Size effect on fracture energy induced by non-locality[END_REF][START_REF] Simone | Incorrect initiation and propagation of failure in non-local and gradient-enhanced media[END_REF][START_REF] Bažant | Nonlocal boundary layer (nbl) model: overcoming boundary condition problems in strength statistics and fracture analysis of quasibrittle materials[END_REF][START_REF] Giry | Stress-based nonlocal damage model[END_REF][START_REF] Grégoire | Failure and size effect for notched and unnotched concrete beams[END_REF]). In particular, considering fixed non-local interaction distances leads to non-physical interactions through holes, cracks, and developing damaged bands. This results in parasite damage spreading across damaged bands, near notches, and free edges in simulations.

Enhancements of the original methods were proposed to face these limitations (e.g., [START_REF] Nguyen | A damage model with evolving nonlocal interactions[END_REF][START_REF] Giry | Stress-based nonlocal damage model[END_REF][START_REF] Krayani | Boundary effect on weight function in nonlocal damage model[END_REF][START_REF] Pijaudier-Cabot | Non local damage model[END_REF]) and to describe more precisely strain localization processes in softening media through considering evolving non-local interactions.

In the Eikonal Non-Local (ENL) formulation [START_REF] Desmorat | Nonlocal models with damagedependent interactions motivated by internal time[END_REF], interaction distances between material points are computed by solving an isotropic time-independent Eikonal equation with a damage dependent metric function. In two-and three-dimensional (2D and 3D) settings, the ENL formulation considers that damage induces a curvature of the Riemannian space in which interaction distances are computed. This mathematical/physical framework allows for directly modeling evolving interactions during the localization process, thus naturally leading to progressive damage -to -fracture transition. In [START_REF] Desmorat | Nonlocal models with damagedependent interactions motivated by internal time[END_REF][START_REF] Jirásek | Localization analysis of nonlocal models with damage-dependent nonlocal interaction[END_REF] (uni-dimensional case, 1D) and [START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF] (two-dimensional case, 2D) it was shown that an isotropic ENL damage formulation allows for representing the main phases of strain localization process in quasibrittle materials, including the progressive transition from diffuse to localized damage. Furthermore, the regularization properties of standard INL formulations are preserved until high damage levels.

In the 1D setting that is considered in this work, finding the interaction distance ˜ xx = ˜ x (x ) between two material points x and x belonging to a domain Ω consists in solving the eikonal equation:

1 -D(η) d dη ˜ x (η) = 1 (1)
˜ x (η = x) = 0 (2)
with D(η) ∈ [0, 1) the isotropic scalar-valued damage field. After integrating between x and x one has:

˜ xx = max (x,x ) min (x,x ) dη 1 -D(η) ≥ xx (3)
where xx = x (x ) = |x -x | ≥ 0 is the Euclidean distance between x and x (i.e., the interaction distance that is considered in INL models) and min(•)/ max(•) are the standard minimum/maximum operators.

In the framework of the Finite Element Method (FEM), one, therefore, has to compute interaction distances between integration points where damage is evaluated. A simple trapezoidal rule was used by Desmorat et al. [START_REF] Desmorat | Nonlocal models with damagedependent interactions motivated by internal time[END_REF] to calculate integral (3), under the assumption that the damage field is constant between adjacent integration points. More recently, Jirásek and Desmorat [START_REF] Jirásek | Localization analysis of nonlocal models with damage-dependent nonlocal interaction[END_REF] proposed to use a slightly modified integration rule based on assuming that the damage field is linear between adjacent integration points. According to these authors, such a choice leads to improved accuracy in integral computation and makes the damage process more gradual, making it easier to follow the structural response when damage localization occurs.

Based on a localization problem into a bar submitted to tensile loading, in this short communication, we show that such a choice induces parasite damage diffusion, leading to loose one of the main advantages of ENL formulations. This makes the use of a standard trapezoidal rule more effective. Moreover, the structural response in the final phases of the test can be described by using a dissipative path-following method based on controlling the non-local strain variation of the localizing element instead of the local strain (as in the articles cited above).

Eikonal non-local damage model

Let us consider a 1D body occupying all points x ∈ Ω = ]0, L[. The mechanical behavior of its constituting material is described through an isotropic continuum damage model [START_REF] Lemaitre | Engineering Damage Mechanics: Ductile, Creep, Fatigue and Brittle Failures[END_REF] with a single scalar variable D. Accordingly, the Cauchy's stress reads σ = σ(x, t) = (1 -D)Eε, where E = E(x) is the Young's modulus, ε = ε(x, t) = du/dx is the infinitesimal strain and u = u(x, t) is the displacement field.

Non-local damage evolution is supposed to occur according to the law, D = g(κ), where κ = κ(x, t) is an internal variable. Its evolution has to fulfill the following non-local Kuhn-Tucker and consistency conditions: f ≤ 0, κ ≥ 0, f κ = 0 and ḟ κ = 0.

In previous definitions, f = f (σ, ε, D) is a damage activation function, and symbol • denotes the rate of variation of the quantity •. We chose f = ε -κ, where ε = ε(x, t) denotes the non-local strain. Alternatively, instead of ε, one can define f in terms of the non-local counterparts of the thermodynamic force associated with damage Y = Y (x, t) = Eε 2 /2 [START_REF] Marigo | Formulation d'une loi d'endommagement d'un materiau élastique[END_REF] or other quantities [START_REF] Jirásek | Nonlocal models for damage and fracture: Comparison of approaches[END_REF].

Integral nonlocal (INL) formulation

In standard INL-type formulations [START_REF] Pijaudier-Cabot | Nonlocal Damage Theory[END_REF], field ε is computed by weighted averaging of field ε over the entire domain:

ε = ε(x, t) = Ω φ(ξ xx )ε(x , t) dx Ω φ(ξ) dx (4)
where ξ xx = ξ(x, x ) denotes the ratio of the interaction distance xx to the characteristic length ( c > 0), and φ(ξ xx ) is a positive monotonic decreasing function of ξ xx . The latter function is chosen such that φ(ξ xx ) = 1 for ξ xx = 0, and φ(ξ xx ) → 0 when ξ xx increases. A common choice to express φ(ξ xx ) consists in using the Gaussian distribution function, i.e., φ(ξ xx ) = exp(-4ξ 2 xx ).

Eikonal non-local (ENL) formulation

According to the ENL approach, the non-local averaging process is the same as in the standard INL formulations. However, non-local interactions change according to the damage field evolution through effective distances. In a 1D setting, the length ratio ξ xx appearing in eq. ( 4) thus reads:

ξ xx = ˜ xx c ≥ |x -x | c (5)
According to eqs. ( 3) and [START_REF] Bourdin | The variational approach to fracture[END_REF],

˜ xx = |x -x | when D(η) = 0 ∀ η ∈ [x, x ] (i.e.,

Numerical formulation in 1D

Consider a domain Ω that is discretized through a FE mesh (Ω h ) comprising an odd number (n el ) of FEs of equal size (h = L/n el ).

A continuous linear interpolation is chosen for the unknown displacement field. Integration points (one per element) are identified according to the ordered sequence

I = {x i } n el i=1 . The discretized damage field {D i } n el i=1 = {D(x i )} n el i=1
, is evaluated at integration points according to an explicit procedure.

Interaction distances ( ˜ i,j ) between points x i and x j ∈ I are computed at the beginning of each time step based on the last converged damage field (at time t n ), {D n i } n el i=1 . They are then used to update the non-local strain field (4) driving damage evolution and kept constant until convergence of the global equilibrium problem (at time t n+1 ).

Interaction distances computation

As illustrated by Jirásek and Desmorat [START_REF] Jirásek | Localization analysis of nonlocal models with damage-dependent nonlocal interaction[END_REF], different integration techniques can be used for computing eq. ( 3), and might influence the obtained results:

(i) the first and natural choice consists in assuming that function √ 1 -D is element-wise constant, and performing integration according to a trapezoidal rule. The effective distance ˜ i,j between points x i and x j > x i thus reads:

˜ i,j = ˜ i,j-1 + h 2 1 1 -D j-1 + 1 1 -D j (6)
where ˜ i,j-1 is the effective distance between x i and x j-1 ;

(ii) an alternative way for computing interaction distance was recently proposed in [START_REF] Jirásek | Localization analysis of nonlocal models with damage-dependent nonlocal interaction[END_REF]. Assuming that the damage field is linear between adjacent integration points, distance ˜ i,j reads:

˜ i,j = ˜ i,j-1 + 2h 1 -D j-1 + 1 -D j (7)
According to the cited work, this formulation shows improved performance from the integration error viewpoint, in particular when damage localization occurs. The influence of such a choice on the damage localization process has not been studied yet, however.

In the following, the ENL formulation with trapezoidal rule (eq. ( 6)) is denoted ENL1, whereas the other one (eq. ( 7)) is named ENL2. For considerations holding for both methods, we will use the generic term ENL.

Considered problem

We consider a bar of length L and constant transverse crosssection S submitted to tensile loading. Displacement is blocked on the left side of the bar (x = 0), whereas an indirectly controlled force is imposed at the right extremity (x = L) to ensure describing structural responses characterized by snap-back phases.

Different FE meshes are considered (n el ranging from 41 to 2001 linear Bar elements). Damage initiation is forced to occur on the central/weak element (the corresponding integration point denoted x w , with w = (n el + 1)/2). For this purpose, the Young's modulus is spatially distributed as

E i = E(x i ) = E 0 (1 -α i /10)
, where E 0 is a reference Young's modulus, and α i = α(x i ) is chosen as a Gaussian function. centered on x w and with length parameter equal to L/10, i.e., α i = exp -400(x i -x w ) 2 /L 2 . This allows for removing the small differences in the structural response (stiffness and load peak) that are observed between different meshes when localization is induced by slightly reducing E(x w ) only.

Damage evolves according to the exponential law

D i = g(κ i ) = 1 -(κ 0 /κ i ) exp (-κ i -κ 0 /(κ c -κ 0 )), where κ i = κ(x i ),
κ 0 is the non-local strain level corresponding to damage activation, κ c is the non-local strain level controlling the shape of the damage evolution function and • denotes the Macaulay operator (given a real number

•, • = • if • ≥ 0 and • = 0 oth- erwise)
. Different evolution models (linear, bi-linear, exponential. . . ) could have alternatively be used without substantially changing the theoretical formulation and its numerical implementation.

A quasi-newton method is used to solve the global equilibrium problem. To follow the structural response in the snap-back phase, the external loading is controlled indirectly based on a simple path-following technique [START_REF] De Borst | Nonlinear Finite Element Analysis of Solids and Structures[END_REF][START_REF] Rastiello | Discontinuity-scale pathfollowing methods for the embedded discontinuity finite element modeling of failure in solids[END_REF].

Two path-following constraint equations are considered:

(i) the first one imposes the local strain variation of the weak element (∆ε w -∆τ = 0, with ∆τ > 0 a user-defined path-step length) between two successive loading steps, as in [START_REF] Desmorat | Nonlocal models with damagedependent interactions motivated by internal time[END_REF][START_REF] Jirásek | Localization analysis of nonlocal models with damage-dependent nonlocal interaction[END_REF];

(ii) the second one imposes the non-local strain variation of the weak element (∆ε w -∆τ = 0, ∆τ > 0). Given the simple considered geometrical and mechanical models, the chosen constraint ensures achieving dissipative responses only since the maximum non-local strain variation is always attained at the weak element.

More details concerning the numerical solution procedure are given in Appendix A.

Results and discussion

Global responses

Let us first consider a finite element mesh with 81 Bar FEs and compare force-displacement (F, u L ) responses obtained by considering the local, INL and ENL damage models (fig. 1), for both integration formulas introduced in section 2.4.

As already observed by Rastiello et al. [START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF], ENL and INL models provide very similar responses for moderate damage levels (non-local interactions are mainly controlled by Euclidean distances in this phase). Then, ENL responses diverge from the INL one and, as expected, tend to the local one when damage grows up to unity on the weak element (D w → 1). The ENL1 and ENL2 formulations give almost identical results until very high damage levels (approximately D w ≈ 0.9), then they exhibit some differences. In particular, the ENL2 formulation provides a smoother failure of mechanical strength than ENL1 and a higher energy dissipation. In both cases, however:

(i) controlling the simulation through imposing ∆ε w (solid lines in fig. 1) makes it possible to properly follow the whole damage localization process;

(ii) conversely, using a constraint equation based on ∆ε w (dashed lines in fig. 1) -or the relative displacement of the nodes of the weak element as in [START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF][START_REF] Jirásek | Localization analysis of nonlocal models with damage-dependent nonlocal interaction[END_REF] -makes it very difficult to control the simulation when D w → 1, in particular for the ENL1 formulation. And this is independent of the chosen ∆τ . By the way, this problem partly motivated Jirásek and Desmorat [START_REF] Jirásek | Localization analysis of nonlocal models with damage-dependent nonlocal interaction[END_REF] to replace the trapezoidal integration rule by eq. ( 7). Notice that for the considered mesh no convergence issues were experienced; however, this can append for finer meshes (e.g., for n el = 161 or higher). See Appendix A for more details.

From here-forth all results are obtained using the constraint equation based on the non-local strain variation. Some additional explanations concerning this choice will be given in the following.

Damage localization

As illustrated by the damage distributions along the bar (fig. 2), both ENL implementations model a progressive damage localization in the central region of the bar.

Based on these results, at any pseudo-time, one can identify three different zones in the bar: (i) an undamaged region (I elastic ), where the material response stays elastic and no damage occurs; (ii) an inactive damage region (I inactive ), where material points are damaged but no more experience damage evolution;

(iii) an active damage region (I active ), where damaging is in progress.

During the fist phases of the test, region I active grows up and I elastic shrinks, whereas I inactive is empty. Then, due to the evolution of interaction distances, points belonging to I active progressively move to I inactive , while I elastic stay constant (i.e., the size of the damaged band does not change).

To better understand such a damage/strain localization process, fig. 3 depicts the evolution of damage and local strain at x w and x w+1 , with respect to ¯ w . Since this quantity is controlled As already mentioned, despite what is provided by INL formulations, during the localization process, damage progressively tends to unity on the weak element (solid lines, fig. 3a) for both ENL models. Some differences between ENL1 and ENL2 formulations are, however, clear. In particular:

(i) the ENL1 formulation clearly models strain localization.

Damage is equal to unity in the central element, whereas it remains constant from a certain instant forward in the neighboring element (dashed blue line). Simultaneously, the non-local strain progressively tends to the local one in the weak element, while it gradually decreases in the neighboring element. As expected, such a process is accompanied by the elastic unloading of the neighboring element, which progressively tends to be undeformed.

It is worth noticing that, during the final damage localization phases, w decreases to ensure a gradual/smooth increase of ¯ w and thus of D w . This response clearly illustrates why a path-following constraint equation based on controlling ∆ w cannot properly work. Indeed, since w is not allowed to decrease in that case, ¯ w abruptly increases, leading to a rapid/uncontrolled increase in D w . This explains the jump (dashed blue line) observed in the global response depicted in fig. 1 (ii) conversely, when using the ENL2 formulation, damage in the neighboring elements continue to increase. Moreover, ¯ w never reaches w . This ensures more numerical robustness for the ENL2 formulation, as strain localization in only one element is significantly delayed.

This behavior can be easily explained in light of eq. ( 7), considering three neighboring integration points (x w , x w+1 , x w+2 ) such that D w+2 < D w+1 < D w and D w → 1 -(same considerations hold for (x w , x w-1 , x w-2 ) due to symmetry). Under these conditions, one can easily show that:

(i) in the case of the ENL1 formulation, condition D w → 1 -implies that ˜ w,w+1 can become larger than the maximum length such that φ( ˜ w,w+1 ) becomes very small (negligible):

˜ w,w+1 ≈ h 2 √ 1 -D w β c ( 8
)
with β ≈ 1.5 for the Gaussian weighting function. As a consequence, point x w+1 (as well as any other integration point) no more interacts with x w , and damage evolution cannot occur in the bar;

(ii) on the contrary, when D w → 1 -, the interaction distance ˜ w,w+1 that is considered in the ENL2 formulation can still be smaller than β c :

˜ w,w+1 ≈ 2h 1 -D w+1 ≤ β c (9)
if D w+1 is small enough (compared to unity). As a consequence, D w+1 can increase due to the increasing strain at x w (this term still contribute to the non-local averaged strain at x w+1 ).

From a theoretical point of view, the only way for preventing non-local interactions between x w and other points to occur, is that also D w+1 → 1 -. Indeed, since ˜ w+1,w+2 = 2h/( 1 -D w+1 + 1 -D w+2 ), one can show that:

˜ w,w+2 = ˜ w,w+1 + ˜ w+1,w+2 ≤ β c ( 10
)
as long as D w+1 and D w+2 are small enough. Once condition D w+1 → 1 -is achieved, ˜ w,w+1 and ˜ w+1,w+2 both become very large:

˜ w,w+1 β c (11) ˜ w+1,w+2 ≈ 2h 1 -D w+1 β c ( 12
)
and no more interactions are possible between points x w , x w+1 and their neighbors.

From a numerical viewpoint, such a condition is however attained when w → +∞ (i.e., for very high sample elongation levels). Notice that when w ≈ 1 = 10 4 × 0 , D w+1 is still close to 0.99 for the ENL2 formulation. Such a damage evolution can be limited numerically through considering different weighting functions (e.g., the bell-shaped polynomial function) such that φ = 0 for ˜ w,w+1 ≥ c . In that case, the maximum attainable damage value at x w+1 will depend on the mesh size (i.e., smaller h higher D w+1 ). Some additional considerations have to be provided concerning the fact that, when using the ENL1 formulation, ε w tends to reduce during the intermediate phases of the test, whereas εw is still increasing (imposed quantity). This response can be easily explained keeping in mind that the non-local strain variation at x w between t n and t n+1 , can be written as:

∆ε w = j (ε n j + ∆ε j )φ n w,j j φ n w,j - j ε n j φ n-1 w,j j φ n-1 w,j = ∆τ (13)
where ε n j = ε n (x j ) is the strain at x j and time t n and, given the chosen explicit procedure for updating interaction distances, φ n w,j and φ n-1 w,j are the weighting functions corresponding to damage distributions at t n and t n-1 . After some algebraic manipulations, eq. ( 13) can also be written as:

∆ε w = ∆ε * w + ∆ε * w + ε * w + ε * w = ∆τ (14)
with: The evolution of contributions appearing in eq. ( 14) is depicted in fig. 4. When damage occurs, integration points x j ∈ I that are far enough from x w start experiencing negative local strain variations. As a consequence, ∆ε * w start decreasing, while other contributions still increase in order to ensure the fulfillment of constraint [START_REF] Giry | Stress-based nonlocal damage model[END_REF]. During damage localization, the number of integration points experiencing negative strain variations increases, leading to ∆ε * w < 0. In the meanwhile, points x j = x w progressively move from I active to I inactive . During this phase, contribution ε * w + ε * w increases (it is mainly driven by ε w ) and become larger than ∆τ . Then, since ε w has become quite large, and negative contributions coming from neighbor points progressively vanish, fulfilling condition [START_REF] Giry | Stress-based nonlocal damage model[END_REF] implies that ∆ε * w < 0. The local strain at x w is thus reduced. Once the strain localization is complete, ∆ε w = ∆ε * w = ∆ε w = ∆τ > 0, whereas all the remaining contributions become negligible, because j φ n-1 w,j = j φ n w,j = 1 and j =w ∆ε j φ n w,j = j =w ε n j φ n w,j = j =w ε n j φ n-1 w,j ≈ 0.

∆ε * w = ∆ε w j φ n w,j (15)
∆ε * w = j =w ∆ε j φ n w,j j φ n w,j (16)
ε * w = ε n w j φ n-1 w,j -j φ n w,j ( j φ n w,j )( j φ n-1 w,j ) (17)
ε * w = j =w ε n j φ n w,j j φ n w,j - j =w ε n j φ n-1 w,j j φ n-1 w,j (18)
INL -ε w INL -ε w+1 ENL1 -ε w ENL1 -ε w+1 ENL2 -ε w ENL2 -ε w+1

Parasite/residual mesh-dependency

Features mentioned above make ENL formulations suitable for representing localization processes characterized by a progressive transition from diffuse damage to fracture (e.g., cracking of quasi-brittle materials). The observed tendency toward strain localization is responsible, however, of the small differences that are observed in the global structural responses computed for different finite element meshes for high damage levels (fig. 5).

For both ENL formulations, when non-local interactions between the weak element and its neighbors vanish (or are lower than a critical value), damage evolution becomes local. The energy still available at the weak element is consumed in a local manner [START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF]. In other words, regularization properties are lost, and parasite mesh dependency issues are observed.

As a consequence, in addition to inducing damage diffusion for high elongation levels, the ENL2 formulation only slightly helps in improving such a feature/limitation which is proper to other non-local damage models with evolving interaction distances (see e.g., [START_REF] Giry | Stress-based nonlocal damage model[END_REF]).

Conclusive remarks

ENL damage formulations exhibit good regularization properties and naturally allow one to model a progressive strain localization and damage-to-fracture transition [START_REF] Desmorat | Nonlocal models with damagedependent interactions motivated by internal time[END_REF][START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF][START_REF] Jirásek | Localization analysis of nonlocal models with damage-dependent nonlocal interaction[END_REF].

Focusing on a 1D formulation, in this contribution, we analyzed the effect of different distances integration techniques to compute interaction distances between points from solving a damage-dependent Eikonal problem. In particular, we pointed out that parasite damage diffusion can occur when assuming that the damage field is element-wise linear in the metric function [START_REF] Jirásek | Localization analysis of nonlocal models with damage-dependent nonlocal interaction[END_REF]. Conversely, this is not the case when considering that damage is element-wise constant [START_REF] Desmorat | Nonlocal models with damagedependent interactions motivated by internal time[END_REF][START_REF] Rastiello | From diffuse damage to strain localization from an Eikonal Non-Local (ENL) Continuum Damage model with evolving internal length[END_REF].

Moreover, controlling the simulation through the local strain variation of the element where damage localization occurs, does not allow to simulate a smooth damage evolution. Conversely, this result can be achieved by using a constraint equation based on the non-local strain variation. Once this choice is made, we showed that the localizing element's local strain has to decrease during the final localization phases to ensure a gradual evolution of the non-local strain and, thus, a smooth damage evolution.

Appendix A. Solving method

A path-following method is used to solve the non-linear equilibrium problem. Accordingly, the external force vector is written as f ext = λ f , where the loading factor λ ∈ R is a novel problem unknown and f is a constant vector providing the direction of the loading.

The problem to be solved consists of computing the nodal displacements vector (d) and the loading factor (λ) that fulfill the equilibrium equation, r(d, λ;

{D i } n el i=1 ) = f int (d; {D i } n el i=1
)f ext (λ) = 0 (with r denoting the residual and f int being the internal forces vector), and a path-following constraint equation, p(d, λ; ∆τ ) = 0.

The unknown nodal displacement variation vector δd k+1 at the k + 1-th iteration is decomposed as [START_REF] Riks | An incremental approach to the solution of snapping and buckling problems[END_REF][START_REF] Crisfield | A fast incremental/iterative solution procedure that handles "snap-through[END_REF][START_REF] De Borst | Nonlinear Finite Element Analysis of Solids and Structures[END_REF]:

δd k+1 = δλ k+1 δd k+1 I + δd k+1 II (A.1)
where δd k+1 I and δd k+1 II are solved from:

K k δd k+1 I = f K k δd k+1 II = -r k (A.2) with K k = K k (d k ; {D k i } n el i=1 ) denoting the global secant stiff- ness matrix at iteration k and r k = r k (d k , λ k ; {D k i } n el i=1 ) the residual vector.
The resolution is performed in two steps. Displacement variations (δd k+1 I , δd k+1 II ) are first computed. Then the loading factor increment δλ k+1 is computed from solving the path-following constraint equation. Details and examples of path-following constraint equations can be found in existing papers (see e.g., [START_REF] Chen | A numerical solution scheme for softening problems involving total strain control[END_REF][START_REF] Lorentz | A new path-following constraint for strainsoftening finite element simulations[END_REF][START_REF] Pohl | Adaptive path following schemes for problems with softening[END_REF][START_REF] Rastiello | Discontinuity-scale pathfollowing methods for the embedded discontinuity finite element modeling of failure in solids[END_REF] among others).

To control the non-local strain variation of the weak element, the constraint equation is written as p = ∆ε k+1 w -∆τ = 0, where ∆ε k+1 w reads as in eq. ( 13). Given the additive displacement decomposition (A.1) and the linear relationship between local strains and nodal displacements (through strain -displacement interpolation matrices), previous constraint equation can be easily rewritten as a linear function of δλ k+1 . One finally obtains: δλ k+1 = ∆τ * j φ n w,jj ∆ε k j φ n w,jj δε k+1 j,II φ n w,j j δε k+1 j,I φ n w,j (A.3) where ∆ε k j is the local strain variation between iteration k and loading step n, (δ εk+1 j,I , δ εk+1 j,II ) are the non-local strain variations corresponding to displacement variations (A.2), and:

∆τ * = ∆τ - j n j φ n w,j j φ n w,j - j n j φ n-1 w,j j φ n-1 w,j (A.4)
Notice that, when using an INL formulation, ∆τ * = ∆τ , and φ n w,j = φ n-1 w,j = φ w,j , with φ w,j being the constant value of the weighting function between points x w and x j .

Solving the constraint equation on the strain variation of the weak element, p = ∆ε k+1 w -∆τ = 0, is even more simple. In that case, one directly obtains: Once δλ k+1 is computed, d k+1 and λ k+1 are updated and the discretized damage field is evaluated according to an explicit procedure. The residual vector is then computed, and convergence checked based on the following error measure, err k+1 = ||r k+1 ||/(Sf t ), with f t denoting the tensile strength of the material.

During computations the stiffness matrix has been updated every 5 global iterations, the Choleski decomposition was used for solving (A.2), and the tolerance was set to 10 -5 . All the results depicted in the paper were obtained considering ∆τ = κ 0 /50 (for both the constraint equations) in order to ensure a gradual damaging process.

One can show that the constraint equations were strictly fulfilled at each pseudo-time step for all the simulations performed in this study. However, convergence properties are quite different depending on the chosen ENL formulation. In both cases, if one refers to results depicted in fig. 1, the number of global iterations to converge (m) is lower for the constraint equation on ∆¯ w (m ∈ [START_REF] Bažant | Nonlocal Integral Formulations of Plasticity and Damage: Survey of Progress[END_REF][START_REF] Ganghoffer | A reappraisal of nonlocal mechanics[END_REF] for both the ENL formulations) than for the constraint equation on ∆ w (m ∈ [START_REF] Bažant | Nonlocal Integral Formulations of Plasticity and Damage: Survey of Progress[END_REF]50] for the ENL1 formulation, and m ∈ [START_REF] Bažant | Nonlocal Integral Formulations of Plasticity and Damage: Survey of Progress[END_REF][START_REF] Krayani | Boundary effect on weight function in nonlocal damage model[END_REF] for the ENL2 model). Notice that for the mesh considered in fig. 1 (n el = 81), convergence is always achieved, but for finer meshes (e.g., for n el = 161) controlling the simulation through ∆ w does not ensure achieving convergence when using the ENL1 formulation. Now, if one refers to results depicted in fig. 5, one can easily show that for the chosen ∆τ , the number of iteration to converge (at a given non-local strain level) is almost independent of the mesh refinement (for both ENL formulations). Conversely, as expected, there is a tendency to increase the number of iterations when ∆τ increases, but such an increase is quite limited and does not strongly affect the convergence features of the solving algorithm.

Figure 1 :

 1 Figure 1: ENL formulations 1 and 2 -Force-displacement responses obtained using local, INL and ENL damage models (finite element mesh with 81 elements). Results obtained controlling the simulation on the non-local strain variation of the weak element and on the local strain variation of the weak element -Numerical parameters: L = 100 mm, S = 5•10 -3 m 2 , E = 100 MPa, κ 0 = 10 -4 , κ c = 10 κ 0 , c = L/5 = 20 mm, ∆τ = κ 0 /50 (for both the constraint equations)

Figure 2 :

 2 Figure 2: ENL formulations 1 and 2 -Damage distributions along the bar at different time-steps (81 finite elements).

2 n

 2 Influence of the integration rule for computing interaction distances and indirect loading control on damage localization x w+1 n el = 81 Non-local strain at x w , εw Damage, D INL -D w INL -D w+1 ENL1 -D w ENL1 -D w+1 ENL2 -D w ENL2 -D w+1 el = 81 Non-local strain at x w , εw Local strain, ε

Figure 3 :

 3 Figure 3: ENL formulations -Damage (a) and non-local strain (b) evolution at x w and x w+1 (81 finite elements).

6 nFigure 4 :

 64 Figure 4: ENL formulation 1 -Contributions to the non-local strain variation (14) at the weak element (finite element mesh comprising 81 elements).

Figure 5 :

 5 Figure 5: ENL formulations -Force-displacement (F, u L ) responses obtained for different mesh refinement levels (41, 81, 161 and 2001 finite elements) considering the ENL1 (a) and ENL2 (b) formulations. The corresponding element size values range from c /8 to c /400 approximately.

Acknowledgements

G. Rastiello was supported by the SEISM Institute (http: //www.institut-seism.fr). The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


D Perconte

K Seurre

V Humbert

C Ulysse

A Sander

J Trastoy

V Zatko

F Godel

P R Kidambi

S Hofmann

X P Zhang

D Bercioux

F S Bergeret

B Dlubak

P Seneor

Javier E Villegas
email: javier.villegas@cnrs-thales.fr

Long-Range Propagation and Interference of d-wave Superconducting Pairs in Graphene

Recent experiments have shown that proximity with high-temperature superconductors induces unconventional superconducting correlations in graphene. Here we demonstrate that those correlations propagate hundreds of nanometer, allowing for the unique observation of d-wave Andreev pair interferences in YBa2Cu3O7-graphene devices that behave as a Fabry-Pérot cavity. The interferences show as a series of pronounced conductance oscillations analogous to those originally predicted by de Gennes-Saint-James for conventional metal-superconductor junctions. The present work is pivotal to the study of exotic directional effects expected for nodal superconductivity in Dirac materials.

The superconducting proximity effect in graphene has attracted much interest since the pioneering experiments [START_REF] Linder | Dirac fermions and conductance oscillations in s-and d-wave superconductor-graphene junctions[END_REF]. This roots down to the graphene's electronic structure, which strongly affects the underlying mechanisms: Andreev reflection and coherent propagation of electron-hole pairs [START_REF] Wei | Directional Tunneling and Andreev Reflection on YBa2Cu3O 7-δ Single Crystals: Predominance of d -Wave Pairing Symmetry Verified with the Generalized Blonder, Tinkham, and Klapwijk Theory[END_REF]. A distinctive feature is the strong dependence of the proximity behavior on the graphene?s doping level, which dramatically changes the Andreev reflection [3,4] as compared to metals. Other unique features include gate [5] or magnetic-field [6,7] driven transitions from bulk to edge transport. Studies on graphene have also paved the way for understanding the proximity effect in other Dirac materials, such as topological insulators [START_REF] Koren | Signature of proximity-induced px + ipy triplet pairing in the doped topological insulator Bi2Se3 by the s-wave superconductor NbN[END_REF][START_REF] Finck | Phase coherence and andreev reflection in topological insulator devices[END_REF][START_REF] Finck | Robust Fabry-Pérot interference in dual-gated Bi2Se3 devices[END_REF][START_REF] Bocquillon | Gapless andreev bound states in the quantum spin hall insulator HgTe[END_REF].

Experiments have mostly focused on conventional low critical temperature (T C ) superconductors with swave pairing. Despite early theoretical studies showing that d-wave (high-T C ) superconductors should lead to novel directional effects [START_REF] Linder | Dirac fermions and conductance oscillations in s-and d-wave superconductor-graphene junctions[END_REF][START_REF] Linder | Tunneling conductance in s-and d-wave superconductor-graphene junctions: Extended Blonder-Tinkham-Klapwijk formalism[END_REF] and exotic pairing [START_REF] Linder | Josephson current in graphene: Role of unconventional pairing symmetries[END_REF][START_REF] Oladunjoye | Probing unconventional superconductivity in proximitized graphene by impurity scattering[END_REF], evidence for unconventional superconductivity in graphene has been found only recently [START_REF] Bernardo | p-wave triggered superconductivity in single-layer graphene on an electron-doped oxide superconductor[END_REF][START_REF] Perconte | Tunable Klein-like tunnelling of high-temperature superconducting pairs into graphene[END_REF]. Scanning electron tunneling microscopy (STM) of graphene on Pr 2-x Ce x CuO 4 (PCCO) revealed a superconducting gap [START_REF] Bernardo | p-wave triggered superconductivity in single-layer graphene on an electron-doped oxide superconductor[END_REF] and spectral features suggesting (p-wave) superconductivity induced in graphene. Experiments on YBa 2 Cu 3 O 7 (YBCO) and chemical-vapor-deposited (CVD) graphene devices by some of us [START_REF] Perconte | Tunable Klein-like tunnelling of high-temperature superconducting pairs into graphene[END_REF] showed trans-FIG. 1. Interferences in a proximitized cavity (a) An electron Andreev-reflected at the SC/cavity interface propagates as a hole back to the other end, where it is normal-reflected towards the SC/cavity interface to undergo again Andreevreflection (AR). The lower AR is the time reversal process of the upper one. (b) The electron can also be normal-reflected at the SC/cavity, travel back to the other cavity interface to be again normal-reflected. This process results in Fabry-Pérot resonances.

parent superconductor-graphene interfaces and clear evidence of Andreev reflection. Interestingly, we also found that the Andreev electron-holes pair transmission can be modulated by a back-gate voltage, through a mechanism analogous to the Klein tunneling [START_REF] Young | Quantum interference and Klein tunnelling in graphene heterojunctions[END_REF]. However, neither the STM experiments nor ours on solid-state devices probed the length scale over which the unconventional correlations penetrate into graphene.

Here we demonstrate the long-range propagation of unconventional superconductivity into graphene via the observation of d-wave Andreev-pair interferences. These manifests themselves as conductance oscillations in devices that allow confining Andreev pairs within a graphene "cavity" whose length L is up to a few hundreds of nanometers. Predicted for proximitized graphene homojunctions [START_REF] Linder | Dirac fermions and conductance oscillations in s-and d-wave superconductor-graphene junctions[END_REF][START_REF] Linder | Tunneling conductance in s-and d-wave superconductor-graphene junctions: Extended Blonder-Tinkham-Klapwijk formalism[END_REF], the oscillations are analogous to the de Gennes-Saint-James [START_REF] De Gennes | Elementary excitations in the vicinity of a normal metal-superconducting metal contact[END_REF] and McMillan-Rowel resonances [START_REF] Rowell | Electron interference in a normal metal induced by superconducting contracts[END_REF] in the electronic density of states of ultrathin normal-metals backed by superconductors. Figure 1(a) displays a cartoon of the underlying mechanism. Electrons injected in the cavity are Andreev-reflected as holes at the interface with the superconductor, retrace their path to the opposite cavity's end where they are normalreflected, and travel bain ck towards the superconductor, where they are again Andreev-reflected (now as electrons). This results in destructive/constructive interferences dictated by the energy-dependent phase accumulated along the loop. This phenomenon shows up as conductance oscillations as a function of the bias voltage V BIAS . Because Andreev pairs stemming from YBCO have d-wave symmetry, they decay over distances comparable to the mean free path l [START_REF] Tsuneto | On dirty superconductors[END_REF][START_REF] Balatsky | Impurityinduced states in conventional and unconventional superconductors[END_REF], the observation of Andreev-pair interferences implies ballistic or quasiballistic transport. Consistently, the Andreev-pair interferences are accompanied by normal-electron resonances that result from normal-reflections at the two cavity's ends [sketch in Fig. 1(b)] and commensurability between L and the electrons' wavelength [START_REF] Young | Quantum interference and Klein tunnelling in graphene heterojunctions[END_REF][START_REF] Liang | Fabry-Pérot interference in a nanotube electron waveguide[END_REF][START_REF] Miao | Phase-coherent transport in graphene quantum billiards[END_REF][START_REF] Campos | Quantum and classical confinement of resonant states in a trilayer graphene Fabry-Pérot interferometer[END_REF][START_REF] Allen | Observation of Electron Coherence and Fabry-Pérot Standing Waves at a Graphene Edge[END_REF]. These different resonances result in oscillations as a function of both the gate V G and bias V BIAS , and can be distinguished from Andreev-pairs related ones by their distinct periodicity and voltage regime in which they dominate.

To our knowledge, the concurrent observation of both resonances had been restricted to junctions between swave superconductors and topological insulators [START_REF] Finck | Phase coherence and andreev reflection in topological insulator devices[END_REF][START_REF] Finck | Robust Fabry-Pérot interference in dual-gated Bi2Se3 devices[END_REF], in which Andreev-pair interferences manifest at very lowtemperatures as a function of the gate voltage V G . Here they can be observed at higher temperatures and show in a rich series of oscillations as a function of V BIAS because of the high-T C and large gap of YBCO (tens of meV) which, contrary to low-temperature superconductors [START_REF] Giazotto | Resonant Transport in Nb /GaAs /AlGaAs Heterostructures: Realization of the de Gennes-Saint-James Model[END_REF], may enclose many orders of interference. The planar devices [Fig. 2(a)] consist of four superconducting (SC) c-axis oriented YBCO 50nm /Au 4nm electrodes disposed within an insulating (INS) YBCO matrix and separated by a gap of length 100 nm< L < 800 nm. As-grown YBCO films showed T C ∼ 90 K while T C ∼ 70 K for lithographed devices, indicating moderate deoxygenation during fabrication. The ultrathin Au interlayer is deposited in situ to preserve the superconducting properties and improve the interfacial transparency [START_REF] Perconte | Tunable Klein-like tunnelling of high-temperature superconducting pairs into graphene[END_REF]. The Au interlayer is crucial: none of the effects reported here could be observed in its absence. Notice that the Au layer thickness it is one order of magnitude thinner than the mean free path in Au, and conse- quently d-wave correlations expectedly propagate undisturbed across it [START_REF] Tsuneto | On dirty superconductors[END_REF][START_REF] Balatsky | Impurityinduced states in conventional and unconventional superconductors[END_REF]. A single-layer CVD grown [START_REF] Piran | The Parameter Space of Graphene Chemical Vapor Deposition on Polycrystalline Cu[END_REF] graphene bridge (G) connects the four electrodes, each of which terminates in a gold pad [1-4 in Fig. 2(b)] used for wire-bonding. The transparency of the electrical contact between each SC electrode and the graphene varied largely, even within a single device. Here we report on effects that arise only at high transparency contacts, which is often the case for just one out of four contacts per device. The structure is covered by a 45 nm thick AlO x layer (gate dielectric) [START_REF] Dlubak | Substrate-assisted nucleation of ultra-thin dielectric layers on graphene by atomic layer deposition[END_REF]. As reported earlier [START_REF] Mzali | Stabilizing a graphene platform toward discrete components[END_REF], with this encapsulation the carrier mobility in the used CVD graphene reaches up to µ ∼ 7000 cm 2 V -1 s -1 . The top-gate is made of Au entirely [Fig. 2(c)] and the gate capacitance as measured over insulating YBCO is C G ∼ 5 × 10 11 e cm -2 V -1 (with e the electron charge). Further details are reported elsewhere [START_REF] Perconte | Tunable Klein-like tunnelling of high-temperature superconducting pairs into graphene[END_REF].

We performed three-probe measurements, with the current I injected e.g. from 1 to contact 4 and the V BIAS measured between 1 and 3. Thus we probe the YBCO-Au-graphene interface in series with the graphene "cavity" formed between the SC electrodes, as sketched in in Fig. 2(d). The wiring as well as Au/YBCO contact resistances are Ohmic and negligible [START_REF] Perconte | Tunable Klein-like tunnelling of high-temperature superconducting pairs into graphene[END_REF]. The differential conductance G(V BIAS ) = dI/dV BIAS is measured using a Keithley 6221 current source coupled to a Keithley 2182.

Examples of low-temperature (3.2 K) differential conductance measurements are shown in Fig. 3(a), 3(d) and Fig. 4(a). The differential conductance G(V BIAS ) shows a low-bias feature that stands out from the high-bias conductance level. That is either a conductance decrease [START_REF] Blonder | Transition from metallic to tunneling regimes in supercon-ducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion[END_REF]. However, here they extend over a bias range δ exp that clearly exceeds the superconducting energy-gap expected for moderately-underdoped to optimally-doped YBCO, which is reported in the range ∆ YBCO ∼ 15-30 meV [START_REF] Maggio-Aprile | Direct Vortex Lattice Imaging and Tunneling Spectroscopy of Flux Lines on YBa2Cu3O 7-δ[END_REF][START_REF] Wei | Directional Tunneling and Andreev Reflection on YBa2Cu3O 7-δ Single Crystals: Predominance of d -Wave Pairing Symmetry Verified with the Generalized Blonder, Tinkham, and Klapwijk Theory[END_REF][START_REF] Dagan | Determination of the superconducting gap in YBa2Cu3O 7-δ by tunneling experiments under magnetic fields[END_REF][START_REF] Rouco | Quasiparticle tunnel electroresistance in superconducting junctions[END_REF]. Furthermore, details such as the sharp zero-bias peak observed for some devices (see e.g. B4D and E3D in the Supplemental Material [START_REF]ulation parameters, V. Table of simulation parameters, VI. Simulation vs. experiment for the set of devices, VII. Induced superconducting correlations: d-wave versus swave[END_REF]) are characteristic of junctions involving d-wave superconductors, and appear when the d-wave nodes form a certain angle with respect to the junction interface [START_REF] Wei | Directional Tunneling and Andreev Reflection on YBa2Cu3O 7-δ Single Crystals: Predominance of d -Wave Pairing Symmetry Verified with the Generalized Blonder, Tinkham, and Klapwijk Theory[END_REF][START_REF] Kashiwaya | Origin of zerobias conductance peaks in high-Tc superconductors[END_REF][START_REF] Sharoni | Proximity Effect in Gold-Coated YBa2Cu3O 7-δ Films Studied by Scanning Tunneling Spectroscopy[END_REF].

In addition to the central features, an oscillation pattern appears in Fig. 3(a) and 3(d), this extends over a bias range well above δ exp . Indeed, two types of oscillations exist: the ones with longer period V long , predominant for |V BIAS | > δ exp , and the other with shorter V short which are more prominent at |V BIAS | < δ exp . Figure 4 show an example of gating effects. Figure 4(a) displays G(V BIAS ) without gate voltage (V G = 0), and shows the general features discussed above. Figure 4(b) is a color plot of a series of G(V BIAS ) at constant temperature and varying V G . The central (red) feature cor-responds to the zero-bias conductance peak. The conductance is periodically modulated by V G : a pattern of oblique lines (light-blue/green) indicates that V G gradually "shifts" the oscillations observed as a function of V BIAS . The oblique lines' slope gradually varies with V G , which results in pronounced curvature over the plot periphery (low V G or high |V BIAS |).

We discuss now a model that explains the main observations, namely (i) the increase or decrease of the zerobias conductance; (ii) the origin and period of conductance oscillations. The model is based on the Blonder-Tinkham-Klapwijk (BTK) formalism [START_REF] Blonder | Transition from metallic to tunneling regimes in supercon-ducting microconstrictions: Excess current, charge imbalance, and supercurrent conversion[END_REF] extended to junctions between d-wave superconductors and normal metals [START_REF] Kashiwaya | Theory for tunneling spectroscopy of anisotropic superconductors[END_REF] and to proximitized graphene homojunctions [START_REF] Linder | Dirac fermions and conductance oscillations in s-and d-wave superconductor-graphene junctions[END_REF][START_REF] Linder | Tunneling conductance in s-and d-wave superconductor-graphene junctions: Extended Blonder-Tinkham-Klapwijk formalism[END_REF]. A scheme is shown in Fig. 2(d). First, we consider the YBCO-Au-graphene interface. Because the Au thickness (∼ 5 nm) is well-below the mean free path l Au ∼ 40 nm [START_REF] Scheer | Proximity effect and multiple an-dreev reflections in gold atomic contacts[END_REF] and the low-T coherence length ξ Au = v F l/(6πKT ) ∼ 30 nm [START_REF] Sharoni | Proximity Effect in Gold-Coated YBa2Cu3O 7-δ Films Studied by Scanning Tunneling Spectroscopy[END_REF], we characterize that interface [black in Fig. 2(d)] via a single Blonder-Tinkham-Klapwijk barrier-strength parameter Z. If Z = 0, the transmission between YBCO and graphene is only mediated by the Andreev reflection. By increasing Z, Andreev reflection turns less dominant, normal reflection is enhanced, and transmission is dominated by tunneling. We model the graphene "channel" where Andreev pairs and normal electrons propagate, considering three different regions A, B, C. A and C correspond to graphene lying on superconducting YBCO/Au [yellow in Fig. 2(d)], but we consider that only contact A has a relatively low Z. The region B of length L [dark in Fig. 2(d)] lies on insulating YBCO. Here the Fermi energy E F,B is expectedly different from E F,A , E F,C as the graphene's doping depends on the substrate's electronic properties [START_REF] Hwang | Fermi velocity engineering in graphene by substrate modification[END_REF]. Within this model, the structure is a graphene homojunction with a gate-tunable Fermi energy step, as sketched in the upper Fig. 2(d), and behaves as a resonant cavity for electrons and Andreev pairs. The conductance is calculated numerically by considering the YBCO-Au-graphene interface [START_REF] Wei | Directional Tunneling and Andreev Reflection on YBa2Cu3O 7-δ Single Crystals: Predominance of d -Wave Pairing Symmetry Verified with the Generalized Blonder, Tinkham, and Klapwijk Theory[END_REF][START_REF] Kashiwaya | Theory for tunneling spectroscopy of anisotropic superconductors[END_REF] and the proximitized graphene homojunction [START_REF] Linder | Dirac fermions and conductance oscillations in s-and d-wave superconductor-graphene junctions[END_REF][START_REF] Linder | Tunneling conductance in s-and d-wave superconductor-graphene junctions: Extended Blonder-Tinkham-Klapwijk formalism[END_REF] in-series (see Supplemental Material for details [START_REF]ulation parameters, V. Table of simulation parameters, VI. Simulation vs. experiment for the set of devices, VII. Induced superconducting correlations: d-wave versus swave[END_REF]). The input parameters are Z, the reduced Fermi energies F,i ≡ E F,i /∆ (with i = A, B, C and ∆ the energy-gap induced in A by proximity effect), the reduced cavity's length Λ ≡ L/λ F,A (with the Fermi wavelength λ F,A , chosen as a reference following the original model [START_REF] Linder | Dirac fermions and conductance oscillations in s-and d-wave superconductor-graphene junctions[END_REF][START_REF] Linder | Tunneling conductance in s-and d-wave superconductor-graphene junctions: Extended Blonder-Tinkham-Klapwijk formalism[END_REF]), and α the effective angle between d-wave anti-nodes and the homojunction interfaces. The key parameters are Λ and F,A : they determine V long and V short . The parameters α and Z determine the conductance background shape, and particularly whether there is a conductance decrease or enhancement around zero bias. We stress that the sharp zero-bias peak observed in some devices cannot be reproduced by the simulation unless we consider α = 0 (see Supplemental Material [START_REF]ulation parameters, V. Table of simulation parameters, VI. Simulation vs. experiment for the set of devices, VII. Induced superconducting correlations: d-wave versus swave[END_REF] Sec. 7). This is a consequence of the d-wave symmetry of the correlations propagating across the different interfaces. The parameters F,B and F,C only change the oscillations' phase and amplitude, but do not affect their period [START_REF] Linder | Tunneling conductance in s-and d-wave superconductor-graphene junctions: Extended Blonder-Tinkham-Klapwijk formalism[END_REF]. To choose the sim-ulation parameters that yield the best agreement with the experimental curves we proceed as detailed in the Supplemental Material in Sec. 2-4 [START_REF]ulation parameters, V. Table of simulation parameters, VI. Simulation vs. experiment for the set of devices, VII. Induced superconducting correlations: d-wave versus swave[END_REF]. Simulations examples are shown in Fig. 3 In order to illustrate the different contributions to the conductance, Figs. 3(c) and 3(f) display the conductance of a proximitized graphene homojunction [START_REF] Linder | Tunneling conductance in s-and d-wave superconductor-graphene junctions: Extended Blonder-Tinkham-Klapwijk formalism[END_REF] alone, i.e. without a finite-Z junction in series. The same F,i , Λ and α as for the full-model calculations are used.

Here the conductance shows a background dependence on V BIAS very different from the experiments, which evidences that it is strongly influenced by the YBCO-Augraphene interface. In particular, the zero-bias "dip" only emerges if a finite-Z junction is considered. However, Figs. 3(c) and 3(f) do show the short and long period oscillations observed experimentally. This demonstrates that they originate within the graphene homojunction. Notice that the ratio between the short and long periods is different from that in the experiments (and full model simulations): in Fig. 3(c) and 3(f), V short is clearly shorter (relative to V long ) than in Fig. 3(a)-3(b) and 3(d)-3(e). This is because in experiments (and full-model) V BIAS is divided between the graphene homojunction and the YBCO-Au-graphene interface, and the strongly non-linear conductance leads to a voltage distribution that varies depending on the V BIAS . In particular, for eV BIAS < ∆ the conductance across the YBCO-Augraphene interface decreases (this is more pronounced for higher Z values) while the homojunction conductance increases [see Figs. 3(c), 3(f)]. Consequently, V short (observed for eV BIAS < ∆) appears "stretched" relative to the V long (observed for eV BIAS > ∆).

The oscillations' physical meaning can be understood from Fig. 5. This displays the short (red circles) and long (black squares) oscillation periods as a function of the inverse of the cavity length L -1 device defined upon devices fabrication. Specifically, we plot the period V = v∆ , where ∆ is estimated as discussed in the Supplemental Material [START_REF]ulation parameters, V. Table of simulation parameters, VI. Simulation vs. experiment for the set of devices, VII. Induced superconducting correlations: d-wave versus swave[END_REF] and v is the period obtained from our model once the contribution of the YBCO-Au-graphene interface has been removed -as in Fig. 3(c) and 3(e) -to avoid the aforementioned V BIAS division artifacts.

The long-period oscillations follow V th, long = hv F /2L device . This is as expected from the interference between electrons travelling back and forth from one cavity side to the other after normal reflections [Fig. 1(b)]. This period results [5,[START_REF] Young | Quantum interference and Klein tunnelling in graphene heterojunctions[END_REF][START_REF] Liang | Fabry-Pérot interference in a nanotube electron waveguide[END_REF][START_REF] Miao | Phase-coherent transport in graphene quantum billiards[END_REF][START_REF] Campos | Quantum and classical confinement of resonant states in a trilayer graphene Fabry-Pérot interferometer[END_REF] from the interference condition 2L device k = 2nπ (with n an integer and k the electron wavevector) and the graphene's linear dispersion, which yields V = v F k [START_REF] Castro Neto | The electronic properties of graphene[END_REF]. These interferences imply normal-electron coherence over ∼ L device . Their )

FIG. 5.
Oscillations period as a function of the cavity length L device . The black and red dashed lines respectively correspond to the period theoretically expected for electron (Fabry-Pérot) and Andreev-pair (De Gennes-Saint James) interferences.

amplitude is greater for 100 nm < L device < 300 nm but drastically diminishes beyond that. This suggests that the mean free path is above 100 nm but clearly below 1 µm, in agreement with a rough estimate based on the carrier mobility [START_REF] Hwang | Acoustic phonon scattering limited carrier mobility in two-dimensional extrinsic graphene[END_REF] which yields l 100 -260 nm . Notice that the scaling of the period with 1/L device rules out Tomasch resonances [START_REF] Tomasch | Geometrical resonance and boundary effects in tunneling from superconducting in[END_REF][START_REF] Nesher | Observation of tomasch oscillations and tunneling-like behavior in oxygen-deficient edge junctions[END_REF][START_REF] Visani | Equal-spin andreev reflection and long-range coherent transport in hightemperature superconductor/half-metallic ferromagnet junctions[END_REF] within YBCO, which should scale with its thickness and should therefore show the same period for all devices.

The short period oscillations, clearly visible for V < ∆, follow V th,short = hv F /4L device . This period is expected from the interference of Andreev pairs [START_REF] De Gennes | Elementary excitations in the vicinity of a normal metal-superconducting metal contact[END_REF][START_REF] Rowell | Electron interference in a normal metal induced by superconducting contracts[END_REF] -scheme in Fig. 1(a). In this case the resonance condition reads 2L device |k 1k 2 | = 2nπ, where k 1 and k 2 are the electron and reflected hole wavevectors, whose difference is established by the applied voltage |k 1k 2 | = 2V / v F [START_REF] Rowell | Electron interference in a normal metal induced by superconducting contracts[END_REF]. These resonances imply that the superconducting coherence is preserved over ∼ L device .

Our model also explains the gating effects. Figure 4(d) displays a simulated conductance as a function of bias and Fermi energy, which reproduces the main experimental features -particularly the pattern of oblique lines that characterize the conductance modulation by by V G [Fig. 4(b)]. The correspondence between the y-axis in simulations and experiments implies that the Fermi energy (and vector k F ) is varied about proportionally to V G , in agreement with earlier experiments [START_REF] Perconte | Tunable Klein-like tunnelling of high-temperature superconducting pairs into graphene[END_REF]. Varying k F produces a periodic modulation of the conductance, via the resonance condition 2L device k F = 2nπ, yielding the oblique lines. Note that the model also reproduces the slope gradual change, and the pronounced curvature over the plot's periphery, which is obtained by including in the model F,B (V G ) = F,A (V G ) = F,C (V G ) to account for a partial pining of graphene Fermi energy on conducting YBCO.

In summary, we have realized ballistic planar devices in which graphene is proximitized by a d-wave cuprate superconductor. We find that the d-wave correlations propagate into graphene over hundreds of nanometers at ∼ 4 K. The confinement of Andreev pairs and electrons in graphene homojunctions of submicrometric lateral dimensions produces quantum interferences that show as conductance oscillations. Those involving Andreev pairs, more prominent below the superconducting gap, are analogous to the De Gennes-Saint James oscillations. Although difficult to observe in graphene/low-temperature superconductor devices, here they are very clear due to the large YBCO gap (tens of mV), which fits many orders of interference within. The unusual d-wave Andreev-pair interferences are accompanied by normal-electron ones. The simultaneous observation of both is rare, and stems from the fact that the coherence length of d-wave correlations is limited by the mean free path [START_REF] Balatsky | Impurityinduced states in conventional and unconventional superconductors[END_REF]. Thus, the necessary condition for both types of resonances is satisfied simultaneously here. Extensions of this study to graphene with lower intrinsic doping and higher carrier mobility should allow exploring specular Andreev reflection regimes [3] rand Josephson effects. The latter would also require enhancing the transparency of the graphene/YBCO interface, possibly by replacing the Au interlayer, in order to routinely obtain low-Z junctions. Further, perspectives include novel directional effects [START_REF] Linder | Josephson current in graphene: Role of unconventional pairing symmetries[END_REF][START_REF] Oladunjoye | Probing unconventional superconductivity in proximitized graphene by impurity scattering[END_REF] linked to the d-wave character of the induced correlations, which could be exploited in more sophisticated devices to realize topological states [START_REF] Lucignano | Topological rf squid with a frustrating π junction for probing the majorana bound state[END_REF]. This work should also further encourage studies in which graphene is replaced by other 2D materials with intrinsically interesting proximity behavior [START_REF] Kuzmanovski | Quantum ground state control in superconductor-silicene structures: 0-π transitions, ϕ0junctions, and majorana bound states[END_REF][START_REF] Linder | Anisotropic andreev reflection and josephson effect in ballistic phosphorene[END_REF], eventually combined with exfoliated d-wave superconductors [START_REF] Li | Bi-2212/1T-TaS2 Van der Waals junctions: Interplay of proximity induced high-T c superconductivity and CDW order[END_REF].

Work at Unité Mixte de Physique CNRS/Thales supported by the ERC Grant No. 647100 "SUSPINTRON-ICS", French ANR Grants No. ANR-15-CE24-0008-01 "SUPERTRONICS" and No. ANR-17-CE30-0018-04 "OPTOFLUXONICS", Labex NanoSaclay No. ANR-10-LABX-0035, European COST action 16218 "Nanocohybri" , and EU Work Programme Graphene Flagship (Grants No. 696656 and 785219). This work was supported by the French RENATECH network (French national nanofabrication platform). X.P.Z., D.B., and F.S.B. acknowledge support from Ministerio de Ciencia e Innovacion (MICINN) through Project No. FIS2017-82804-P. F. S. B. acknowledges Grupos Consolidados UP-V/EHU del Gobierno Vasco (Grant No. IT1249-19). We thank Professor. J. Santamaría for discussions.

I. THEORETICAL MODEL

We reproduce the experimental device conductance via numerical simulations based on combination of the models of Refs. [START_REF] Linder | Dirac fermions and conductance oscillations in s-and d-wave superconductor-graphene junctions[END_REF] and [START_REF] Wei | Directional Tunneling and Andreev Reflection on YBa2Cu3O 7-δ Single Crystals: Predominance of d -Wave Pairing Symmetry Verified with the Generalized Blonder, Tinkham, and Klapwijk Theory[END_REF]. The model of Ref. [START_REF] Linder | Dirac fermions and conductance oscillations in s-and d-wave superconductor-graphene junctions[END_REF] applies to a graphene homojunction proximitized by a d-wave superconductor.

The junction is divided into three different sections A, B, C connected in series. Each of them has a different Fermi energy and together form an energy quantum well which is tuneable upon application of a gate voltage (in the case studied here, we assume that the different gate capacitance on insulating and metallic YBCO lead to different gating effects in A, C and B). In order to allow for Andreev reflection at the A/B interface, the model [START_REF] Linder | Dirac fermions and conductance oscillations in s-and d-wave superconductor-graphene junctions[END_REF] assumes that the electronic density of states in region A presents a superconducting energy-gap. The homostructure conductance is obtained by matching the electronic wave functions in the A/B and C/B interfaces. Given the transmission T = 1 -|r| 2 and reflection R = |r A | 2 coefficient, the full expressions for r and r A are cumbersome and can be found in the appendix of Ref. [3]. The conductance across the A/B/C structure is given by:

G 1 (V 1 ) = π 2 -π 2 dθ 1 -|r| 2 cos (θ) + |r A | 2 cos (θ A ) . (SE1) With θ A = arcsin [(eV 1 + E f ) / (eV 1 -E f ) sin (θ)
] and θ the electron angle of incidence with respect to the interfaces.

Both r and r A and thus the conductance G 1 (V 1 ) depend on the following parameters: the reduced channel length Λ ≡ L/λ F,A (with λ F,A the Fermi wavelength in A), the reduced graphene Fermi energy in the different regions ( F,i ≡ E F,i /∆, with i=A, B, C) and ∆ the amplitude of the superconducting energy-gap, and the angle between the d-wave nodes of the superconducting order parameter and the interface α. The model makes a full wave matching calculation such that effects coming from different angles of incidence after several reflections are taken into account.

As discussed in the manuscript, this first block of the model explains the conductance oscillations as a function of gate and bias voltage. It is important to notice that the above result applies to electrons that propagate ballistically within the B-region.

The second block of the model takes into account the finite transmission of the YBCO/Au/graphene interface that is model as a d-wave/metal interface as studied in Refs. [START_REF] Wei | Directional Tunneling and Andreev Reflection on YBa2Cu3O 7-δ Single Crystals: Predominance of d -Wave Pairing Symmetry Verified with the Generalized Blonder, Tinkham, and Klapwijk Theory[END_REF] and [4]. Specifically, the conductance is given by:

G 2 (V 2 ) = π 2 -π 2 dθ 16 1 + |Γ + | 2 cos(θ) 4 + 4Z 2 1 -|Γ + Γ -| 2 cos(θ) 2 |4 cos(θ) 2 + Z 2 (1 -Γ + Γ -) | 2 , ( SE2
)
where

Γ ± = (eV 2 /|∆ (θ ± ) |) -(eV 2 /|∆ (θ ± ) |) -1 and θ + = θ, θ -= π -θ . While ∆(θ) = ∆ cos [2 (θ -α)]
and α is the angle between the superconducting order parameter and the interface. The conductance G 2 (V 2 ) thus depends on the parameters α and on the BTK barrier strength Z. As discussed in the main text this block of the model leads to a proper description of the background conductance vs. the voltage bias.

Using I = G (V ) dV we obtain G 1 (I) and G 2 (I) and from these, we calculate the conductance when the two building blocks are connected in-series:

G (I) = 1 G 1 (I) -1 + G 2 (I) -1 (SE3)
Finally, using V = 1 G(I) dI, we obtain G(V ). best agreement is found for the combination F,A = 20, Λ = 10.

Finally, F,B and F,C are adjusted in order to modulate the oscillations amplitude. As shown earlier for this model [3], the two latter parameters have no impact on the oscillations period, only on their amplitude and phase. This is because the Fermi velocity mismatch modifies the transmission/reflection across the interface. This is illustrated

IV. QUANTITATIVE EVALUATION OF THE SIMULATION PARAMETERS

The simulations are based on dimensionless parameters whose values for each device are given in a table below. It is very important that their consistency can be quantitatively evaluated using fundamental relations and the actual device dimensions. As explained, above the two parameters that determine the oscillations period are Λ are F,A .

Specifically, what defines the oscillations period is the ratio Λ/ F,A .

If we take the definition of the parameters Λ ≡ L/λ F,A and F,A ≡ E F,A /∆, and we consider the dispersion relationship in graphene E F = v F k F = hv F /λ F (with k F and v F = 10 6 m • s -1 the Fermi vector and velocity, and h the Planck constant), it is straightforward to find the following expression for the superconducting gap ∆:

∆ = Λ F,A hv F L -1 (SE4)
Equation (SE4) allows evaluating the simulations parameters quantitatively. If we input in Eq. (SE4) the actual graphene cavity length in the studied devices L device (distance between YBCO electrodes defined lithographically)

and the used Λ F,A , we obtain a value for superconducting gap, which is similar for most of the devices varies between 15 meV and 20 meV. This is well within the values expected for the same YBCO films when moderately underdoped [5] and materials in proximity with YBCO [6].

That a strong proof of the consistency of the analysis, from the quantitative point of view.

V.

VI. SIMULATION VS. EXPERIMENT FOR THE SET OF DEVICES

We present in the following the comparison between the experimental curves measured at 3.2 K (in blue) and the numerical simulation (in orange) based on the model detailed above. In principle and accordingly to [7], we do not expect the Au layer to significantly disturb the d-wave pairing, because it is one order-of-magnitude thinner than the mean free path in this material.

In our experiments, the key evidence for the propagation of the d-wave correlations into graphene is the observation of spectral features that cannot be explained based on s-wave superconductivity. In particular, the sharp zero-bias peak in some others [Fig. 3 This is illustrated by the example shown in Fig. S12, which corresponds to device B4D. The best agreement between theory and experiment is obtained with a barrier strength Z = 2.5 and an effective angle between d-wave nodes and junction interfaces α = 0.39.

In Fig. S13, we show simulations based on s-wave superconductivity both at the YBCO/Au-graphene interface and across the proximitized graphene homojunction. One can see that this assumption does not allow reproducing the experimental behaviour, particularly the sharp zero-bias conductance peak, even if the barrier strength Z is diminished to zero. Notice also that if the barrier strength Z is finite, a "dip" emerges around zero bias where s-wave quasiparticle tunnelling is strongly supressed. This is just opposite to the behaviour observed in the experiments.

Finally, we show in Fig. S14 simulations which consider d-wave superconductivity in the YBCO/Au-graphene interface and s-wave superconducting correlations across the graphene homojunction. One can see that under this assumption one cannot reproduce the behaviours observed experimentally.

The next example corresponds to device A5U, results in Fig. S15. We assume s-wave correlations across every junction interface. Again, the zero-bias conductance peak cannot be reproduced by considering only s-wave correlations-wave propagating the junction interfaces.

If we consider d-wave correlation only for the YBCO/Au-graphene interface and s-wave correlations across the homojunctions, the simulations do not allow reproducing the zero bias sharp zero bias conductance peak either, regardless of the barrier strength Z, see Fig. S16.

In summary, the experimental findings are well explained when it is assumed that d-wave superconducting correla- tions propagate into the graphene homojunction, and exclude the s-wave pairing component.

Notice finally that the effective angle between the d-wave nodes and the junction interfaces changes randomly from device to device, which implies that it is not fixed by the device geometry. This suggests that the the contact between the graphene and the superconducting electrodes might be inhomoegeneous and occur at edges of the SC electrodes or at exposed a-axis facets in the YBCO surface [6]. Local probes such as Scaning Tunneling microscopy would be necessary to investigate this question.

FIG. 2 .

 2 FIG. 2. Micrograph of a graphene/YBCO superconducting device. The scale bar is 10 µm. (a) shows a zoom on the junction: Superconducting electrodes (SC) are connected by graphene bridges (G) on top of insulating YBCO (INS). The scale bar is 100 µm. (c) Micrograph of a device with the top gate. (d) Scheme of the graphene homojunction model. The doping depends on whether graphene lies on insulating or superconducting YBCO, which results in three graphene regions A B, C with different Fermi energy. This creates the cavity where interferences occur. The lower scheme shows the equivalent electrical circuit. The YBCO-Au-graphene interface, characterized by a barrier strength Z, is measured in series with the cavity.

FIG. 3 .

 3 FIG. 3. Conductance versus VBIAS for two different devices, B3U (a) and B4U (b) measured at 3.2 K. Vertical lines point out the series of oscillations, and the horizontal double-headed arrows indicate the periods V long (black) and V short (red). The shaded area indicates the width 2δ of the superconducting-gap related feature. (b) and (e) show corresponding full-model calculations of the conductance with the Z and indicated in the legend (c) and (f) shows simulations for a proximitized graphene cavity alone, without finite-Z junction in series.

[

  Fig. 3(d)] or an enhancement [Fig. 4(a)] within a typical bias range |V BIAS | < δ exp ∼ 20 -70 mV highlighted by the grey shade. Those features are reminiscent of the conductance decrease/increase observed for e|V BIAS | < ∆ in superconducting/normal-metal junctions depending on the interface transparency

FIG. 4 .

 4 FIG. 4. (a) Conductance versus VBIAS for device A5U, measured at T = 4 K. (b) Same device conductance (color scale) as function of bias voltage (horizontal axis) and gate voltage (vertical axis). (c) Full-model calculation of the conductance with Z and α indicated in the legend. (d) corresponding fullmodel calculation of the conductance as a function of bias and the Fermi energy variation in the cavity.

  (b), 3(e) and Fig. 4(c) for the experiments in Fig. 3(a), 3(d) and Fig. 4(a), respectively. The simulations closely reproduce the main experimental features [further examples and the simulations' parameters shown in the Supplemental Material [36] in Secs. 5 and 6].

  FIG. S2. Differential conductance experimental (blue) and theoretical (orange) as a function of the bias for α = 0.23 and various values of the Z parameter: Z = 0 (top left), Z = 0.5 (top right), Z = 1 (bottom left) and Z = 2 (bottom right).

  FIG. S4. Differential conductance experimental (blue) and theoretical (orange) as a function of the bias for F,A = 20 and various values of the Λ parameter: Λ = 8 (top left), Λ = 10 (top right), Λ = 13 (bottom left) and Λ = 16 (bottom right).

FIG. S9 .

 S9 FIG. S9. Differential conductance experimental (blue) and theoretical (orange) as a function of the bias for different sets of fitting parameters: F,A = 12, Z = 0.25, α = 0.39 (top left), F,A = 16, Z = 0.25, α = 0.39 (top right), F,A = 14, Z = 0, α = 0.39 (central left), F,A = 14, Z = 0.5, α = 0.39 (central right), F,A = 14, Z=0.25, α = 0.23 (bottom left), F,A = 14, Z = 0.25, α = 0.53 (bottom right).

  II. Comparison between the experimental curves measured at 3.2 K (in blue) and the numerical simulation (in orange) based on the model detailed above for all the devices analyzed. VII. INDUCED SUPERCONDUCTING CORRELATIONS: d-WAVE VERSUS S-WAVE

  FIG. S12. Device B4D: Differential conductance experimental (blue) and theoretical (orange) as a function of the bias for Z = 2.5, α = 0.39.

  FIG. S13. Differential conductance experimental (blue) and theoretical (0.39) as a function of the bias for different sets of fitting parameters for the device B4D assuming a s-wave proximization: Z = 1 (top left), Z = 0.5 (top right), Z = 1 (bottom left), Z = 2 (bottom right).

  FIG. S15. Differential conductance experimental (blue) and theoretical (0.39) as a function of the bias for different sets of fitting parameters for the device A5U assuming a s-wave proximization: Z = 1 (top left), Z = 0.5 (top right), Z = 1 (bottom left), Z = 2 (bottom right).

TABLE OF

 OF

						SIMULATION PARAMETERS
	Device L device			Input parameters		∆
									from Eq. SE4
		nm	Z	α	Λ	F,A	F,B	F,C	meV
	A3D	100	2	0	10	20	400 17.5	20.7
	A3U	100	2	0	38	40	42.5 15.5	39.3
	A5U	150	0.25	0.39	10	14	150	9.5	19.7
	B3U	300	1	0.008	25	20	150 11.5	17.2
	B4U	300	2	0	21	20	400	15	14.4
	B4D	300	2.5	0.39	13	20	400	10	8.9
	E3D	800	2.5	0.15	80	25	80	20	16.5

TABLE I .

 I List of the fitting parameters for the devices analyzed.

Supplemental Material:

Long-Range Propagation and Interference of d-wave Superconducting Pairs in Graphene D. Perconte, 1, 2 K. Seurre, 1 V. Humbert, 1 C. Ulysse, 3 A. Sander, 1 J. Trastoy, 1 V. Zatko, 1 F. Godel, 1 P. R. Kidambi, 4, S. Hofmann, 5 X. P. Zhang, 6, 7 D. Bercioux, 7, 8 F. S. Bergeret, 6, 7 B. Dlubak, 1 P. Seneor, 1 and Javier E. Villegas We first find Z and α that allow mimicking the conductance background. These parameters respectively characterize the YBCO/Au-graphene interface and the d-wave nature of the superconducting correlations.

Let us first assume that the normalized cavity length Λ → 0 (so that there are no visible oscillations), fix α = 0 and vary Z, the results are presented in the various panels of Fig. S1. We see that one needs relatively large Z = 2 to better reproduce the V-shaped "dip" around zero bias. Now we explore the effect of varying α in Figs. S2 andS3. With these examples, one quickly realizes that it is not possible to use α = 0 in the present case, regardless of the choice of Z, because α = 0 yields to a sharp zero/bias conductance peak that is absent for this particular device. Considering the above, we would fix Z = 2 and α = 0 for the A3D device. Now we look for Λ and F,A that yield oscillations of the same reduced period ∆V /δ as in the experiments. We show in Figs. S4 and S5 some examples that illustrate how these parameters influence the oscillations period. First, we show examples with fixed F,A and varying Λ (the red vertical lines mark minima in the experimental curve, allowing for a reference of experimental high-bias period): We see that increasing the cavity length Λ diminishes the oscillations period, as one would naturally expect. As shown in the Figs. S4 and S5, if one fixes Λ and varies F,A the period is also modified, which is also expected. Indeed, from the definition of the parameters and the graphene's linear dispersion relation [3], the oscillations period is the same for any constant Λ/ F,A . Note however that the choice of F,A is not entirely free because this parameter modifies may introduce an asymmetry in the curve. Globally the