Matthieu Bellec
email: matthieu.bellec@inphyni.cnrs.fr

Charles Poli

Ulrich Kuhl

Fabrice Mortessagne
email: fabrice.mortessagne@inphyni.cnrs.fr

Henning Schomerus
email: h.schomerus@lancaster.ac.uk

Observation of supersymmetric pseudo-Landau levels in strained microwave graphene

Using an array of coupled microwave resonators arranged in a deformed honeycomb lattice, we experimentally observe the formation of pseudo-Landau levels in the whole crossover from vanishing to large pseudomagnetic field strength. This is achieved by utilizing an adaptable set-up in a geometry that is compatible with the pseudo-Landau levels at all field strengths. The adopted approach enables us to observe fully formed flat-band pseudo-Landau levels spectrally as sharp peaks in the photonic density of states, and image the associated wavefunctions spatially, where we provide clear evidence for a characteristic nodal structure reflecting the previously elusive supersymmetry in the underlying low-energy theory. In particular, we resolve the full sublattice polarization of the anomalous 0th pseudo-Landau level, which reveals a deep connection to zigzag edge states in the unstrained case.

I. INTRODUCTION

Topological states enjoy intense attention as they equip quantum systems with desirable robust properties. Much of the early focus rested on their unique spectral positions as isolated or dispersive states in a band gap, as well as their spatial localization at edges and interfaces [START_REF] Hasan | Colloquium: Topological insulators[END_REF], or more recently also corners [START_REF] Benalcazar | Quantized electric multipole insulators[END_REF]. More fundamental characterizations, on the other hand, often invoke a third, somewhat deeper feature of topological states, which is connected to the anomalous expectation values of the underlying symmetry operators [START_REF] Nielsen | The Adler-Bell-Jackiw anomaly and Weyl fermions in a crystal[END_REF][START_REF] Callan | Anomalies and fermion zero modes on strings and domain walls[END_REF]. In momentum space, this feature underpins, e.g., the unidirectional chiral currents around the edges of topological insulators, whilst in real space it manifests itself, e.g., in the sublattice polarization of defect states in bipartite lattice systems [START_REF] Poli | Selective enhancement of topologically induced interface states in a dielectric resonator chain[END_REF], as has been exploited in recent topological lasers or non-linear limiters based on photonic Su-Schrieffer-Heeger structures [START_REF] St-Jean | Lasing in topological edge states of a onedimensional lattice[END_REF][START_REF] Zhao | Topological hybrid silicon microlasers[END_REF][START_REF] Parto | Edge-mode lasing in 1d topological active arrays[END_REF][START_REF] Reisner | Self-shielded topological receiver protectors[END_REF].

An additional attractive aspect of these anomalous features are that they tie topological effects together that are often seen as separate, due to the varied nature of the specific encountered spectral and spatial features that first come into focus. A prime example are flat bands, which have been observed in recent experiments focussing on Lieb lattices [START_REF] Guzmán-Silva | Experimental observation of bulk and edge transport in photonic Lieb lattices[END_REF][START_REF] Mukherjee | Observation of a localized flat-band state in a photonic Lieb lattice[END_REF][START_REF] Diebel | Conical diffraction and composite Lieb bosons in photonic lattices[END_REF][START_REF] Poli | Partial chiral symmetry-breaking as a route to spectrally isolated topological defect states in two-dimensional artificial materials[END_REF] and one-dimensional counterparts [START_REF] Baboux | Bosonic condensation and disorder-induced localization in a flat band[END_REF], as well as suitably deformed graphene [START_REF] Levy | Strain-induced pseudo-magnetic fields greater than 300 tesla in graphene nanobubbles[END_REF] and analogous quantum [START_REF] Gomes | Designer Dirac fermions and topological phases in molecular graphene[END_REF] and classical systems [START_REF] Rechtsman | Strain-induced pseudomagnetic field and photonic Landau levels in dielectric structures[END_REF][START_REF] Wen | Acoustic Landau quantization and quantumhall-like edge states[END_REF]. In the latter case they constitute pseudo-Landau levels arising from a synthetic magnetic field [START_REF] Guinea | Energy gaps and a zero-field quantum Hall effect in graphene by strain engineering[END_REF]. In particular, signatures of photonic pseudo-Landau levels have been detected by probing the edges of a honeycomb array of optical waveguides [START_REF] Rechtsman | Strain-induced pseudomagnetic field and photonic Landau levels in dielectric structures[END_REF]. A second example are a class of helical edge states in reciprocal systems, as observed, e.g., in zigzag terminated graphene [START_REF] Fujita | Peculiar localized state at zigzag graphite edge[END_REF][START_REF] Salerno | Propagating edge states in strained honeycomb lattices[END_REF]. While these bulk and edge phenomena do not naturally fall into the scope of standard topological band structure theory [START_REF] Hasan | Colloquium: Topological insulators[END_REF], they are still intimately linked to wavefunctions with a characteristic sublattice polarization. This provides a promising perspective from which one can seek to develop very general unifying descriptions (see, e.g., Ref. [START_REF] Kunst | Lattice models with exactly solvable topological hinge and corner states[END_REF] for a recent approach utilizing this perspective).

In this work, we demonstrate experimentally for the case of photonic graphene-like systems that the anomalous edge and bulk phenomena tied to sublattice polarization are in fact directly linked. This is achieved by tracing the formation of pseudo-Landau levels all the way from vanishing to large pseudomagnetic field strengths. In particular, we report the direct observation of the spatially resolved sublattice polarization in the 0th pseudo-Landau level of strained photonic graphene, and trace it back to the unstrained case, where the system only possesses edge states. By observing a characteristic nodal structure for the higher-order levels, we can then establish a direct link to the supersymmetric Hamiltonian of the underlying low-energy theory [START_REF] Jackiw | Fractional charge and zero modes for planar systems in a magnetic field[END_REF]. Thereby, our observations connect a broad variety of topological phenomena to a unifying principle.

II. RESULTS

A. Microwave set-up and optimal strain geometry

Our experimental set-up is illustrated in Fig. 1. The unstrained system forms a honeycomb lattice with nearest-neighbor spacing a 0 = 13.9 mm, combining two triangular lattices of A and B sites, where each vertex denotes the position of a dielectric microwave resonator with bare frequency ω 0 = 6.653 GHz, while adjacent resonators are coupled at strength t 0 = 21.5 MHz (for details see the materials and methods section). This gives rise to a standard graphene-like photonic band structure [START_REF] Bellec | Topological transition of Dirac points in a microwave experiment[END_REF][START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF], with two Dirac cones at the K and K points in the Brillouin zone. Around these two so-called valleys, at relative momentum q, the low-energy dispersion ω(q) ∼ ω 0 ± v|q| resembles massless relativistic particles moving in two dimensions at velocity v = 3a 0 t 0 /2. 1) and the text for details). b) In the experiment, the lattice is composed of 196 identically designed cylindrical dielectric resonators that are coupled through the evanescent field of the fundamental TE mode. The structure is placed inside a microwave cavity made of two metallic plates (top plate only partially shown). A loop antenna, mounted on a scanning system (white arrows) crossing the top plate and connected to a vectorial network analyzer (VNA), is used to generate and collect the microwave signal both spectrally and spatially resolved, which allows to obtain the local density of states in the system.

In the deformed system, the indicated couplings t l depend on the distances to the three neighboring resonators, which we can utilize to create a pseudomagnetic field corresponding to that in strained graphene. Such a field arises when the resonators are displaced nonuniformly, where the positions are selected to give a triaxial spatial coupling profile [START_REF] Guinea | Energy gaps and a zero-field quantum Hall effect in graphene by strain engineering[END_REF][START_REF] Schomerus | Parity anomaly and Landau-level lasing in strained photonic honeycomb lattices[END_REF]

t l = t 0 [1 -(β/2a 2 0 )ρ l • r], (1)
where r refers to the positions of the links between the coupled resonators in the unstrained system and the bond vectors ρ l are pointing along these coupling directions. This coupling profile ensures a constant pseudomagnetic field of strength β throughout the whole system. Theoretically, the system is well described in a coupledmode theory with nearest-neighbour couplings t l as given above, so that the eigenfrequencies and mode profiles can be obtained from an effective Hamiltonian H. As any such bipartite system, it then displays a chiral symmetry relative to the central frequency,

Σ z (H -ω 0 )Σ z = -(H -ω 0 )
, where the Pauli-like matrix Σ z acts on the sublattice degree of freedom, hence keeps the amplitudes on A sites fixed but inverts those on the B sites. The balance of zero modes on the A and B sublattices is given by the signature of this operator [START_REF] Poli | Partial chiral symmetry-breaking as a route to spectrally isolated topological defect states in two-dimensional artificial materials[END_REF][START_REF] Sutherland | Localization of electronic wave functions due to local topology[END_REF],

#(A zero modes) -#(B zero modes) = tr Σ z = #(A sites) -#(B sites). (2)
These zero modes have frequency ω 0 , and as indicated are localized on a given sublattice. Furthermore, the chiral symmetry dictates that all non-zero modes occur in spectral pairs at symmetric positions ω 0 ± δω, and have equal intensity on both sublattices. Importantly, and in contrast to earlier experimental work, we select a triangular geometry and terminate the system with zigzag edges. This ensure a number of beneficial features [START_REF] Poli | Degeneracy doubling and sublattice polarization in strain-induced pseudo-Landau levels[END_REF][START_REF] Rachel | Straininduced Landau levels in arbitrary dimensions with an exact spectrum[END_REF]. Particularly relevant for us, the boundary conditions are then compatible with the bulk pseudo-Landau levels at all field strengths; furthermore, a consistent coupling profile can be maintained even at maximal field strength, whose description requires to go beyond the conventional low-energy analogy to magnetic fields with opposite signs in the two valleys. Setting the pseudomagnetic field strength to a fixed value, we then see that the values of the couplings dictated by the profile (1) drop to zero exactly at the terminating edges of a zigzag terminated triangle. The corresponding maximal field strength is β M = 4/L, where the size parameter L counts the terminating A sites along each edge.

In the experiments we realize these conditions in a system with 196 resonators, corresponding to a triangle with L = 14 resonators along each terminating edge (see Fig. 1). Of these, 105 resonators are on the A sublattice while 91 are on the B sublattice. This allows us to realize field strength up to β = 0.2, where the extremal couplings still exceed the homogeneous resonator linewidth γ = 1.7 MHz, and sufficiently close to β M to clearly demonstrate the detailed features of well-formed pseudo-Landau levels. To obtain the analogous orbital effects for an electron in graphene, a magnetic field of 42000 T would have to be applied.

B. Formation of pseudo-Landau levels

Figure 2 shows our main experimental results. The panels on the left show the density of states for the resonator geometry on the right, where each row corresponds to a different value of the pseudomagnetic field strength β. The color density plot overlayed with the resonator lattice depicts the local density of states integrated over the central peak, situated at the bare resonator frequency ω 0 .

In the pristine system (β = 0), this peak arises from the zigzag edge states, which are localized on the terminating resonators. Note that these resonators all occupy the same sublattice of A sites. On this sublattice, the zigzag states decay into the bulk of the system, whilst they maintain a vanishing density on the B sublattice. The simple rule (2) can be exploited to count the number of these zigzag states: As a consequence of the chiral symmetry of the system, this number is expected as the difference 105 -91 = 14 = L of A and B sites in the system. Away from the central peak, the density of states displays a broad continuum, with fluctuations arising from the finite-size quantization of bulk graphene-like states, including the states near the Dirac cones. The formation of the 0th pseudo-Landau level proceeds continuously by transforming zigzag edge states into bulk states, whilst retaining the degeneracy (spectral weight) and sublattice polarization (remaining confined on the A sublattice).

As the pseudomagnetic field strength is increased, we observe that the spectral weight from this continuum reorganises into a sequence of well-defined peaks, of similar width and weight as the central peak, whose weight and position remains essentially unchanged. The spatial profile of the associated zero-modes, however, changes significantly, in that they move into the bulk, where they form the desired 0th pseudo-Landau level.

The observed spectral positions of the emerging higherorder pseudo-Landau levels conform well with the characteristic square-root dependence on relativistic Landau levels [START_REF] Jackiw | Fractional charge and zero modes for planar systems in a magnetic field[END_REF][START_REF] Castro Neto | The electronic properties of graphene[END_REF], revisited later in the text and depicted by the gray dashed lines in Fig. 2. The same applies to the observed spectral weights. As the transformation from the zigzag edge states to the 0th pseudo-Landau level is continuous, this bulk level retains the same spectral degeneracy, hence here consists of L = 14 modes. The coupled-mode theory predicts that the nth pseudo-Landau level encompasses 14 -|n| states, which explains the gradual drop of the observed spectral weights of the peaks moving outwards from the central peak. Note that this implies an important difference to Landau levels arising from a magnetic field, for which the degeneracy is dictated by the sample area ∝ L 2 , but not the linear size L, as observed here and underpinned by general theory [START_REF] Poli | Degeneracy doubling and sublattice polarization in strain-induced pseudo-Landau levels[END_REF][START_REF] Rachel | Straininduced Landau levels in arbitrary dimensions with an exact spectrum[END_REF].

C. Supersymmetric nodal profiles

As anticipated, a notable feature in the formation of the 0th pseudo-Landau level verified in the experiment is the observation that across the whole transition, the associated modes remain localized on the A sublattice. In contrast, the modes in the emerging higher-order levels, whilst also localized in the bulk, are anticipated to display an equal weight on both sublattices. We analyze this distinction in detail in Fig. 3. The top rows display the experimental local density of states in the spectral range of the pseudo-Landau levels for the indices n = 0, 1, 2, 3, 4, all taken at the large pseudomagnetic field strength β = 0.2. As above, the level with n = 0, shown in the left panels, is located on the A sublattice. The higher-order levels shown in the other panels indeed display an approximately equal weight on both sublattices. Furthermore, they have an intensity profile that increasingly seeps into the corner areas of the triangle, which remains in line with the predictions of coupledmode theory where the states form a complete basis.

The other two rows in Fig. 3 show the local density of states along the edges of the system separately for the A and B sublattice, where we averaged over the three edges and indicate the range of observed values by the shaded areas. For the higher-order levels, we furthermore include the levels with index -n into the average. Along each edge, we observe standing-wave patterns

ψ (A, edge) r ∝ sin[(|n| + 1)πr/(L + 1)] (3)
for resonator index r = 1, 2, 3, . . . , L on the A sublattice, and a corresponding pattern

ψ (B, edge) r ∝ sin(|n|πr/L) (4)
for resonator index r = 1, 2, 3, . . . , L -1 on the B sublattice.

Note that the mode index (|n| + 1) vs |n| in these patterns for the two sublattices is offset by one. This observation finds its counterpart in the low-energy theory of the pseudo-Landau levels relative to the central frequency ω 0 [START_REF] Jackiw | Fractional charge and zero modes for planar systems in a magnetic field[END_REF][START_REF] Schomerus | Parity anomaly and Landau-level lasing in strained photonic honeycomb lattices[END_REF][START_REF] Castro Neto | The electronic properties of graphene[END_REF], which are described by an effective Hamiltonian

H 0 = 0 π † π 0 , ( 5
) LL 0 A sites LL 1 LL 2 LL 3 LL 4 B sites 0 1 A sites 1 r max 0 1 B sites 1 r max 1 r max Edge site index r 1 r max 1 r max

Intensity (arb. units)

FIG. 3. Supersymmetric node structure. Experimentally determined spatially resolved mode intensities associated with the pseudo-Landau levels of index n = 0 to 4 (from left to right) for a large pseudomagnetic field strength β = 0.2. Top rows: the area of the circles correspond to the intensity on the A sites (above) and B sites (below). The 0th level is almost entirely localized on the A sublattice, while all other levels have almost-equal overall weight on both sublattices. Bottom rows: intensity on A sites (above) and B sites (below), averaged along the three extreme edges at resonator position r. Note that rmax = 14 (resp. 13) for A (resp. B) sites. For the higher-order levels, the average is performed combining the indices ±n. The shaded area corresponds to the standard deviation. The two y-axis ticks indicated intensity values of 0 and 1 in arbitrary, but uniformly applied units. The pseudo-Landau levels display a clear nodal structure offset by 1 mode index, as further discussed in the text.

where the Landau-level creation operators fulfill [π, π † ] = 2v 2 β/a 2 0 (we assume β > 0; for opposite deformation the operators π and π † interchange their roles of creation and annihilation operators). This corresponds to the Hamiltonian of a relativistic electron in a magnetic field, which is formally identical to a supersymmetric harmonic oscillator [START_REF] Jackiw | Fractional charge and zero modes for planar systems in a magnetic field[END_REF][START_REF] Thaller | The Dirac Equation[END_REF]. The link to the standing-wave patterns (3) and ( 4) along the edges arises from the fact that H 2 0 = diag(π † π, ππ † ) factorizes into two equidistant level sequences E 2 n = 2v 2 nβ/a 2 0 on the two sublattices, where n = 0, 1, 2, 3 . . . on the A sublattice while on the B sublattice n = 1, 2, 3 . . . so that the sequence is offset by one. While this low-energy analogy does not capture the difference in degeneracies and the importance of boundary conditions described above, it is notable that the offset between the level sequences persists in the experimental observations, and is furthermore reflected in the mentioned nodal wave-patterns.

III. DISCUSSION

We achieved the direct observation of the formation of pseudo-Landau levels in deformed honeycomb systems, both spectrally as well as in terms of their key spatial features. In particular, adopting a flexible dielectricresonator array design with a purposefully selected geometry allowed us to follow the formation of these levels in the transition from vanishing to large pseudomagnetic field strengths. In this way, we could observe how the 0th pseudo-Landau level originates from the transformation of zigzag edge states into bulk states, whilst maintaining its characteristic anomalous polarization on only one of the two sublattices in the system. Extending these considerations to the higher-order levels allowed us to reveal a characteristic nodal structure of the pseudo-Landau level sequence that reflects the supersymmetric structure of the underlying low-energy description. These features underline the general usefulness to account for anomalous expectation values to provide a more general perspective on topological states.

Resolving the reported features in an experiment poses a significant challenge. In electronic systems such as graphene, spectral imaging techniques do not provide the required atomistic resolution, so that the information can only be extracted indirectly, for example from the Fourier transformation of form factors [START_REF] Daimon | Structure factor in photoemission from valence band[END_REF][START_REF] Bostwick | Quasiparticle dynamics in graphene[END_REF]. A better resolution is offered by photonic systems [START_REF] Rechtsman | Strain-induced pseudomagnetic field and photonic Landau levels in dielectric structures[END_REF], which however so far could not access the characteristic spatial features of the pseudo-Landau levels in the bulk of the system, and furthermore considered a geometry that does not enable to observed fully formed flat-band pseudo-Landau levels [START_REF] Poli | Degeneracy doubling and sublattice polarization in strain-induced pseudo-Landau levels[END_REF]. Drawing on an acoustic analogue [START_REF] Abbaszadeh | Sonic Landau levels and synthetic gauge fields in mechanical metamaterials[END_REF], a recent experiment [START_REF] Wen | Acoustic Landau quantization and quantumhall-like edge states[END_REF] managed to excite a compacton-like state in the 0th Landau level, demonstrating its characteristic sublattice polarization in the bulk [START_REF] Schomerus | Parity anomaly and Landau-level lasing in strained photonic honeycomb lattices[END_REF]. Here, we exploited an adaptable dielectric microwave-resonator array geometry to provide a complete characterization of the system from the unstrained to the fully strained case. This allows us to reveal how zigzag edge states transform into the bulk states of the anomalous 0th pseudo-Landau level, where they retain their characteristic sublattice polarization. Furthermore, this perspective dictates a natural geometry in which maximal pseudomagnetic fields can be attained, and for which the pseudo-Landau levels remain compatible with the boundary conditions at all field strengths, in contrast to the previous experiments. This is required to obtain pseudo-Landau levels that are flat, which we evidence spectrally by the observation of sharp peaks in the photonic density of states, and enables to reveal the supersymmetric signatures in the nodal structure.

Our results extend to a wide variety of flat-band systems, such as the rich physics arising from higher-order resonator modes in deformed honeycomb lattices, as reported for exciton polaritons in Ref. [START_REF] Jamadi | Direct observation of photonic Landau levels and helical edge states in strained honeycomb lattices[END_REF]. Practically, the observed spatial features should help to pave the way to applications such as flat-band lasers [START_REF] Poli | Partial chiral symmetry-breaking as a route to spectrally isolated topological defect states in two-dimensional artificial materials[END_REF][START_REF] Schomerus | Parity anomaly and Landau-level lasing in strained photonic honeycomb lattices[END_REF][START_REF] Smirnova | Parity anomaly laser[END_REF], as well as sublattice-dependent sensors where one could exploit that chiral-symmetry breaking perturbations equip the 0th Landau level with a finite weight on the opposite sublattice.

IV. MATERIALS AND METHODS

A. Experiment

The experimental setup is designed to realize a microwave system that is well approximated by a nearestneighbour tight-binding description [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF]. The sites of the lattice are occupied by dielectric microwave resonators with a cylindrical shape and made of ZrSnTiO ceramics (Temex-Ceramics, E2000 series: 5 mm height, 8 mm diameter and a refractive index n ≈ 6) sandwiched between two metallic plates at a distance h = 16 mm. Each resonator supports a fundamental TE mode of bare fre-quency ω 0 = 6.653 GHz, which corresponds to the onsite energy of atoms in a tight-binding model. As the resonance frequency is below the cut-off frequency of the first TE mode defined by the two plates the adjacent resonators are coupled through evanescent wave components leading to an approximately exponential decay of the coupling strength t with the distance between the resonators [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF]. The system is excited via a loop antenna fixed in the movable top plate thus allowing to scan spatially the magnetic field B z , which is the only magnetic field component for this mode [START_REF] Reisner | Microwave limiters implemented by coupled dielectric resonators based on a topological defect mode and CT-symmetry breaking[END_REF]. From the reflection measurements performed by a vector network analyzer (ZVA 24 from Rohde & Schwarz) the local density of states can be extracted (for details see Ref. [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF]) and finally, by integrating over space, the density of states. In all these experiments, we face an intrinsic on-site disorder of ∼ 0.15% in the values of ω 0 .

B. Modelling

For the design of the system, we modelled finite strained photonic honeycomb lattices with a range of system sizes and boundary geometries within coupled-mode tight-binding theory, in which we accounted for the fundamental TE mode and incorporated the experimental distance dependence of the coupling strengths following Ref. [START_REF] Bellec | Tight-binding couplings in microwave artificial graphene[END_REF]. Exact diagonalization gives us access to the resonant modes and their spatial intensity distribution. The modelling confirmed that the triaxial coupling profile (1) can be attained by suitable positioning of the resonators, resulting in an excellent match with the continuum theory predictions for the lowest Landau levels in systems as used in the experiments. The modelling further confirmed that only for zigzag boundaries one can attain the full sublattice polarization of the zeroth Landau level, and only for a triangular shape the Landau levels become maximally degenerate.

VI. AUTHOR CONTRIBUTIONS

MB and FM carried out the experiments and analyzed the data with input from UK. CP modelled the system with input from HS, who provided the theory and initiated the project. All authors contributed substantially to this work, including the interpretation of the results and the preparation of the manuscript.

FIG. 1 .

 1 FIG. 1. Experimental set-up. a) Sketch of the unstrained (gray) and strained (black) honeycomb lattice geometry, forming a zigzag terminated triangle of size L = 14. The green arrows indicates the directions ρ l , l = 1, 2, 3, of the effective triaxial strain. Inset: the underlying lattice is composed of two sublattices A (red) and B (blue). The local coupling strengths are denoted t l . They depend on the location of the bonds in the lattice [see Eq. (1) and the text for details). b) In the experiment, the lattice is composed of 196 identically designed cylindrical dielectric resonators that are coupled through the evanescent field of the fundamental TE mode. The structure is placed inside a microwave cavity made of two metallic plates (top plate only partially shown).A loop antenna, mounted on a scanning system (white arrows) crossing the top plate and connected to a vectorial network analyzer (VNA), is used to generate and collect the microwave signal both spectrally and spatially resolved, which allows to obtain the local density of states in the system.

2 FIG. 2 .

 22 FIG.2. Landau-level formation. Experimentally determined density of states (left panels) and spatially resolved mode intensity associated with the central peak centered at ω0 = 6.653 GHz (right panels, obtained by integrating the local density of states over the peak). The area of the circles corresponds to the intensity on the A sites (red) and B sites (blue). From top to bottom, the pseudomagnetic field strength β varies from 0 to 0.2. The gray dashed lines in the total density of states for β = 0.2 depict the expected pseudo-Landau levels frequencies from coupled-mode theory. The formation of the 0th pseudo-Landau level proceeds continuously by transforming zigzag edge states into bulk states, whilst retaining the degeneracy (spectral weight) and sublattice polarization (remaining confined on the A sublattice).

  and FM are grateful to G. Salerno and I. Carusotto for stimulating discussions during their visit to Nice. The authors acknowledge funding by EPSRC via Grants No. EP/J019585/1, EP/P010180/1, and Programme Grant No. EP/N031776/1, and CNRS through a visiting fellowship of HS to Nice.


Matthieu Bellec

Claire Michel

Disordering a superfluid of light

Keywords:

come

However, low-dimensionality is actually not always a restriction. 1D or 2D photonic devices can lead to rich and unexpected phenomena that may offer applications in both classical and quantum optics. For instance, topological photonics, a rapidly growing research field, is steeped in the 2D quantum Hall solid-state physics effect. Even the simpler topological photonic insulator 1D counterpart recently led to intense investigations. Feasible applications in robust optical transport, topological lasing and quantum communication have been proposed and realized [1].

Likewise, the Berezinskii-Kosterlitz-Thouless (BKT) transition [2] is another example showing that low dimensional condensed matter physics yields rich phenomenon. This nonconventional phase transition from disordered to quasi-ordered states in 2D systems, mediated by topological defects underlying dynamics, has been observed in quantum systems. Nonetheless, this transition mainly involves classical nonlinear wave dynamics. Now writing in Nature Photonics, Situ and Fleischer reported the first observation of the BKT transition in optics, using a continuous laser in a crystal exhibiting a large optical nonlinear response, creating a 2D 'fluid of light' evolving along the propagation axis [3]. In their work, the associated topological defects are optical vortices, defined as phase singularities in the transverse optical wave field.

The science of structured light with singularities has already led to strong basic and applied scientific interests. Thus, besides the fundamental importance of this observation, the disclosed results, involving complex optical vortex dynamics, naturally appeal to potential applications in singularity-based photonics 2D devices [4].

An important property of low-dimensional physical systems (1D or 2D) is the absence of long range order, like e.g. crystalline, ferromagnetic or superfluid states, due to inherent thermal fluctuations at finite temperature. When short range interactions are additionally considered, Berezinskii and Kosterlitz/Thouless independently predicted in the 70s, that a phase transition between ordered and disordered states may occur in 2D systems [2]; it's usually said to be an unconventional transition, being continuous, mediated by topological defects, and without symmetry breaking. While 2D crystals or magnets are rather complicated to implement and study experimentally, two-dimensional fluids of interacting particles, possibly obeying Bose-Einstein statistics, are much more accessible. Examples of such systems range from liquid Helium and superconductors to cold atoms and, more recently, quantum fluids of light. The order is then characterized by a high degree of phase coherence at long range, which is also called the superfluid state.

The role of the vortices in the BKT transition can be intuitively understood as follows. Consider the lowest energy state in a quantum fluid which has constant amplitude with flat phase. Adding a single vortex, defined by a wavefront spinning around a vanishing amplitude point, affects the field in a global way. At high temperatures, random motion of such vortices, which proliferate as the temperature increases, destroys large-scale phase coherence. If now two vortices spinning in opposite direction are paired together, which may occur as the temperature is decreased, the phase pattern is altered only locally and their circulation actually cancels, thus allowing long range order to appear. More precisely, the spatial phase correlations decrease exponentially in the disordered (normal) state and algebraically (power law) in the quasi-ordered (superfluid) state. At the BKT transition, which occurs at a critical temperature 𝑇 ! , a universal jump in the superfluid density is predicted. The first experimental measurements of this jump took place soon after the BKT predictions in helium films [5]. In 2006, the transition was observed in ultracold atoms and first attempts to probe vortices dynamics were reported [6]. More recently, quantum fluids of light have been proposed as a rich alternative playground to quantitatively study advanced quantum hydrodynamics features [7], including superfluidity [7,8] and the BKT transition. Algebraic correlation functions were reported with exciton-polaritons in 2D semiconductor microcavities [START_REF] Roumpos | Proc. Natl. Acad. Sci[END_REF]. Inspired by a pioneering theoretical study [START_REF] Small | [END_REF], Situ and Fleischer characterize the optical BKT transition in a nonlinear optical setup allowing investigation of both the phase coherence features of the fluid of light, and the complex vortices underlying dynamics of the transition [3].

The experimental setup consists of a continuous laser beam propagating into an optical crystal with a large nonlinear photorefractive response (see Figure 1). As in other nonlinear media such as atomic vapors [11] or thermal liquids [12], under the paraxial approximation, the transverse evolution of the optical beam is analogous to a weakly interacting Bose-Einstein condensate (BEC). The propagation axis plays the role of time, while nonlinearities and refractive index modulation act as effective interactions and external potential. Such nonlinear optical dynamics is usually reformulated in terms of quantum hydrodynamics where the laser intensity and wavefront gradient correspond to the fluid density and velocity respectively. To deepen the analogy, a uniform wavefront scenario is typically analogous to a zero temperature BEC. Adding a controlled amount of random-phase components in the beam, achieved experimentally by using a spatial light modulator, allows increasing the condensate's effective temperature. This imposes density and phase fluctuations of the fluid of light, eventually leading to the generation of vortices, and enabling triggering of the BKT transition. In principle, noise alone should provoke the generation of the vortices, but in practice it is not sufficient. Here, a lattice is photo-induced in the photorefractive crystal, to help trigger vortices, as discussed below.

For each parameter set, varying the effective temperature and the strength and sign of the interaction, the output complex wavefunction (i.e. amplitude and phase of the optical field) is reconstructed via a standard interferometric technique. Two main quantities are then extracted to characterize the BKT transition: i) the number of new vortices that have been generated along the beam path along with their pairing status (vortices are considered as paired if their separation is smaller or equal to three vortex cores); ii) the radial correlation function. Its cumulative integral can actually be expressed, for all temperatures, 𝑇, simply as Ω "#(%) where Ω is the integrating area and 𝜅 is an exponent to be fitted. In theory, by looking at 𝜅 as a function of 𝑇, the phase transition can be identified as a jump from 𝜅 = 0.875 at just below a critical temperature 𝑇 ! , to 𝜅 = 0.5 just above 𝑇 ! , corresponding to the algebraic (exponential) decay of the spatial correlations discussed above.

The presented experimental results show that the correlation functions obey the exact scaling behavior predicted by the BKT theory. The values of 𝜅 are directly linked to the (universal) jump in the superfluid density. It's worth mentioning that the expected step function at 𝑇 ! was not observed; instead a smooth crossover was seen, an effect attributable to the finite size of the experimental system, which is an effect also observed in ultracold atoms experiments [6].

Above the critical temperature, the number of single free vortices that proliferate drastically increases. Counting them as function of 𝑇 allows to quantitatively obtain the value of 𝑇 ! . Interestingly, the researchers observe that it actually depends on two parameters. Typical values of the normalized critical temperature range from 0.6 and 0.95 in accordance with reported results. First, 𝑇 ! is slightly shifted forward when the nonlinear strength is larger, whatever its sign. This is expected as the role of the density fluctuations is smaller in this case. Second, they measure a linear dependence of 𝑇 ! with the lattice spacing. A thermodynamics approach tells that smaller lattice spacings give more sites available for occupation (per a fixed area), increasing the entropy and indeed lowering the critical temperature needed to overcome the vortex energy. For a homogeneous system (i.e. when lattice spacing approaches zero) the transition should occur at finite 𝑇 ! . What the researchers observe is, however, a wave condensation phenomenon with no evidence of vortex proliferation [13]. This raises the important issue of the role of the lattice in these experiments and the interplay between continuous and discrete dynamics. Situ and Fleischer found that the experimental parameter ranges are limited and the vortices can't be trigged by noise only (temperature fluctuations) in homogeneous fluids of light, thus requiring an additional lattice.

The results reported by Situ and Fleischer are of fundamental importance and constitute a starting point for exploring the full phase space of non-equilibrium and vortex-coherence dynamics in fluids of light. An interesting aspect to follow-up on may be the role of the lattice arrangement (with e.g. topological or disorder features) on the transition. Fluids of lightthanks to the various and complementary existing photonic platforms -offer a rich playground to quantitatively study advanced quantum hydrodynamics features. Although restricted to 2D geometries, this work also shows that photonic systems generally exhibit complex dynamics with new light transport properties that may potentially lead to important applications in photonics.

Matthieu Bellec and Claire Michel Institut de Physique de Nice, Université Côte d'Azur & CNRS, Nice, France

Emails : matthieu.bellec@inphyni.cnrs.fr, claire.michel@inphyni.cnrs.fr


Lei Wang

Dr Yi-Ming Chen
email: chenym@ysu.edu.cn

Gang Cheng

Thierry Barriere

Thierry Barrière

P R China

to polymeric visco-elastic Euler-Bernoulli beam under quasi-static loads

Keywords: Fractional calculus, Fractional partial differential equation, Fractional rheological models, Visco-elastic properties, Euler-Bernoulli beam, Shifted Chebyshev polynomials Keywords:Fractional calculus, Fractional partial differential equation, Fractional rheological models, Visco-elastic properties, Euler-Bernoulli beam, Shifted Chebyshev polynomials

Numerical analysis of fractional partial differential equations applied to polymeric visco-elastic Euler-Bernoulli beam under quasistatic loads

Introduction

The polymeric material has gradually received the researchers' attentions due to its elastic and viscous behavior and high-efficiency vibration damping effect. As a vibration attenuating material, the polymer absorbs the vibration energy during mechanical vibration and converts it into heat, electric, magnetic energy to be consumed, thereby effectively reducing the structural vibration and noise. Nowadays, the polymeric material has been widely used in the fields of aerospace, construction and mechanical engineering.

*Manuscript Click here to view linked References

Poly (ether ether ketone) (PEEK), a special engineering polymeric material, is a semicrystalline thermoplastic polymer with excellent mechanical properties including high-temperature resistance, high-rigidity, high fracture toughness, dimensional stability and stable chemical properties. It plays an increasingly important role in electronics, aerospace, automotive and manufacturing field. The visco-elastic effects in PEEK and its composites have been specifically studied by D'Amore et al. [START_REF] Pompo | Viscoelastic effects in poly (ether ether ketone) (PEEK) and PEEK-based composites[END_REF] .

High-density polyethylene (HDPE) is a highly crystalline and non-polar thermoplastic polymer with high-strength to weight ratio and excellent properties such as water resistance and dielectric strength. It is widely used as the basic material in the packaging, piping and other fields. The visco-elastic behavior of HDPE has been intensively investigated in the long-term loading (tensile and compressive creep tests) by Elleuch and Taktak [START_REF] Elleuch | Viscoelastic behavior of HDPE polymer using tensile and compressive loading[J][END_REF] . Lai and Bakker [START_REF] Lai | Analysis of the non-linear creep of high-density polyethylene[END_REF] effectuated the tensile creep tests at different stress levels from 2 to 6 MPa and physical aging from 10 0 to 10 4 s.

The long-term creep behavior was predicted by using a non-linear creep model including the physical aging effects. Guo and Bradshaw [START_REF] Guo | Isothermal physical aging characterization of Polyether-ether-ketone (PEEK) and Polyphenylene sulfide (PPS) films by creep and stress relaxation[END_REF] studied the effect of transient creep and stressrelaxation of PEEK. Based on the HDPE and PEEK creep experimental data, Xu et al. [START_REF] Xu | Creep constitutive models for viscoelastic materials based on fractional derivatives[END_REF] proposed finite element, fractional Maxwell, fractional Kelvin-Voight and fractional Poynting-Thomson models to describe their visco-elastic behavior.

Fractional order has better memory than integer order, which permits to describe the polymeric material's visco-elastic behavior more accurately. Gement [START_REF] Gemant | On fractional differences[END_REF] firstly established a viscoelastic constitutive model using fractional derivatives. This rheological model can be used to well describe the mechanical properties with fewer experimental parameters.

A significant amount of work was effectuated to investigate the visco-elastic properties of the polymeric materials by using constitutive rheological models containing fractional derivatives [START_REF] Bagley | On the fractional calculus model of viscoelasticity behavior [J][END_REF] . For example, Leung et al. [START_REF] Leung | Steady state response of fractionally damped nonlinear viscoelastic arches by residue harmonic homotopy[END_REF] studied the non-linear steady-state vibrations of visco-elastic polymeric arches with fractional derivatives control. Lewandowski et al. [START_REF] Lewandowski | Dynamic analysis of frames with viscoelastic dampers modelled by rheological models with fractional derivatives[END_REF] explored the dynamics of polymeric damper frames based on the fractional derivative rheological model. Kang et al. [START_REF] Kang | A fractional non-linear creep model for coal considering damage effect and experimental validation[END_REF] established a fractional non-linear model to describe the creep behavior by taking into account the visco-elastic-plastic characteristics and damage effects. Gioacchino et al. [START_REF] Gioacchino | On the dynamics of non-local fractional viscoelastic beams under stochastic agencies[END_REF] used a fractionalorder state-variable expansion to establish a constitutive equation, and implement a numerical method to calculate the random response of a non-local fractional model under Gaussian white noise.

The visco-elastic behavior of polymeric materials and their constitutive rheological behavior modelling have gradually received the researchers' attention. A large number of fractional models have been intensively investigated to describe their visco-elastic properties. Lewandowski et al. [START_REF] Lewandowski | Nonlinear vibration of viscoelastic beams described using fractional order derivatives[END_REF] used a fractional Zener model to describe the non-linear vibration of a visco-elastic composite beam composed by the elastic and visco-elastic layers. The effect of harmonic force on non-linear vibration was investigated. Baum et al. [START_REF] Lewandowski | Dynamic characteristics of multilayered beams with viscoelastic layers described by the fractional Zener model[END_REF] applied a four-parameter rheological model with fractional derivatives to describe the mechanical properties of the multi-layer composite beam. This innovative method was validated in order to determine the dynamic properties of the composite beam. Cortes et al. [START_REF] Fernando | Finite element formulations for transient dynamic analysis in structural systems with viscoelastic treatments containing fractional derivative models[END_REF] used a five-parameter fractional derivative model to study the transient dynamics of a cantilever beam with different damping values under various loading conditions. Bahraini et al. [START_REF] Bahraini | Large deflection of viscoelastic beams using fractional derivative model[END_REF] established the fractional four-parameter derivative model to calculate the large deflection of the visco-elastic beams and the corresponding non-linear analysis was realized by finite element method. Paola et al. [START_REF] Di | Fractional visco-elastic Euler-Bernoulli beam[END_REF] examined the response of fractional visco-elastic Euler-Bernoulli beam under quasi-static and dynamic loads. He et al. [START_REF] He | Large amplitude vibration of fractionally damped viscoelastic CNTs/fiber/polymer multiscale composite beams[END_REF] used a fractionalorder time derivative damping model to study the large-scale free and forced vibration response of the carbon nanotubes/polymer laminated multi-scale composite beams.

The numerical algorithm to solve the vibration analysis equations based on these fractional rheological models with visco-elastic properties are well documented in the literature reviews. The determination of the numerical solution of the governing equations is more intractable, due to the complexity of the fractional derivatives. Laplace or Fourier transforms have been widely used to convert the time-domain problem into frequency-domain problem. For example, Lewandowski et al. [START_REF] Lewandowski | Dynamic characteristics of multilayered beams with viscoelastic layers described by the fractional Zener model[END_REF] used the virtual work principle and the Laplace transform to derive the frequency-domain motion equation for the dynamic analysis on the polymeric multi-layer beam. The other operational methods in fractional calculus including the finite element method, multi-scale method, Galerkin method and the variational iteration method are also found from the recent literature. Chang et al. [START_REF] Chang | Vibration and stability of an axially moving Rayleigh beam[END_REF] used the variable-domain finite element method to derive the motion equation to study the vibration and stability of the axial moving beam. Friswell et al. [START_REF] Friswell | Vibration analysis of beams with non-local foundations using the finite element method[END_REF] used the finite element method to analyze the dynamic characteristics of polymeric beams under different boundary conditions. Demir et al. [START_REF] Demir | Linear dynamical analysis of fractionally damped beams and rods[END_REF] established a general model in the dynamics of beams with fractional derivatives and solved the motion equation by multi-scale method. The conclusion is that the coefficient of the fractional derivative affects the stability of natural frequencies and vibration amplitudes. Permoon et al. [START_REF] Permoon | Application of radial basis functions and sinc method for solving the forced vibration of fractional viscoelastic beam[END_REF] used the Galerkin method to discretize the motion equation into a set of linear ordinary differential equations and then studied the forced vibration of beams. Martinez-Agirre et al. [START_REF] Martinez-Agirre | Higher order eigensensitivities-based numerical method for the harmonic analysis of viscoelastically damped structures[END_REF] studied the harmonic response of the constrained layer damped cantilever beam and analyzed the damping structure system by using the complex modal superposition method. Based on the constitutive relation in the form of genetic integral algorithm, Martin [START_REF] Martin | A modified variational iteration method for the analysis of viscoelastic beams[END_REF] established the mathematical model for the dynamic analysis of beams with viscoelastic properties by using Galerkin and variational iteration methods for both quasi-static and dynamic analysis.

The polynomial algorithm is quite suitable for solving the fractional equation because of its rapidity, high-efficiency and high-accuracy. The polynomial approximation method, such as Legendre polynomial method [START_REF] Meng | Numerical solutions of nonlinear fractional differential equations by alternative Legendre polynomials[END_REF][START_REF] Chen | Numerical solution of fractional partial differential equations with variable coefficients using generalized fractional-order Legendre functions[END_REF] , Chebyshev wavelet method [START_REF] Xie | A two-dimensional Chebyshev wavelets approach for solving the Fokker-Planck equations of time and space fractional derivatives type with variable coefficients[END_REF][START_REF] Chen | Numerical solution of nonlinear fractional integral differential equations by using the second kind Chebyshev wavelets[END_REF] , Bernoulli wavelet method [28] , Bernstein polynomial method [START_REF] Chen | Numerical study of a class of variable order nonlinear fractional differential equation in terms of Bernstein polynomials[END_REF][START_REF] Chen | Numerical solution for the variable order linear cable equation with Bernstein polynomials[END_REF] and shifted Chebyshev polynomials (SCPs) method [START_REF] Fernanda | Extrapolating for attaining high precision solutions for fractional partial differential equations[END_REF] are usually used to find the numerical solutions of the fractional equations. The SCPs method can be used to approximate the unknown function on the extended interval, which makes it easier to solve the fractional differential equations with different physical mechanisms governing and historical background. The SCPs is used as a basis to derive its integer order and fractional order differential operator matrices. The governing equations are transformed into a system of algebraic equations, which can be solved easily.

In this paper, the different fractional constitutive rheological models of the polymeric viscoelastic Euler-Bernoulli beam based on the fractional derivative element (FDE) model and the fractional derivative Kelvin-Voigt (FDKV) model have been established. These models have been implemented in the vibration governing equations of the visco-elastic beam, which are solved using specific fractional derivatives and especially SCPs approximation method. The numerical solutions are calculated directly form the governing equation in the time-space domain. The error analysis is effectuated and a correction method is proposed to obtain a high-resolution numerical solution. Euler-Bernoulli beams with different material properties for two specific thermoplastic polymers (PEEK and HDPE) have been studied as benchmark numerical examples. The numerical solutions concerning the displacement, the deformation and the stress of the beams have been obtained under quasi-static load and various boundary conditions. This paper is structured as follows, in Section 2, some preliminaries including the basic definition of fractional differential operators, the fractional rheological constitutive equations and the properties of SCPs are described. In Section 3, FDE and FDKV models have been used to establish the governing equation for different polymeric beam under different load and boundary conditions. In Section 4, the SCPs solving algorithm is described. An error analysis with the correction method is performed in Section 5. In Section 6, the numerical results of the displacement, the deformation and the stress of the Euler-Bernoulli beams are obtained and discussed to show the advantage of the proposed approach. The research work is concluded in Section 7.

Preliminaries and notations

The basic definition of fractional differential operator

Definition1. The Riemann-Liouville fractional differential operator ( 𝐷 𝛼 𝑅𝐿 𝑓) of order α is defined by [START_REF] Lin | Finite difference/spectral approximations for the time-fractional diffusion equation[END_REF] :

( 𝐷 𝛼 𝑓 𝑅𝐿 )(𝑦) = { 1 Γ(𝑚-𝛼) 𝑑 𝑚 𝑑𝑦 𝑚 ∫ 𝑓(𝑠) (𝑦-𝑠) 𝛼-𝑚+1 𝑑𝑠, α > 0, 𝑚 -1 ≤ α < 𝑚, 𝑦 0 𝑑 𝑚 𝑓(𝑦) 𝑑𝑦 𝑚 , 𝛼 = 𝑚. (1)
where 𝑚 ∈ 𝑁, 𝑓(𝑦) is continuous over interval (0, +∞) and is integrable over any subinterval [0, +∞), α is fractional derivative order, 𝑦 > 0, α and 𝑦 are real.

The gamma function, denoted by Γ(. ) is defined as Γ(𝑧) = ∫ 𝑒 -𝑦 𝑦 𝑧-1 𝑑𝑦 ∞ 0 for complex arguments with positive real part. Definition2. The Caputo definition of fractional differential operator ( 𝐷 𝛼 𝐶 𝑓) of order α is given by [START_REF] Taghvafard | Phase and anti-phase synchronization of fractional order chaotic systems via active control[END_REF] :

( 𝐷 𝛼 𝐶 𝑓)(𝑦) = { 1 Γ(𝑚-𝛼) ∫ 𝑓 (𝑚) (𝜏) (𝑦-𝜏) 𝛼-𝑚+1 𝑑𝜏, α > 0, 𝑚 -1 ≤ α < 𝑚, 𝑦 0 𝑑 𝑚 𝑓(𝑦) 𝑑𝑦 𝑚 , 𝛼 = 𝑚. (2)
The Caputo fractional derivative of distributed order is defined as:

𝐷 𝛼 𝐶 𝑦 𝜒 = { 0, 𝑓𝑜𝑟 𝜒 ∈ 𝑁 0 𝑎𝑛𝑑 𝜒 < ⌈𝛼⌉, 𝛤 (𝜒+1)
𝛤(𝜒+1-𝛼) 𝑦 𝜒-𝛼 , 𝑓𝑜𝑟 𝜒 ∈ 𝑁 0 𝑎𝑛𝑑 𝜒 ≥ ⌈𝛼⌉ 𝑜𝑟 𝜒 ∉ 𝑁 0 𝑎𝑛𝑑 𝜒 > ⌊𝛼⌋.

(3)

Properties of the SCPs

SCPs are considered as useful tools to solve the fractional equation in the physical problems. The well-known Chebyshev polynomials satisfy the following three-term recurrence relation:

𝐹 𝑖+1 (𝑐) = 2𝑐𝐹 𝑖 (𝑐) -𝐹 𝑖-1 (𝑐), 𝑖 = 1, 2, ⋯ (4)
where 𝐹 0 (𝑐) = 1 and 𝐹 1 (𝑐) = 𝑐. 𝑐 is defined on the interval [-1, 1] and 𝑖 = 1, 2, ⋯ On the interval 𝑚 ∈ [0, 𝐻], where 𝐻 is a non-negative real number, the SCPs are defined by the change of variable 𝑐 = 2𝑚 𝐻 -1. Let the SCPs 𝐹 𝑖 ( 2𝑚 𝐻 -1) denote by 𝐺 𝑖 (𝑚), which can be obtained as follows:

𝐺 𝑖+1 (𝑚) = 2 ( 2𝑚 𝐻 -1) 𝐺 𝑖 (𝑚) -𝐺 𝑖-1 (𝑚), 𝑖 = 1, 2, ⋯ (5)
where

𝐺 0 (𝑚) = 1, 𝐺 1 (𝑚) = 2𝑚 𝐻 -1.
The analytic form of 𝐺 𝑖 (𝑚) of i-degree is given by:

𝐺 𝑖 (𝑚) = 𝑖 ∑ (-1) 𝑖-𝑘 (𝑖+𝑘-1)!2 2𝑘 (𝑖-𝑘)!(2𝑘)!𝐻 𝑘 𝑖 𝑘=0 𝑚 𝑘 , 𝑖 = 1, 2, ⋯ (6)
where 𝐺 𝑖 (0) = (-1) 𝑖 and 𝐺 𝑖 (𝐻) = 1.

The orthogonally condition is

∫ 𝐺 𝑗 (𝑚) 𝐻 0
𝐺 𝑘 (𝑚)𝑊 𝐻 (𝑚)𝑑𝑚 = ℎ 𝑘 (7) where 𝑊 𝐻 (𝑚) = The operational matrix is defined by:

Φ 𝑛 (𝑚) = [𝐺 0 (𝑚), 𝐺 1 (𝑚), … , 𝐺 𝑛 (𝑚)] 𝑇
The following equation can be obtained:

Φ 𝑛 (𝑚) = 𝐴 𝑛 𝑍 𝑛 (𝑚) (8)
where 𝑍 𝑛 (𝑚) = [1, 𝑚, … , 𝑚 𝑛 ] 𝑇 , and 𝐴 𝑛 is the SCPs coefficient matrix given as follows:

𝐴 𝑛 = [ 𝑃 0,0 0 ⋯ 0 𝑃 1,0 𝑃 1,1 ⋯ 0 ⋮ ⋮ ⋱ ⋮ 𝑃 𝑛,0 𝑃 𝑛,1 ⋯ 𝑃 𝑛,𝑛 ] (9)
where

{ 𝑃 0,0 = 1, 𝑃 𝑖,𝑗 = 2 ( 2 𝐻 𝑃 𝑖-1,𝑗-1 -𝑃 𝑖-1,𝑗 ) -𝑃 𝑖-2,𝑗 𝑃 𝑖,𝑗 = 0, 𝑓𝑜𝑟 𝑖 < 𝑗 𝑜𝑟 𝑖 < 0 𝑜𝑟 𝑗 < 0.
, Obviously, 𝐴 𝑛 is full rank and reversible.

Definition of the visco-elastic constitutive relations based on fractional derivative models

The general form of a one-dimensional generalized fractional constitutive equation describing the stress-strain relationship of a visco-elastic material is [START_REF] Wilhelm | Viscoelasticity[END_REF] :

∑ 𝜇 𝜅 𝑑 𝑝 𝜅 𝜎 𝑑𝑡 𝑝 𝜅 𝑛 𝜅=0 = ∑ 𝜂 𝜅 𝑑 𝑞 𝜅 𝜀 𝑑𝑡 𝑞 𝜅 𝑛 𝜅=0 (10)
where 𝑑 𝜅 𝑑𝑡 𝜅 ⁄ uses the Caputo type fractional differential definition, 𝑡 is the time, 𝜎 is the stress, 𝜀 is the strain, and 𝜇 𝜅 , 𝜂 𝜅 are the material constants. 𝜅 is a positive integer. 𝑝 𝜅 , 𝑞 𝜅 are real numbers corresponding to fractional order of the time derivative.

Du et al. [START_REF] Du | Measuring memory with the order of fractional derivative[END_REF] proved that the physical meaning of fractional order is an indicator of memory, which is suitable to describe the memory phenomena in different disciplines. It is more accurately to describe the visco-elastic behavior of the polymeric material by using the fractional derivative constitutive equations.

When 𝑝 𝜅 , 𝑞 𝜅 = 0, the constitutive equation becomes ideal elastic behavior law or Hooke's law. When 𝑝 𝜅 = 0, 𝑞 𝜅 = 1 , the constitutive equation becomes ideal viscous behavior law or Newton's law. When 0 < 𝑞 𝜅 < 1, the constitutive equation could be used to describe the physical behavior of a viscoelastic material.

In the current visco-elastic constitutive fractional derivative models, the first order derivatives 𝑑 𝑑𝑡 ⁄ are replaced by the fractional derivatives 𝑑 𝑞 𝜅 𝑑𝑡 𝑞 𝜅 ⁄ with 0 < 𝑞 𝜅 < 1. The simplest visco-elastic fractional derivative model is the FDE model [START_REF] Xu | Creep constitutive models for viscoelastic materials based on fractional derivatives[END_REF] as shown in Fig. 1(a). The stress-strain relationship could be presented by the following equation:

σ = 𝜂 𝑑 𝛼 𝜀 𝑑𝑡 𝛼 (11)
where 𝜂, 𝛼 represent the material constants, 𝛼 is a real number within (0, 1). Here, 𝑑 𝛼 𝑑𝑡 𝛼 is the Riemann-Liouville fractional differential operator of order 𝛼.

The two-dimensional FDE model can be expressed as:

𝜎(𝑥, 𝑡) = 𝜂 𝐷 𝛼 𝐶 𝜀(𝑥, 𝑡) (12)
where 𝑥 is the position and 𝐷 𝛼 𝐶 is Caputo fractional differential operator, defined in Eq. ( 2). When two fractional derivative elements are arranged in parallel, as shown in Fig. 1(b), it is the FDKV model [START_REF] Xu | Creep constitutive models for viscoelastic materials based on fractional derivatives[END_REF] , in which the constitutive equation is defined as follows:

σ = 𝜂 1 𝑑 𝛽 𝜀 𝑑𝑡 𝛽 + 𝜂 2 𝑑 𝛾 𝜀 𝑑𝑡 𝛾 (13)
where 𝜂 1 , 𝜂 2 , 𝛽, 𝛾 are the material constants, and 𝛽, 𝛾 are the real number within (0, 1).

The two-dimensional FDKV model can be expressed as:

𝜎(𝑥, 𝑡) = 𝜂 1 𝐷 𝛽 𝐶 𝜀(𝑥, 𝑡) + 𝜂 2 𝐷 𝛾 𝐶 𝜀(𝑥, 𝑡) (14)
where

𝐷 𝛽 𝐶 , 𝐷 𝛾 𝐶
is Caputo fractional differential operator.

The governing equations of the visco-elastic Euler-Bernoulli beam

Beam governing equation with the FDKV model

A visco-elastic Euler-Bernoulli beam, fixed at both ends, is considered in this study. A distributed load is applied on the vertical direction of the beam. The beam is treated as a viscoelastic material. The bending deformation occurs on the beam in the vertical direction, as shown in Fig. 2, in which 𝜔(𝑥, 𝑡) is the beam deflection, 𝑓(𝑥, 𝑡) is the distributed load, 𝑙 is the length of the beam, 𝑏 is the width of a rectangular cross section and ℎ is the height of a rectangular cross section.

FIGURE 2

The bending deformation of the visco-elastic Euler-Bernoulli beam under the distributed load.

Based on Fig. 2, the most important displacement of the beam is the deflection. The beam deflection varies in function of 𝑥 and 𝑡. The differential equation of the transverse oscillation has been studied by Martin [START_REF] Martin | Stability approach to the fractional variational iteration method used for the dynamic analysis of viscoelastic beams[END_REF] and the integral differential equation of the beam is given by:

𝜌𝐴 𝜕 2 𝜔(𝑥,𝑡) 𝜕𝑡 2 + 𝜕 2 𝑀(𝑥,𝑡) 𝜕𝑥 2 = 𝑓(𝑥, 𝑡) ( 15
)
where 𝜌 is the density of the material and 𝐴 is the square cross-sectional area. The bending moment 𝑀(𝑥, 𝑡) of the beam is written as:

𝑀(𝑥, 𝑡) = ∫ 𝑧 𝐴 𝜎(𝑥, 𝑡)𝑑𝐴 (16)
where 𝜎(𝑥, 𝑡) is the normal stress on the cross-section and z represent the transverse coordinate.

Based on the Euler-Bernoulli beam theory, the relation between strain and displacement can be expressed as:

𝜀(𝑥, 𝑡) = 𝑧 𝜕 2 𝜔(𝑥,𝑡) 𝜕𝑥 2 (17)
When the FDKV model is used to describe the visco-elastic behavior of the beam, the stress of the beam can be obtained by using Eqs. ( 14) and (17):

𝜎(𝑥, 𝑡) = 𝜂 1 𝐷 𝛽 𝐶 𝑧 𝜕 2 𝜔(𝑥,𝑡) 𝜕𝑥 2 + 𝜂 2 𝐷 𝛾 𝐶 𝑧 𝜕 2 𝜔(𝑥,𝑡) 𝜕𝑥 2 (18)
Based on Eqs ( 16) and ( 18), the beam bending moment 𝑀(𝑥, 𝑡) is equal to:

𝑀(𝑥, 𝑡) = 𝜂 1 𝐼 𝐷 𝛽 𝐶 𝜕 2 𝜔(𝑥,𝑡) 𝜕𝑥 2 + 𝜂 2 𝐼 𝐷 𝛾 𝐶 𝜕 2 𝜔(𝑥,𝑡) 𝜕𝑥 2 (19)
where the moment inertia is 𝐼 = ∫ 𝑧

Beam governing equation under the FDE model

When the FDE model is used to describe the visco-elastic behavior of the beam, the stress of the beam is proposed:

𝜎(𝑥, 𝑡) = 𝜂 𝐷 𝛼 𝐶 𝑧 𝜕 2 𝜔(𝑥,𝑡) 𝜕𝑥 2 (21)
where 𝛼 ∈ (0,1).

Based on Eqs. ( 16) and ( 21), the beam bending moment 𝑀(𝑥, 𝑡) can be rewritten as:

𝑀(𝑥, 𝑡) = 𝜂𝐼 𝐷 𝛼 𝐶 𝜕 2 𝜔(𝑥,𝑡) 𝜕𝑥 2 (22)
The bending vibration equation of the visco-elastic Euler-Bernoulli beam with FDE model is obtained as follows:

𝜌𝐴 𝜕 2 𝜔(𝑥,𝑡) 𝜕𝑡 2 + 𝜂𝐼 𝐷 𝛼 𝐶 𝜕 4 𝜔(𝑥,𝑡) 𝜕𝑥 4 = 𝑓(𝑥, 𝑡) (23)

Numerical algorithm

The approximation function based on the families of SCPs of 𝑥 and 𝑡 was applied to replace the real displacement function. The collocation method was used to discretize the variables 𝑥 and 𝑡 to transform the fractional beam governing equations into a set of algebraic equations. An error analysis was preformed to estimate the calculation accuracy.

The most important contribution of this numerical algorithm is that the numerical solution of the beam governing equation is successfully obtained with high accuracy and the obtained displacement and stress of beam is directly in time domain. The numerical algorithm is summarized in Fig. 3: FIGURE 3 A schematic illustration of the proposed numerical algorithm.

Approximation of the beam deflection function

The beam deflection function 𝜔(𝑥, 𝑡) ∈ 𝐿 2 ([0, 𝐻] × [0, 𝑇]) can be expanded as the following formula:

𝜔(𝑥, 𝑡) = ∑ ∑ 𝜔 𝑖𝑗 𝐺 𝑖 (𝑥)𝐺 𝑗 (𝑡) ∞ 𝑗=0 ∞ 𝑖=0 (24)
where

𝜔 𝑖𝑗 = 1 ℎ 𝑖 ℎ 𝑗 ∫ ∫ 𝜔(𝑥, 𝑡)𝐺 𝑖 (𝑥)𝐺 𝑗 (𝑡)𝜔 𝐻 (𝑥)𝜔 𝑇 (𝑡)𝑑𝑥𝑑𝑡, 𝑖, 𝑗 𝑇 0 𝐻 0 = 0, 1, 2, ⋯.
If the truncated series are considered in Eq. ( 24), then it can be rewritten as:

𝜔(𝑥, 𝑡) ≈ ∑ ∑ 𝜔 𝑖𝑗 𝐺 𝑖 (𝑥)𝐺 𝑗 (𝑡) 𝑛 𝑗=0 𝑛 𝑖=0 = Φ 𝑛 𝑇 (𝑥)𝑈Φ 𝑛 (𝑡) (25)
where

Φ 𝑛 (𝑥) = [𝐺 0 (𝑥), 𝐺 1 (𝑥), … , 𝐺 𝑛 (𝑥)] 𝑇 , Φ 𝑛 (𝑡) = [𝐺 0 (𝑡), 𝐺 1 (𝑡), … , 𝐺 𝑛 (𝑡)] 𝑇 , 𝑈 = {𝑢 𝑖𝑗 } 𝑖,𝑗=0 𝑛,𝑛 .

SCPs differential operator matrix of displacement function

First order differential operator matrix of SCPs

The first order derivative of Φ 𝑛 (𝑥) for 𝑥 can be expressed as follows by using Eq. ( 8):

𝑑 𝑑𝑥 Φ 𝑛 (𝑥) = (𝐴 𝑋 𝑍 𝑛 (𝑥)) ′ = 𝐴 𝑋 𝑍 𝑛 ′ (𝑥) = 𝐴 𝑋 [ 1 ′ 𝑥 ′ ⋮ (𝑥 𝑛 ) ′ ] = 𝐴 𝑋 [ 0 1 ⋮ 𝑛𝑥 𝑛-1 ] = 𝐴 𝑋 𝑉 (𝑛+1)×(𝑛+1) 𝑍 𝑛 (𝑥) (26)
where 𝐴 𝑋 is obtained by replacing 𝑛 by 𝑋 in Eq. ( 9), and 𝑉 (𝑛+1)×(𝑛+1) = [ 0 0 0 ⋯ 0 0 1 0 0 ⋯ 0 0 0 2 0 ⋯ 0 0 ⋮ ⋮ ⋮ ⋱ ⋮ ⋮ 0 0 0 ⋯ 0 0 0 0 0 ⋯ 𝑛 0 ] 𝑍 𝑛 (𝑥) is obtained by using Eq. ( 8):

𝑍 𝑛 (𝑥) = 𝐴 𝑋 -1 Φ 𝑛 (𝑥) (27)
Based on Eqs. ( 26) and ( 27), the first order derivative of Φ 𝑛 (𝑥) for 𝑥 can be rewritten as:

𝑑 𝑑𝑥 Φ 𝑛 (𝑥) = 𝐴 𝑋 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑋 -1 Φ 𝑛 (𝑥) = 𝑃 𝑥 1 Φ 𝑛 (𝑥) (28)
where 𝑃 𝑥 1 = 𝐴 𝑋 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑋 -1 is called the first order differential operator matrix of the SCPs for 𝑥.

Based on Eqs. ( 25) and ( 28), the following equations could be obtained: where 𝐴 𝑇 is obtained by replacing 𝑛 by 𝑇 in Eq. ( 9), and 𝑃 𝑡 1 = 𝐴 𝑇 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑇 -1 is the first order differential operator matrix of the SCPs for 𝑡.

𝜕𝜔(𝑥,𝑡) 𝜕𝑥 ≈ 𝜕(Φ 𝑛 𝑇 (𝑥)𝑈Φ 𝑛 (𝑡)) 𝜕𝑥 = Φ 𝑛 𝑇 (𝑥)(𝑃 𝑥 1 ) 𝑇 𝑈Φ 𝑛 (𝑡) = Φ 𝑛 𝑇 (𝑥)(𝐴 𝑋 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑋 -1 ) 𝑇 𝑈Φ 𝑛 (𝑡) (29) 𝜕 4 𝜔(𝑥,𝑡) 𝜕𝑥 4 ≈ 𝜕 4 (Φ 𝑛 𝑇 (𝑥)𝑈Φ 𝑛 (𝑡)) 𝜕𝑥 4 = Φ 𝑛 𝑇 (𝑥)((𝑃 𝑥 1 ) 𝑇 ) 4 𝑈Φ 𝑛 (𝑡) = Φ 𝑛 𝑇 (𝑥) ((𝐴 𝑋 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑋 -1 ) 𝑇 ) 4 𝑈Φ 𝑛 (𝑡) (30)

Fractional order differential operator matrix of the SCPs

The fractional derivative of Φ 𝑛 (𝑡) for 𝑡 could be calculated by the following equation according to Eq. ( 25):

𝐷 𝛼 𝐶 Φ 𝑛 (𝑡) = 𝐷 𝛼 𝐶 (𝐴 𝑇 𝑍 𝑛 (𝑡)) = 𝐴 𝑇 𝐷 𝛼 𝑍 𝑛 (𝑡) = 𝐶 𝐴 𝑇 [ 0 Γ(2) Γ(2-𝛼) 𝑡 1-𝛼 ⋮ Γ(𝑛+1) Γ(𝑛+1-𝛼) 𝑡 𝑛-𝛼 ] = 𝐴 𝑇 𝑉 (𝑛+1)×(𝑛+1) 𝛼 𝐴 𝑇 -1 Φ 𝑛 (𝑡) = 𝑃 𝑡 𝛼 Φ 𝑛 (𝑡) (32)
where 𝑃 𝑡 𝛼 = 𝐴 𝑇 𝑉 (𝑛+1)×(𝑛+1) 𝛼 𝐴 𝑇 -1 is called the fractional order differential operator matrix of the SCPs for 𝑡 and

V (n+1)×(n+1) α = [ 0 0 ⋯ 0 0 Γ(2) Γ(2-α) t -α ⋯ 0 ⋮ ⋮ ⋱ ⋮ 0 0 ⋯ Γ(n+1) Γ(n+1-α) t -α ] (33)
Based on Eq. ( 32), the following equations could be found:

where

𝑉 (𝑛+1)×(𝑛+1)
𝛾 is obtained by replacing 𝛼 by 𝛾 in Eq. ( 33).

The bending vibration equation of the visco-elastic Euler-Bernoulli beam with FDKV model, Eq. ( 20), could be transformed into the following form:

𝜌𝐴Φ 𝑛 𝑇 (𝑥)𝑈(𝐴 𝑇 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑇 -1 ) 2 Φ 𝑛 (𝑡) + 𝜂 1 𝐼Φ 𝑛 𝑇 (𝑥) ((𝐴 𝑋 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑋 -1 ) 𝑇 ) 4 𝑈𝐴 𝑇 𝑉 (𝑛+1)×(𝑛+1) 𝛽 × 𝐴 𝑇 -1 Φ 𝑛 (𝑡) + 𝜂 2 𝐼Φ 𝑛 𝑇 (𝑥) ((𝐴 𝑋 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑋 -1 ) 𝑇 ) 4 𝑈𝐴 𝑇 𝑉 (𝑛+1)×(𝑛+1) 𝛾 𝐴 𝑇 -1 Φ 𝑛 (𝑡) = 𝑓(𝑥, 𝑡) (36)
Similarly, the bending vibration equation of the visco-elastic Euler-Bernoulli beam with FDE model, Eq. ( 23), could be transformed into the following form:

Based on the collocation method, the reasonable match points 𝑥 𝑖 = 2𝑖-1 2(𝑛+1) 𝐻, 𝑖 = 0, 1, 2, ⋯ , 𝑛, 𝑡 𝑗 = 2𝑗-1 2(𝑛+1) 𝑇, 𝑗 = 0, 1, 2, ⋯ , 𝑛. have been used to discretize the variable (𝑥, 𝑡) to (𝑥 𝑖 , 𝑡 𝑗 ). Eqs. ( 36) and ( 37) are transformed into a set of algebraic equations. The coefficient 𝜔 𝑖𝑗 ( 𝑖 = 0, 1, 2, ⋯ , 𝑛; 𝑗 = 0, 1, 2, ⋯ , 𝑛) is determinate by using Matlab platform and least square method. The numerical solution of the fractional derivative equations can be obtained.

Error Analysis

Error correction

The equation to be solved is defined as the operator 𝐿[ * ]:

𝐿[𝜔(𝑥, 𝑡)] = 0 (38)
where 𝜔(𝑥, 𝑡) is the exact solution of the beam deflection.

where 𝑒 𝑛 (𝑥, 𝑡) can be approximated by 𝑒 * (𝑥, 𝑡) , which is calculated by using the algorithm proposed in the previous section. 𝐿[𝑒 𝑛 (𝑥, 𝑡)] is the approximate error function in Eq. ( 44). Therefore, the corrected solution 𝜔 * (𝑥, 𝑡) can be obtained as follows:

𝜔 * (𝑥, 𝑡) = 𝜔 𝑛 (𝑥, 𝑡) + 𝑒 * (𝑥, 𝑡)

Furthermore, correction error function 𝐸(𝑥, 𝑡) can be expressed as:

𝐸(𝑥, 𝑡) = 𝑒 𝑛 (𝑥, 𝑡) -𝑒 * (𝑥, 𝑡) = 𝜔(𝑥, 𝑡) -𝜔 𝑛 (𝑥, 𝑡) -𝑒 * (𝑥, 𝑡) (46)

Absolute error bound of correction solution method

Theorem1.

Let 𝜔(𝑥, 𝑡) and 𝜔 * (𝑥, 𝑡) are the exact and correction solutions of the governing equations Eq. ( 20 where 𝜔 𝑇,𝑀 (𝑥, 𝑡) is 𝑀 𝑡ℎ order Taylor series expansion of 𝜔(𝑥, 𝑡) in the neighborhood of (𝑥 0 , 𝑡 0 ), and 𝑅 𝑀 represents the remainder of the Taylor series expansion.

Proof: If the function 𝜔(𝑥, 𝑡) is continuous in a neighborhood around (𝑥 0 , 𝑡 0 ), and (𝑥 0 + ℎ, 𝑡 0 + 𝑘) is any point in the neighborhood. Since 𝜔(𝑥, 𝑡) has an 𝑀 + 1 order differential, 𝜔(𝑥, 𝑡) can be expressed as a Taylor series expansion:

𝜔(𝑥, 𝑡) = ∑ (ℎ 𝜕 𝜕𝑥 + 𝑘 𝜕 𝜕𝑡 ) 𝑛 𝜔(𝑥 0 , 𝑡 0 ) 𝑀 𝑛=0 + 𝑅 𝑀 (48)
where 𝑅 𝑀 is the remainder of the Taylor series expansion:

𝑅 𝑀 = 1 (𝑀+1)! (ℎ 𝜕 𝜕𝑥 + 𝑘 𝜕 𝜕𝑡 ) 𝑀+1 𝜔(𝑥 0 + 𝜃ℎ, 𝑡 0 + 𝜃𝑘), 0 < θ < 1 (49)
Furthermore, the 𝑀 𝑡ℎ Taylor series expansion 𝜔 𝑇,𝑀 (𝑥, 𝑡) of 𝜔(𝑥, 𝑡) in the neighborhood of (𝑥 0 , 𝑡 0 ) is obtained:

𝜔(𝑥,

Numerical analysis 6.1 Mathematical example

A general mathematical example is proposed, associated with its exact solution. The general mathematical example is solved by the proposed SCPs numerical algorithm. The numerical solution is compared with the exact solution to verify the accuracy and efficiency of the proposed algorithm. The general mathematical example is similar with the governing equation with FDE model, therefore the parameters in the example are arbitrary values and have no physical meanings. The general mathematical example is given by: 𝑡 1.5 .

The boundary and initial conditions are:

{ 𝜔(0, 𝑡) = 𝜔(2, 𝑡) = 0 𝜔 ′ (0, 𝑡) = 𝜔 ′ (2, 𝑡) = 0 𝜔(𝑥, 0) = 𝜔 ′ (𝑥, 0) = 0 (53)
The exact solution of the Eq. ( 52) is 𝜔(𝑥, 𝑡) = 𝑥 2 (2 -𝑥) 2 𝑡 2 . The preceding proposed method, defined in paragraph 4.2.2, was applied to find the solution of Eq. (52). The SCPs numerical algorithm with the number of terms equal to 4 (𝑛 = 4) was used to calculate the numerical solution 𝜔 𝑛 (𝑥, 𝑡) of Eq. ( 52). The comparison between its numerical solution and its exact solution for different 𝑡 and 𝑥 are shown in Fig. 4. The numerical results seem to be consistent with the exact solutions, which proves a high accuracy of the proposed algorithm in the approximation of the numerical solution of the fractional derivative equations. The absolute error function, defined in Eq. ( 42), is used to calculate the difference between the numerical and exact solutions. The absolute error at different t: t = 0, t = 0.2, t = 0.4, t = 0.6, 𝑡 = 0.8, 𝑡 = 1 are shown in Fig. 5. The absolute errors are around 10 -14 , which confirm that a high-precision has been achieved by using the proposed algorithm. The algorithm can effectively solve the fractional partial differential equations. The number of terms of the SCPs may affect the calculate accuracy during the approximation of the exact solution of the fractional partial differential equations. Different number of terms of the SCPs 𝑛 = 4, 𝑛 = 5, 𝑛 = 6 has been proposed to solve the fractional equations to investigate its influence on the numerical calculation efficiency. The absolute and relative errors of the numerical and exact solutions at same (𝑥, 𝑡) with the terms of the SCPs equal to 4, 5 and 6 are summarized in Table 1. It shows that a high-precision of numerical calculation is obtained even with a small number of terms of SCPs. The absolute error between the numerical and the exact solutions decreases when the number of terms of the SCPs increases. When 𝑛 = 6, the absolute error can reach 10 -17 . The concept of relative error is introduced to describe the precision between the numerical and the exact solutions, the accuracy performance of the algorithm is shown in this numerical example. .4910 -13 4.6810 -16 2.8610 -14 3.1210 -17 1.9110 -15 1.0 5.9310 -15 1.4810 -13 3.0510 -15 7.6310 -14 1.5810 -15 3.9610 -14 1. 6 3.8910 -15 2.3710 -13 3.2810 -15 2.0010 -13 1.8910 -15 1.1510 -13 0.6 0.4 2.4610 -15 1.6710 -14 1.8810 -15 1.2810 -14 1.7610 -15 1.1910 -14 1.0 1.9510 -14 5.4310 -14 2.2210 -15 6.1710 -15 2.1110 -15 5.8610 -15 1.6

3.7110 -14 2.5110 -13 1.3310 -15 9.0310 -15 8.8810 -16 6.0210 -15

1.0 0.4 1.4410 -14 3.5310 -14 1.2310 -16 2.9910 -16 7.4710 -16 1.8210 -15 1.0 5.0810 -15 5.0810 -15 1.7810 -15 1.7810 -15 5.2110 -16 5.2110 -16 1.6

1.0810 -14 2.6310 -14 1.6610 - Based on the Section 5, the correction error between the corrected solution and the exact solution can be obtained. Fig. 6 shows the correction error is smaller than the absolute error between the numerical solution and the exact solution. That is, the correction solution is closer to the exact solution. When the number of terms of the SCPs 𝑛 = 4, the error precision can be improved from 10 -14 to 10 -16 . The proposed correction solution improves the accuracy of the algorithm.

FIGURE 7

The CPU time of the numerical solution of the equation is obtained by the SCPs method when 𝑛 takes a different number.

The calculate time with different number of terms of the SCPs is shown in Fig. 7. An i5-8250U processor with the ram of 8GB has been used to effectuate these calculates. Fig 7 shows that the calculate time increase with the increase of the number of terms of SCPs. Therefore, the numerical solution can be efficiently approximated in a relatively short time. The proposed numerical algorithm shows excellent calculating efficiency and could solve the fractional partial differential equations quickly.

Numerical solution of the deformation of visco-elastic beam

The material parameters in FDE model for two visco-elastic polymeric materials (HDPE and PEEK) are defined in Table 2 and 3 [START_REF] Xu | Creep constitutive models for viscoelastic materials based on fractional derivatives[END_REF] : The length of the beam is 𝑙 = 5 m and its cross-section area is 𝐴 = 0.04 m 2 . The moment of inertia of the beam is 𝐼 = (0.2) [START_REF] Guo | Isothermal physical aging characterization of Polyether-ether-ketone (PEEK) and Polyphenylene sulfide (PPS) films by creep and stress relaxation[END_REF] 12 m 4 . The beam is fixed at both ends, so the boundary conditions are as follows:

𝜔(𝑥, 𝑡) = 0, 𝜕𝜔(𝑥,𝑡) 𝜕𝑥 = 0, 𝑥 = 0, 𝑙 (55)
The initial conditions of the system are:

𝜔(𝑥, 0) = 0, 𝜕𝜔(𝑥,𝑡) 𝜕𝑥 = 0 (56)

HDPE beam with FDKV model

In this section, the governing equation with FDKV model was used to describe the viscoelastic behavior of the HDPE beam. The deflections, strain and stress of HDPE beam are obtained by using the proposed SCPs algorithm. The deflections of the beam under different load conditions (Heaviside in cases Based on Eq. ( 18), the stress of the beam under different load conditions (Heaviside, harmonics and linear loads) is illustrated in Fig. 9. The stress increases with the loading time for all kinds of load conditions. At the same loading time, the value of the stress is minimum when near 𝑥 = 2.5 m, and the force resisting external load is the weakest, and the displacement is the smallest. The value of the stress increase when the position is close to the fixed end of the beam, and the greater the force against external load, the smaller the displacement. The deformation of the beam is symmetric with respect to the axis 𝑥 = 2.5 m, which is similar with the observation in the deflection. It can be seen that the values of deflection and stress are successfully obtained by using the proposed numerical algorithm. The numerical results are consistent with the experimental investigations of the beam vibration deformation [START_REF] Martin | Stability approach to the fractional variational iteration method used for the dynamic analysis of viscoelastic beams[END_REF][START_REF] Yu | A numerical method for solving fractionalorder viscoelastic Euler-Bernoulli beams[END_REF] .

Comparison of HDPE and PEEK beams

In this section, the governing equation with FDE model was used to describe the visco-elastic behavior of the HDPE and PEEK beam. The material parameters of the polymers in the FDE model were put into the Eq. ( 23), which is solved by using the proposed numerical algorithm. The deflection of the HDPE and PEEK beam in function of 𝑥 under different load conditions is shown in Fig. 10. Based on Fig. 10, the deflection of the PEEK beam is less than that of the HDPE beam under the same load condition. The difference of the deflection of these two polymers becomes more important with the loading time. One can conclude that the PEEK exhibits better bending resistance than the HDPE due to its smaller deflection under the same load conditions. This is consistent with the material properties, because the elastic modulus of PEEK is more important than that of HDPE in the same thermal and mechanical conditions. The proposed numerical algorithm seems to be efficient to predict the vibration deformation of the visco-elastic polymers under various load conditions.

Comparison of FDE and FDKV models

In this section, the governing equations with FDE and FDKV model were used to describe the visco-elastic behavior of the HDPE beam. The objective of this work is to investigate the sensibility of the visco-elastic models in the prediction of beam deflection under different load conditions. The deflection of the HDPE beam, which is calculated by the governing equation based on FDE model and FDKV model is shown in Based on Fig. 11, the deflection of the HDPE beam with the FDE model is slightly more important than that with the FDKV model under the same load, especially in the positions far from the fixed ends. But the difference between the values of defection with these two models is quite small relative to the length of the beam. One can conclude that the displacements of the beam calculated by using the governing equations with FDE and FDKV models are approximately coincide under the same load, which confirms the efficiency and accuracy of the proposed SCPs method in the numerical solution of the fractional partial differential equations.

Conclusion

In this paper, the governing equations of visco-elastic Euler-Bernoulli beams have been created by using the fractional rheological constitutive models. An effective numerical algorithm for solving the governing equations is proposed in the time-domain. The fractional rheological model containing the fractional derivatives is used to analyze the inherent laws of the dynamic performance of visco-elastic damping materials, which can provide a theoretical basis for the research, development and performance prediction of vibration damping materials. Based on the visco-elasticity of the FDE and FDKV constitutive models, the SCPs operator matrix is deduced and the numerical solution of the governing equations of the visco-elastic beam under different loads is obtained. The displacement solutions of the HDPE beam based on two constitutive models are compared. Different visco-elastic beams have been studied to verify the calculation accuracy and efficiency of the algorithm. 1) Two governing equations based on the FDE and FDKV models are established to analyze the deformation of the visco-elastic Euler-Bernoulli beam. The stress and strain of the beam are obtained by using the proposed algorithm.

2) The proposed method has been validated by the numerical solution of the governing equation. Based on the comparison of the numerical solutions with the FDKV and FDE model, the deflections of the beam are almost completely coincident at different loading time, which permits to verify the accuracy of the algorithm.

3) The displacement solutions of the visco-elastic beam of HDPE and PEEK under uniform load, harmonic load and linear load are compared to investigate the material sensibility of the proposed algorithm. The conclusion is that the displacement of PEEK beam is smaller than that of HDPE beam, when the beam is under the same load condition. The greater the damping of the corresponding visco-elastic material and the better the bending resistance.

Declaration of interests

√ 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.

☐ The authors declare the following financial interests/personal relationships which may be considered as potential competing interests:

There are no conflicts of interest.

*Declaration of Interest Statement

SAVE ICMJE Form for Disclosure of Potential Conflicts of Interest

Chen 1

The purpose of this form is to provide readers of your manuscript with information about your other interests that could influence how they receive and understand your work. The form is designed to be completed electronically and stored electronically. It contains programming that allows appropriate data display. Each author should submit a separate form and is responsible for the accuracy and completeness of the submitted information. The form is in six parts.

1.

2.

3.

4.

5.

Identifying information.

The work under consideration for publication.

This section asks for information about the work that you have submitted for publication. The time frame for this reporting is that of the work itself, from the initial conception and planning to the present. The requested information is about resources that you received, either directly or indirectly (via your institution), to enable you to complete the work. Checking "No" means that you did the work without receiving any financial support from any third party --that is, the work was supported by funds from the same institution that pays your salary and that institution did not receive third-party funds with which to pay you. If you or your institution received funds from a third party to support the work, such as a government granting agency, charitable foundation or commercial sponsor, check "Yes".

Relevant financial activities outside the submitted work.

This section asks about your financial relationships with entities in the bio-medical arena that could be perceived to influence, or that give the appearance of potentially influencing, what you wrote in the submitted work. You should disclose interactions with ANY entity that could be considered broadly relevant to the work. For example, if your article is about testing an epidermal growth factor receptor (EGFR) antagonist in lung cancer, you should report all associations with entities pursuing diagnostic or therapeutic strategies in cancer in general, not just in the area of EGFR or lung cancer.

Report all sources of revenue paid (or promised to be paid) directly to you or your institution on your behalf over the 36 months prior to submission of the work. This should include all monies from sources with relevance to the submitted work, not just monies from the entity that sponsored the research. Please note that your interactions with the work's sponsor that are outside the submitted work should also be listed here. If there is any question, it is usually better to disclose a relationship than not to do so.

For grants you have received for work outside the submitted work, you should disclose support ONLY from entities that could b e perceived to be affected financially by the published work, such as drug companies, or foundations supported by entities that could be perceived to have a financial stake in the outcome. Public funding sources, such as government agencies, charitable foundatio ns or academic institutions, need not be disclosed. For example, if a government agency sponsored a study in which you have been involved and drugs were provided by a pharmaceutical company, you need only list the pharmaceutical company.

Intellectual Property.

This section asks about patents and copyrights, whether pending, issued, licensed and/or receiving royalties.

Relationships not covered above.

Use this section to report other relationships or activities that readers could perceive to have influenced, or that give the appearance of potentially influencing, what you wrote in the submitted work.

1 √𝐻𝑚

 1

FIGURE 1

 1 FIGURE 1 Schematic representation of visco-elastic model: (a) FDE model and (b) FDKV model.

  𝑇 (𝑥)𝑈(𝑃 𝑡 1 ) 2 Φ 𝑛 (𝑡) = Φ 𝑛 𝑇 (𝑥)𝑈(𝐴 𝑇 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑇 -1 ) 2 Φ 𝑛 (𝑡)(31)

  𝜌𝐴Φ 𝑛 𝑇 (𝑥)𝑈(𝐴 𝑇 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑇 -1 ) 2 Φ 𝑛 (𝑡) + 𝜂𝐼Φ 𝑛 𝑇 (𝑥) ((𝐴 𝑋 𝑉 (𝑛+1)×(𝑛+1) 𝐴 𝑋 -𝑛 (𝑡) = 𝑓(𝑥, 𝑡)

  ). If the continuous function 𝜔(𝑥, 𝑡) has 𝐾 order partial derivative on the interval Λ ∈ [0, 𝑋] × [0, 𝑇], where 𝐾 = 0, 1, ⋯ , 𝑀, then: |𝜔(𝑥, 𝑡) -𝜔 * (𝑥, 𝑡)| ≤ |𝑅 𝑀 | + |𝜔 𝑇,𝑀 (𝑥, 𝑡) -𝜔 * (𝑥, 𝑡)| (47)

  𝑡) -𝜔 𝑇,𝑀 (𝑥, 𝑡) = 𝑅 𝑀 (50) Based on Eq. (50) and the triangle inequality, the absolute error bound of the corrected solution could be obtained: |𝜔(𝑥, 𝑡) -𝜔 * (𝑥, 𝑡)| = |𝜔(𝑥, 𝑡) -𝜔 * (𝑥, 𝑡) + 𝜔 𝑇,𝑀 (𝑥, 𝑡) -𝜔 𝑇,𝑀 (𝑥, 𝑡)| ≤ |𝜔(𝑥, 𝑡) -𝜔 𝑇,𝑀 (𝑥, 𝑡)| + |𝜔 𝑇,𝑀 (𝑥, 𝑡) -𝜔 * (𝑥, 𝑡)| = |𝑅 𝑀 | + |𝜔 𝑇,𝑀 (𝑥, 𝑡) -𝜔 * (𝑥, 𝑡)| (51)

FIGURE 4

 4 Comparison of exact and numerical solutions of 𝜔(𝑥, 𝑡) in the equation for different values: (a) 𝑡; (b) 𝑥.

FIGURE 5

 5 The evolution of absolute error between the numerical solution of the equation when 𝑡 takes different values: (a) 𝑡 = 0; (b) 𝑡 = 0.2; (c) 𝑡 = 0.4; (d) 𝑡 = 0.6; (e) 𝑡 = 0.8; (f) 𝑡 = 1. Since the relative error can better explain the validity of the method, the relative error 𝑒(𝑥, 𝑡) of the equation is introduced when solving the example equation. The relative error 𝑒(𝑥, 𝑡) can be expressed as: 𝑒(𝑥, 𝑡) = | 𝜔 𝑛 (𝑥,𝑡)-𝜔(𝑥,𝑡)

FIGURE 6

 6 FIGURE 6 When 𝑛 = 4, the absolute and correction errors between the numerical solution and the exact solution of the equation solved by SCPs method are compared: (a) Absolute error 𝑒 𝑛 (𝑥, 𝑡); (b) Correction error 𝑒 * (𝑥, 𝑡).

FIGURE 8 8 .FIGURE 9

 889 FIGURE 8 Numerical solutions of deflection of the HDPE beam under different load conditions: (a) 𝑓(𝑥, 𝑡) = 10𝐻𝑒𝑎𝑣𝑖𝑠𝑖𝑑𝑒(𝑡) ; (b) 𝑓(𝑥, 𝑡) = 30𝐻𝑒𝑎𝑣𝑖𝑠𝑖𝑑𝑒(𝑡) ; (c) 𝑓(𝑥, 𝑡) = 10𝜋 sin 𝑡 ; (d) 𝑓(𝑥, 𝑡) = 10 + 2𝑥.

FIGURE 10

 10 FIGURE 10 Comparison of the deflection of the HDPE and PEEK beams under different load conditions: (a) 𝑓(𝑥, 𝑡) = 10𝐻𝑒𝑎𝑣𝑖𝑠𝑖𝑑𝑒(𝑡) ; (b) 𝑓(𝑥, 𝑡) = 30𝐻𝑒𝑎𝑣𝑖𝑠𝑖𝑑𝑒(𝑡) ; (c) 𝑓(𝑥, 𝑡) = 10𝜋 sin 𝑡; (d) 𝑓(𝑥, 𝑡) = 10 + 2𝑥.

Fig. 11 .FIGURE 11

 1111 FIGURE 11 Comparison of the deflection of the HDPE beam with FE and FDKV models under different load conditions: (a) 𝑓(𝑥, 𝑡) = 10𝐻𝑒𝑎𝑣𝑖𝑠𝑖𝑑𝑒(𝑡) ; (b) 𝑓(𝑥, 𝑡) = 30𝐻𝑒𝑎𝑣𝑖𝑠𝑖𝑑𝑒(𝑡) ; (c) 𝑓(𝑥, 𝑡) = 10𝜋 sin 𝑡; (d) 𝑓(𝑥, 𝑡) = 10 + 2𝑥.

  The numerical solution of the beam deflection is 𝜔 𝑛 (𝑥, 𝑡). If 𝜔 𝑛 (𝑥, 𝑡) is brought into Eq. (38) there must be a residual term called the residual function 𝑅 𝑛 (𝑥, 𝑡):

	then, 𝜔 𝑛 (𝑥, 𝑡) can be expressed as:
	𝐿[𝜔 𝑛 (𝑥, 𝑡)] =	𝜕 2 𝜔 𝑛 (𝑥,𝑡) 𝜕𝑡 2	+ 𝐷 𝐶	𝛽 𝜕 4 𝜔 𝑛 (𝑥,𝑡) 𝜕𝑥 4	+ 𝐷 𝐶	𝛾 𝜕 4 𝜔 𝑛 (𝑥,𝑡) 𝜕𝑥 4	-𝑓(𝑥, 𝑡) = 𝑅 𝑛 (𝑥, 𝑡)	(41)
	The absolute error function is defined as 𝑒 𝑛 (𝑥, 𝑡), which satisfies:
	𝑒 𝑛 (𝑥, 𝑡) = 𝜔(𝑥, 𝑡) -𝜔 𝑛 (𝑥, 𝑡)			(42)
	Based on the definition of 𝐿[ * ]:
	𝐿[𝑒 𝑛 (𝑥, 𝑡)] = 𝐿[𝜔(𝑥, 𝑡)] -𝐿[𝜔 𝑛 (𝑥, 𝑡)] = -𝑅 𝑛 (𝑥, 𝑡)	(43)
	Hence, Eq. (43) is rewritten as:
	𝐿[𝑒 𝑛 (𝑥, 𝑡)] =	𝜕 2 𝑒 𝑛 (𝑥,𝑡) 𝜕𝑡 2	+ 𝐷 𝐶	𝛽 𝜕 4 𝑒 𝑛 (𝑥,𝑡) 𝜕𝑥 4	+ 𝐷 𝐶	𝛾 𝜕 4 𝑒 𝑛 (𝑥,𝑡) 𝜕𝑥 4	-𝑓(𝑥, 𝑡) = -𝑅 𝑛 (𝑥, 𝑡)
	𝐿[𝜔 𝑛 (𝑥, 𝑡)] = 𝑅 𝑛 (𝑥, 𝑡)					(39)
									(40)

So:

𝐿[𝜔(𝑥, 𝑡)] -𝐿[𝜔 𝑛 (𝑥, 𝑡)] = -𝑅 𝑛 (𝑥, 𝑡)

TABLE 1

 1 Absolute and relative errors of the equation with different values of 𝑛.

	t	x	absolute	n 	4	relative	absolute	n 	5	relative	absolute	n 	6	relative
			error			error	error			error	error			error
		0.4	2.4410 -15
	0.2

TABLE 2

 2 The simulation parameters of HDPE and PEEK in FDE constitutive model.

	Material	𝜌 (kg 𝑚 3 ⁄ )	𝛼	𝜂
	HDPE	960	0.1603	3.34110 5
	PEEK	1290	0.2341	5.5010 6

TABLE 3

 3 The simulation parameters of HDPE in FDKV constitutive model.

	Material	𝛽	𝛾	𝜂 1	𝜂 2
	HDPE	0.3320	0.1088	2.87410 5	1.55810 5

TABLE 4 Nomenclature

 4

	Symbol				Explanation
	PEEK				Poly (ether ether ketone)
	HDPE				High density polyethylene
	SCPs				Shifted Chebyshev polynomials
	FDE model		Fractional derivative element model
	FDKV model		Fractional derivative Kelvin-Voigt model
	𝛼, 𝛽, 𝛾, 𝑝 𝜅 , 𝑞 𝜅		Fractional derivative order
	𝜎						Normal stress on the cross section
	𝜀						Strain
	𝜂, 𝜂 1 , 𝜂 2				Material constants in FDE and FDKV models
	𝐶	𝐷 𝛼	, 𝐷 𝐶	𝛽	, 𝐷 𝐶	𝛾	Fractional order differential operator
	Φ 𝑛 (𝑥), Φ 𝑛 (𝑡)		Family of shifted Chebyshev polynomials
	𝐹 𝑖 (𝑥), 𝐺 𝑖 (𝑥), 𝐺 𝑘 (𝑥), 𝐺 𝑗 (𝑥), 𝐺 𝑗 (𝑡)	Shifted Chebyshev polynomials
	𝐴 𝑛 , 𝐴 𝑋 , 𝐴 𝑇 , 𝑈		Coefficient matrix
	𝑃 𝑥 1 , 𝑃 𝑡 1				First order differential operator matrix
	𝑃 𝑡 𝛼					Fractional order differential operator matrix
	𝑍 𝑛 (𝑥)				Family of basical polynomials
	𝜇 𝜅 , 𝜂 𝜅				Material constants in one-dimensional visco-elastic model
	𝑎, 𝑦, 𝑚, 𝐻, 𝑇, 𝜅		Non-negative real number
	𝑐						Real number on the interval [-1, 1]
	𝑏						Width of a rectangular cross section
	ℎ						Height of a rectangular cross section
	𝑥						Position
	𝑡						Time
	𝑓(𝑥, 𝑡)				Load
	𝑙						Length of beam
	𝐴						Cross-sectional area
	𝐼						Moment of inertia
	𝜌						Density of material
	𝑀(𝑥, 𝑡)				Beam bending moment
	𝐿[ * ]					Operator for error correction
	𝑅 𝑛 (𝑥, 𝑡)				Residual function
	𝐸(𝑥, 𝑡)				Correction error function
	𝜔(𝑥, 𝑡)				Exact solution of displacement
	𝜔 𝑛 (𝑥, 𝑡)				Numerical solution of displacement
	𝜔 * (𝑥, 𝑡)				Correction solution of displacement
	𝑅 𝑀					Taylor series remainder
	𝜔 𝑇,𝑀 (𝑥, 𝑡)			𝑀 𝑡ℎ Taylor series expansion of ω(𝑥, 𝑡)
	𝑒 𝑛 (𝑥, 𝑡)				Absolute error function
	𝑒 * (𝑥, 𝑡)				Correction error function
	𝑒(𝑥, 𝑡)				Relative error function

ACKNOWLEDGEMENT

This work is supported by the Natural Science Foundation of Hebei Province (A2017203100) in China and Le Studium research professorship award of Centre-Val de Loire region in France.

Definitions.

Entity: government agency, foundation, commercial sponsor, academic institution, etc.

Grant:

A grant from an entity, generally [but not always] paid to your organization Personal Fees: Monies paid to you for services rendered, generally honoraria, royalties, or fees for consulting , lectures, speakers bureaus, expert testimony, employment, or other affiliations Non-Financial Support: Examples include drugs/equipment supplied by the entity, travel paid by the entity, writing assistance, administrative support, etc.

Other:

Anything not covered under the previous three boxes Pending: The patent has been filed but not issued Issued: The patent has been issued by the agency Licensed: The patent has been licensed to an entity, whether earning royalties or not Royalties: Dr. Chen has nothing to disclose.

Evaluati on and Feedback

Are there other relationships or activities that readers could perceive to have influenced, or that give the appearance of potentially influencing, what you wrote in the submitted work?

Yes, the following relationships/conditions/circumstances are present (explain below):

No other relationships/conditions/circumstances that present a potential conflict of interest At the time of manuscript acceptance, journals will ask authors to confirm and, if necessary, update their disclosure statements. On occasion, journals may ask authors to disclose further information about reported relationships.

Based on the above disclosures, this form will automatically generate a disclosure statement, which will appear in the box below.

Generate Disclosure Statement

Please visit http://www.icmje.org/cgi-bin/feedback to provide feedback on your experience with completing this form.

Section 5

Relationships not covered above . ✔