Jiawei Cao

Dr Yiming Chen
email: chenym@ysu.edu.cn

Yuanhui Wang
email: wangyuanhui527@163.com

Gang Cheng

Thierry Barriere

Prof Abdon Atangana

Cao Jiawei

Keywords: Fractional calculus, Fractional partial differential equation, Dynamic analysis, Polymethyl methacrylate, Viscoelastic model, Shifted Legendre polynomial Fractional calculus, Fractional partial differential equation, Dynamic analysis, Polymethyl methacrylate, Viscoelastic model, Shifted Legendre polynomial

Shifted Legendre polynomials algorithm used for the dynamic analysis of PMMA viscoelastic beam with an improved fractional model

Manuscript ID:CHAOS-D-20-02466

We greatly appreciate the editor's effort in handling our manuscript and the positive and constructive comments from anonymous reviewers.

In the revised version, we have addressed all the comments point by point and concerns to our best understanding and capability. The detailed "Response to the reviews" and "Summary of changes" are presented below.

The valuable comments of the reviewers enable us to improve the quality of our manuscript. We believe our response to reviewers will be satisfactory and the revised version is suitable for publication in the prestigious Chaos, Solitons & Fractals. In the manuscript the comments for Reviewer 1, Reviewer 2 and Reviewer 4 are all highlighted in blue, red and cyan.

Response to the review

Reviewer 1:

Comment 1: In the abstract line 12, "determinated" need to be changed to "determined."

Response: We are very sorry for the writing mistakes. We have revised and marked it in blue in the revised manuscript.

Comment 2: Introduction part: From line 15, there is a repetition of the word "the". From line 37, check the spelling of "mathacrvlate" and also add "an" before amorphous.

From line 55, replace described with describe. From lines 65 and 66, do you mean "time-domain"? From lines 78 and 79, please rephrase the sentence, what do you mean by "preformed"?

Response: We are very sorry for these writing mistakes. We have carefully checked the possible typos/spellings and grammar errors to improve quality of the writing in the whole paper. The improvements are marked in blue in the revised manuscript.

Comment 3: Add hyphen to the word "fractional order" in lines 18, 108, 113, 116 & 152. Response: We have added hyphen for all "fractional order" and marked it in blue in the revised manuscript.

Comment 4: Subscript "T" used in equations 26, 28, 43 and 48 should be defined.

Response: In mathematical expression, the superscript "T" represents the transpose of the matrix, there is no need to define it.

Comment 5: Check the word "discretize" in line 144.

Response:

We have revised and marked in blue in the revised manuscript.

The variable   , x t is discretized by the reasonable match points

    2 1 2 1 , 1, 2, , , 1, 2, , 2 1 2 1
i j i j x H i n t S j n n n          
based on the collection method.

Comment 6: Can you buttress your eq. 50 with more explanation of how you obtained decimal number for partial derivative subscript (i.e. 4.1), subscript t (i.e. 0.1 and 1.9)?

Response: In order to prove the accuracy of the algorithm proposed in this paper, the numerical example is consistent with the governing equation form of the viscoelastic beam, and the parameters are different. The parameters of the numerical example are chosen to be arbitrary values. In order to compare the error between the algebraic solution and the numerical solution, we can derive the expression of   , f x t based on the algebraic solution of the equation, so as to find the numerical solution to compare with the algebraic solution.

Tab. 1 for your table then you can use Tab. 1 not Table 1 for your explanation and so on)

Response: We are very sorry for the difference between the table name and the explanation. We have revised the table name and explanation uniformly and marked it in blue in the revised manuscript.

Comment 10: Check your Figures naming and your explanation of these figures, consistency is very important (i.e. if you are using figure 1 for your diagram then use figure 1, not fig. 1 for your explanation. Also, if you are using fig. 1 for your table then you can use fig. 1 not figure 1 for your explanation and so on)

Response: We are very sorry for the difference between the Figures naming and explanation. We have revised the naming and explanation of the figure and marked it in blue in the revised manuscript.

Comment 11: Check the statement below eq. ( 47) line 3, should be written as "Legendre polynomials. Afterwards, the considered ..."

Response: We are very sorry for the error in segmentation. We have made revised and marked it in blue in the revised manuscript.

Comment 12: Check line 141, "... is determinate ..." This statement needs to be corrected.

Response:

We have changed "is determinate" to "can be obtained".

Comment 13: Check lines 156 & 165, "... determinated ..." This word is not in the dictionary rather "determined".

Response: We are very sorry for these writing mistakes. We have carefully checked the possible typos/spellings and grammar errors to improve quality of the writing in the whole paper. We have changed "determinated" to "determined".

Comment 14: Below table 2, this statement " After to have ..." needs to be reconstructed.

Response:

We have revised this sentence and marked it in blue in the revised manuscript.

Reviewer 2:

Comment 1: Please check the entire manuscript for grammatical errors/typos.

Response:We have carefully checked the possible typos/spellings and grammar errors to improve quality of the writing in the whole paper.

Comment 2: In Eq. 1,

m m       1 , 0 should be m m       1 , 0 .
Response: We are very sorry for the error in Eq. 1. We have modified Eq.

1.

Comment 3: Just for verification, the authors can re-verify if the absolute value of x=1.6,t=4 and n=5 is indeed in Table 1.

Response:

We are very sorry for this error, the absolute error value should be 9 0.0232 10   Comment 4: The author needs to improve the introduction and talk also about analytical methods of solving fractional PDEs. For this purpose, the authors can add the following references to enrich the introductory section.

Response:

We have revised the format of references and added relevant documents suggested by reviewers.

Reviewer 4:

Comment 1: The highlights are effective, however, the last highlight may be reworded.

Response:

The last highlight has been revised to "Dynamic analysis of PMMA beam under different load conditions and temperatures". Response: The stress is uniaxial normal stress and the strain is normal strain. They are not tensors, because we are interested only the stress and strain in the uniaxial direction of the beam. PMMA used in this paper is isotropic viscoelastic material, which exhibits identical material properties in all directions at every given point.

Comment 4: The authors may use the complex elastic modulus instead of the complex modulus in all of the manuscript.

Response:

We have revised all "complex modulus" to "complex elastic modulus" and marked them in cyan.

Comment 5: The proof of the equation-11 and the equation-12 must be illustrated.

Response:

The proof process of the formula has been given and is marked in cyan in the revised manuscript.

Comment 6: In the most literatures, the symbol used for displacement is Y rather than  and the authors are encouraged to use this.

Response:

We have accepted the reviewer's suggestion to replace all   , x t  in the paper and figures with   , y x t .

Comment 7: The authors used the same symbol T for two different things.

In the equations from 19 and after that in Legendre polynomial and before that for kinetic energy. Also, the same thing is repeated with the symbols A which was used for three different aspects. The authors must reword the used symbols and use different symbols for different things.

Response

Response:

We have added references suggested by reviewers.

Comment 11: In the figure-2, the authors may used the traditional symbol for the frequency or if they mean this is the angular frequency, they have to change the unit.

Response: This paper used the traditional frequency and we have revised the symbol of the frequency.

Comment 12: In a lot of experimental works and engineering application the function f (x, t) = A sin(100<pi>t) instead of f (x, t) = A sin(0.5t) when the unit of the time is second. The authors should explain why they used this form. Also, in the same function, there is no need for the unit and the unit should be put after the values of A in the same subsection.

Response: According to the reviewer's comment, we modified the form of the simple harmonic load in the paper, and made corresponding changes to the relevant pictures.

Comment 13: The authors may add more details about expanding the method to other fluids and how it can be applied.

Response: In the conclusion, we added a paragraph to illustrate the significance and importance of the algorithm proposed in this paper.

Not only the unsuited contents have been revised one by one according to the reviewer's comments, but also the careful polishing has been done to the whole manuscript. Any changes haven been marked in the revised manuscript.

Thanks again for your help for our manuscript.

Jiawei Cao

Summary of changes 1. Page 1, in the abstract line 12, "determinated" has been revised to "determined".

2. Page 1, 4, 21, all "time domain" have been revised to "time-domain".

3. Page 2, 6, 13, 14, 22, all "fractional order" have been revised to "fractional-order". 7. Page 3, line 45, "mathacrvlate" has been revised to "methacrylate" and "an" has been added before "amorphous".

8. Page 3, line 63, "described" has been revised to "describe". 9. Page 4, "Akinyemi and Lyiola [START_REF] Akinyemi | Exact and approximate solutions of timefractional models arising from physics via Shehu transform[END_REF] used the Shehu transform method to decompose the nonlinear term to solve the time fractional differential equations in physics and engineering. Senol et al. [START_REF] Iyiola | Efficient analytical techniques for solving time-fractional nonlinear coupled Jaulent-Miodek system with energy-dependent Schrödinger potential[END_REF] proposed an analytical approximate method based on residual power series and q-homotopy analysis method [START_REF] Akinyemi | q-Homotopy analysis method for solving the seventh-order time-fractional Lax's Korteweg-de Vries and Sawada-Kotera equations[END_REF] to solve the nonlinear time fractional equations. Akinyemi et al. [START_REF] Akinyemi | Iterative methods for solving fourth and sixth order time-fractional Cahn-Hillard equation[END_REF] developed a new iterative method with q-homotopy analysis to solve the fourth and sixth order time fractional

Cahn-Hilliard equations. The simplicity and accuracy of the methods were confirmed by the provided examples. Al-Raeei and El-Daher [START_REF] Al-Raeei | Numerical simulation of the space dependent fractional Schrdinger equation for London dispersion potential type[END_REF] developed numerical-integral methods to find numerical solutions of the fractional Schrodinger equation without using Laplace transform." has been added. [START_REF] Martin | Stability approach to the fractional variational iteration method used for the dynamic analysis of viscoelastic beams[END_REF]. Page 4, line 96, "In this paper, the dynamic analysis of the viscoelastic beam is preformed under various loading conditions" has been revised to "In this paper, the dynamic analysis of viscoelastic beam under various loading conditions is studied".

11. Page 4, line 99, "models are" has revised to "model is".

12. Page 5, "

0, 1 m m       " has revised to " 0, 1 m m       ".
13. Page 6, Section 2.2, "constitution" has revised to "constitutive". 15. Page 7, "complex modulus" has revised to "complex elastic modulus".

16. Page 7, Eq. 9 has been rewritten.

Page 7, "

A is the cross-sectional area" has revised to " x A is the cross-sectional area".

    2 1 2 1 , 1, 2, , , 1, 2, , 2 1 2 1
i j i j x H i n t S j n n n           have been used to discrete the variable   , x t to   , i j x t ." has revised to "The variable   , x t
is discretized by the reasonable match points

    2 1 2 1 , 1, 2, , , 1, 2, , 2 1 2 1
i j i j x H i n t S j n n n          
based on the collection method."

22. Page 14, line 177, "is determinate" has revised to "can be obtained".

23. Page 14, line 178, "using Matlab platform and least square method" has been revised to "using the numerical algorithm in MATLAB platform.

The fractional-order equation is transformed to algebraic equations. The least square method is used to calculate the solution of the algebraic equations.".

24. Page 14, line 186, " in which the left term is similar to the governing equation of the viscoelastic beam (defined in Eq. ( 17)) and   , f x t is derived from algebraic solution." has been added. [START_REF] Chen | Numerical study of a class of variable order nonlinear fractional differential equation in terms of Bernstein polynomials[END_REF]. Page 15, line 189, "determinated" has revised to "determined".

26. Page 15, Table 1, " ".

32. Page 15, "After to have implemented these identified parameters in Eq. ( 11) and ( 12), the values of the simulative modulus E and E   are in excellent agreements with the experimental modulus, see Figure 2" has deleted.

33. Page 16, Table . 2, "isotropic" has been added.

34. Page 17, line 216, "is" has been revised to "are".

35. Page 17, line 219, "to" has been revised to "on".

36. Page 17, before Eq. ( 53), "are" has been revised to "are as".

Highlights (for review)

Shifted Legendre polynomials algorithm used for the dynamic analysis of PMMA viscoelastic beam with an improved fractional model

Introduction

Fractional calculus has been widely used in the fields of science and engineering in recent years. The fractional-order models of the materials are integrated into the integrals and differential equations to improve their behaviour modelling efficiency and accuracy. Fractional calculus was used to investigate the viscoelastic behaviour of the fluid in the pipes [START_REF] Wang | Analytical and numerical study of electroosmotic slip flows of fractional second grade fluids[END_REF] and the dynamic behaviour of the active elastomers in magnetic field [START_REF] Nadzharyan | Fractional rheological models of dynamic mechanical behavior of magnetoactive elastomers in magnetic fields[END_REF]. Fractionalorder Constitutive Equation (FCE) was proposed to describe the physical behaviour of non-Newtonian fluid [START_REF] Sun | A space fractional constitutive equation model for non-Newtonian fluid flow[END_REF]. The results showed that the proposed equations could successfully capture the observed increasing of shear stress for different velocity gradients. Meng et al. [START_REF] Meng | Parameter study of variable order fractional model for the strain hardening behavior of glassy polymers[END_REF] proposed a new variableorder fractional derivative viscoelastic model to describe the strain hardening behaviour of amorphous glassy polymers, in which the order function was assumed to be linearly varied with time. The comparison between the model prediction and experimental data confirmed the close relationships of order change and strain hardening. A variable-order Fractional Differential Equation (FDE) was developed to modelling the shape memory behaviour of the polymers [START_REF] Li | A variable-order fractional differential equation model of shape memory polymers[END_REF], which proves more suitable than the conventional constant order FDE. A fractional model composed parallel fractional Maxwell elements was used to describe the mechanical behaviour of the polymers [START_REF] Lei | A fractional model with parallel fractional Maxwell elements for amorphous thermoplastics[END_REF]. It was able to describe the evolution of the master curves of storage modulus and loss modulus during the stress relaxation experiments. A viscoelastic model employing fractional-order derivatives was applied to describe the dielectric properties of materials [START_REF] Wharmby | Modifying Maxwell's equations for dielectric materials based on techniques from viscoelasticity and concepts from fractional calculus[END_REF]. Henriques et al. [START_REF] Henriques | Viscoelastic behavior of polymeric foams: Experiments and modeling[END_REF] used the fractional derivative model to describe the viscoelastic behaviour of two polymeric foams. Excellent correlation between the experimental data and model predictions was observed for the shear storage modulus. Fractional calculus has shown great advantages in modelling the physical constitutive behaviours of materials, especially the viscoelastic properties of polymeric materials.

The viscoelastic beam is an important structural element to resist load in engineering. Many fractional material models have been proposed to analyze the dynamic response of the beam under various loading conditions. A fractional Kelvin-Voigt model was used to describe the viscoelastic properties of the beam. The numerical results revealed that the increase of the derivative order could cause a decrease in the vibration amplitudes of the beam [START_REF] Freundlich | Transient vibrations of a fractional Kelvin-Voigt viscoelastic cantilever beam with a tip mass and subjected to a base excitation[END_REF]. The dynamic analysis of a simple supported viscoelastic beam with a fractional Zener model was performed by using a modified variational iteration method [START_REF] Martin | Stability approach to the fractional variational iteration method used for the dynamic analysis of viscoelastic beams[END_REF]. The numerical examples were effectuated to study the influences of fractional derivative on the dynamic response of the structure. The comparison with the classical Zener model proved the efficiency of the proposed algorithm for solving the fractional governing equations. It is considered as an important issue to propose the efficient fractional constitutive viscoelastic models to describe accurately the physical constitutive behaviour of the beams in mechanical engineering.

Polymethyl methacrylate (PMMA) is an amorphous thermoplastic polymer, which is widely used in the manufacturing of medical devices due to its excellent mechanical and optical properties. The mechanical property of PMMA material above its glass transition temperature is greatly affected for further applications. Its mechanical properties above the glass transition temperature are greatly affected for further applications in the life prediction and design of engineering structures. Jo et al. [START_REF] Jo | Constitutive modeling for mechanical behaviour of PMMA microcellular foams[END_REF] described the nonlinear tensile behaviours of PMMA foam by different integer viscoelastic models.

The proposed constitutive equations were developed in terms of strain, strain rate, elastic modulus, relative density of foam and relaxation time constant. The nonlinear tensile stress-strain behaviours were well described. Varghese

and Batra [START_REF] Varghese | Constitutive equations for thrmomechanical deformations of glassy polymers[END_REF] utilized the modified viscoplastic constitutive equations to simulate the mechanical behaviour of PMMA at high strain rates. The proposed models were found to be well coherent with the experimental results available in the literature. Cheng et al. [START_REF] Cheng | Physical modelling and identification of polymer viscoelastic behaviour above glass transition temperature and application to the numerical simulation of the hot embossing process[END_REF] used a generalized Maxwell model to describe the viscoelastic behaviour of PMMA lightly about its glass transition temperature. The generalized Maxwell model parameters were determined via dynamic mechanical analysis. The proposed model was applied to describe the relaxation modulus to achieving the numerical simulation of the hot embossing process [START_REF] Cheng | Physical modelling, numerical simulation and experimental investigation of microfluidic devices with amorphous thermoplastic polymers using a hot embossing process[END_REF]. However, there is little research work concerning the fractional constitutive viscoelastic model for PMMA above its glass transition temperature.

A suitable numerical algorithm has to be proposed to solve the fractional governing equation of the beam established based on the viscoelastic material models. The ordinary numerical algorithms are developed by using the multi-scale method [START_REF] Chang | Vibration and stability of an axially moving Rayleigh beam[END_REF], Galerkin method [START_REF] Permoon | Application of radial basis functions and sinc method for solving the forced vibration of fractional viscoelastic beam[END_REF], finite element method [START_REF] Demir | Linear dynamical analysis of fractionally damped beams and rods[END_REF] and so on. One of the difficulties with these methods is that the timedomain solutions could not be obtained directly. The Laplace transform has to be used to transform time-domain fractional equations into frequency domain fractional equations. The time-domain solution is then obtained by applying the inverse Laplace transform to the frequency domain solutions.

The complexity of Laplace and inverse transform makes it difficult to solve effectively the fractional equations. Akinyemi and Lyiola [START_REF] Akinyemi | Exact and approximate solutions of timefractional models arising from physics via Shehu transform[END_REF] used the Shehu transform method to decompose the nonlinear term to solve the time fractional differential equations in physics and engineering. S ¸enol et al. [START_REF] Iyiola | Efficient analytical techniques for solving time-fractional nonlinear coupled Jaulent-Miodek system with energy-dependent Schrödinger potential[END_REF] proposed an analytical approximate method based on residual power series and q-homotopy analysis method [START_REF] Akinyemi | q-Homotopy analysis method for solving the seventh-order time-fractional Lax's Korteweg-de Vries and Sawada-Kotera equations[END_REF] to solve the nonlinear time fractional equations. Akinyemi et al. [START_REF] Akinyemi | Iterative methods for solving fourth and sixth order time-fractional Cahn-Hillard equation[END_REF] developed a new iterative method with q-homotopy analysis to solve the fourth and sixth order time fractional Cahn-Hilliard equations. The simplicity and accuracy of the methods were confirmed by the provided examples. Al-Raeei and El-Daher [START_REF] Al-Raeei | Numerical simulation of the space dependent fractional Schrdinger equation for London dispersion potential type[END_REF] developed numerical-integral methods to find numerical solutions of the fractional Schrödinger equation without using Laplace transform. This is the reason why the algorithms based on the polynomials are proposed to solve fractional differential equations. The polynomial algorithms can solve directly the fractional differential equations in the time-domain. They exhibit many other advantages such as high accuracy and low complexity [START_REF] Wang | Numerical solutions for systems of fractional order differential equations with Bernoulli wavelets[END_REF]. The commonly used polynomial algorithms include the Legendre polynomial [START_REF] Meng | Numerical solutions of nonlinear fractional differential equations by alternative Legendre polynomials[END_REF], the Bernstein polynomial [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 the Chebyshev polynomial [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]. Therefore, the polynomial algorithm is proposed for solving the fractional governing equation for PMMA viscoelastic beams.

In this paper, the dynamic analysis of viscoelastic beam under various loading conditions is studied. An improved fractional viscoelastic model is proposed to describe the constitutive behaviour of PMMA above its glass transition temperature. The proposed model is characterized in different temperatures via dynamic mechanical analysis. The fractional governing equations of the beam are established by integrating the fractional viscoelastic model. A numerical algorithm based on the polynomial and wavelet method is proposed to solve directly in the time-domain. The shifted Legendre polynomials retained can resolve effectively and accurately the complex fractional differential equations over a larger interval. The displacement of the viscoelastic beam versus time and position are obtained.

The fractional derivative, the fractional viscoelastic constitutive model and the governing equation of the viscoelastic beam are introduced in Section 2. The shifted Legendre polynomials are presented in Section 3. The numerical example is proposed to verify the accuracy of the algorithm in Section 4. The dynamic analysis of the PMMA beam is performed under various loading conditions and temperatures in Section 5. The research work is concluded in Section 6. A list of the symbols used in the paper is given in Appendix A.

Preliminary knowledge and definitions

Caputo fractional derivative operator Definition 2.1 The Caputo fractional derivative operator D α

x of order α is defined as [START_REF] Chen | Numerical solution of fractional partial differential equations with variable coefficients using generalized fractionalorder Legendre functions[END_REF][START_REF] Yi | Wavelet operational matrix method for solvig fractional differential equations with variable coefficients[END_REF]

D α x u(x) =          1 Γ(m -α) x 0 (x -τ ) -α+m-1 u (m) (τ ) dτ, α > 0, m -1 < α < m d (m) u(x) dx m , α = m (1)
where x ≥ 0, and m ∈ N (N denotes positive integer).

The gamma function, denoted by Γ (•) is defined as Γ (z) = ∞ 0 e -t t z-1 dt for complex arguments with positive real part.

Several significant properties of the Caputo fractional derivative are shown as follows:

D α x C = 0 ( 2
)
where C is constant.

The Caputo fractional derivative is a linear operator.

D α x (λ(u(x))) = λD α x u(x) (3)
where λ is constant. The Caputo fractional derivative of the power function satisfies

D α x x m =    0, m = 0 Γ(m + 1) Γ(m + 1 -α) x m-α , m = 1, 2, 3, • • • (4)

Fractional constitutive viscoelastic model

The general form of conventional viscoelastic stress strain relationship is expanded in Taylor series as follows [START_REF] Wilhelm | Viscoelasticity[END_REF]:

n r=0 p r d r σ(t) dt r = n r=0 q r d r ε(t) dt r (5)
where σ(t) is stress, ε(t) is strain, r ∈ N (N denotes positive integer).

The traditional differential equations established by integer order operators are limited to describe the viscoelastic behaviour of the material in mechanical engineering. The fractional-order differential operators are widely used in fractional viscoelastic constitutive models thanks to their memory-dependent properties. Several fractional models such as Maxwell model [START_REF] Shen | Unsteady MHD flow and heat transfer of fractional Maxwell viscoelastic nanofluid with Cattaneo heat flux and different particle shapes[END_REF], Kelvin-Voigt [START_REF] Xu | Creep constitutive models for viscoelastic materials based on fractional derivatives[END_REF] and Zenner [START_REF] Mokhtari | Aeroelastic analysis of sandwich cylinder with fractional viscoelastic core described by Zener model[END_REF] models are applied to describe the viscoelastic properties. The current integer order viscoelastic model (Eq. ( 5

where α r , β r ∈ R (R denotes real number). For complex materials, choosing one or two items does not describe the properties of the material well. Therefore, one item is selected on the left side of Eq. ( 6), and three items are selected on the right side to represent the constitutive relationship of the viscoelastic material, which includes integer order and fractional-order derivatives. The more terms are taken, the closer the constitutive equation is to the actual situation. The following formula is selected as the constitutive equation in this paper.

σ(t) = q α D α t ε(t) + q 0 ε(t) + q 1 D 2 t ε(t) (7)
where q α , q 0 , q 1 ∈ R, α ∈ (0, 1).

Taking Fourier transform of both sides of Eq. ( 7) yields σ(w) = q α (iw) α ε(w) + q 0 ε(w) + q 1 (iw) 2 ε(w) [START_REF] Henriques | Viscoelastic behavior of polymeric foams: Experiments and modeling[END_REF] where w is frequency.

The complex elastic modulus E * is expressed as:

E * = σ(w) ε(w) = q α (iw) α + q 0 + q 1 (iw) 2 = q α w α cos( π 2 α) + iq α w α sin( π 2 α) + q 0 -q 1 w 2 (9)
E * = E + iE ( 10
)
where E is storage modulus, corresponding the real part of E * , E is loss modulus, corresponding the imaginary part of E * . So, E and E are expressed as:

E = q α w α cos( π 2 α) + q 0 -q 1 w 2 (11)
E = q α w α sin( π 2 α) (12)

Fractional governing equations of viscoelastic beam

The deformation diagram of the viscoelastic beam is shown in Figure 1. The differential equation of motion of the viscoelastic beam is established according to Hamilton's principle [START_REF] Jinkyu | Extended framework of Hamilton's 415 principle for continuum dynamics[END_REF]. The kinetic energy [START_REF] Wang | Shifted-Chebyshev-polynomial-based numerical algorithm for fractional order polymer visco-elastic rotating beam[END_REF] of the beam is expressed as

E k = 1 2 H 0 ρA x ∂ 2 y (x, t) ∂t 2 dx ( 13
)
where ρ is the material density, A x is the cross-sectional area, H is the length 140 of the beam and y(x, t) is the displacement. The potential energy [START_REF] Wang | Shifted-Chebyshev-polynomial-based numerical algorithm for fractional order polymer visco-elastic rotating beam[END_REF] is expressed as

V = 1 2 H 0 q 0 I ∂ 2 y (x, t) ∂x 2 + q α I ∂ 2+α y (x, t) ∂x 2 ∂t α + q 1 I ∂ 4 y (x, t) ∂x 2 ∂t 2 dx ( 14
)
where I is the moment of inertia, q α , q 0 , q 1 , α are the identified parameters of the constitutive viscoelastic model. The work done by the external forces is expressed as

W = 1 2 H 0 f (x, t)y (x, t) dx (15)
where f (x, t) is the external transverse axial load applied on the beam.

According to the Hamilton's principle

δ t 2 t 1 (E k -V )dt + t 2 t 1 δW dt = 0 ( 16
)
The fractional governing equation of motion of the viscoelastic beam is obtained:

I(q α ∂ 4+α y(x, t) ∂x 4 ∂t α + q 0 ∂ 4 y(x, t) ∂x 4 + q 1 ∂ 6 y(x, t) ∂x 4 t 2 ) + ρA x ∂ 2 y(x, t) ∂t 2 = f (x, t) (17)
The beam is fixed at both ends and the boundary conditions are:

y(x, 0) = 0 y(0, t) = 0 y(H, t) = 0 y x (0, t) = 0 y x (H, t) = 0 (18) The shifted Legendre polynomial of degree n in [0, 1] is defined as [37]

Numerical algorithm of fractional governing equations

l n,i (x) = n-i k=0 (-1) k n -i k n + i + k + 1 n -i x i+k ( 19
)
where i = 0, 1, • • • , n, x ∈ [0, 1]. Then ϕ n (x) is formulated as ϕ n (x) = [l n,0 (x), l n,1 (x), • • • , l n,n (x)] T ϕ n (x) = AT n (x) ( 20
)
where

T n (x) = [1, x, • • • , x n ] T , A = [a ij ] n i,j=0 , a ij =      0, 0 ≤ k ≤ i (-1) k n -i k n + i + k + 1 n -i , k ≤ i ≤ n
The shifted Legendre polynomial of degree n in [0, H] is formulated as

L n,i (x) = n-i k=0 (-1) k n -i k n + i + k + 1 n -i x H i+k = n-i k=0 (-1) k n -i k n + i + k + 1 n -i 1 H i+k x i+k , ( 21
)
where i = 0, 1, • • • , n, x ∈ [0, H].
Then Φ n (x) is formulated as

Φ n (x) = ALT n (x) (22)
where

L = [l ij ] n i,j=0 , l ij = 0, i = j H -i , i = j
The shifted Legendre polynomial of degree n in [0, S] is formulated as

Ln,i (t) = n-i k=0 (-1) k n -i k n + i + k + 1 n -i t S i+k = n-i k=0 (-1) k n -i k n + i + k + 1 n -i 1 S i+k t i+k , (23)
where Then φ n (t) is defined as

i = 0, 1, • • • , n, t ∈ [0, S].
φ n (t) = AM T n (t) (24)
where

M = [m ij ] n i,j=0 , m ij = 0, i = j S -i , i = j

Function approximation

A continuous function y(x) in the domain [0, H] can be expanded in terms of shifted Legendre polynomials as y(x) = lim n→∞ n i=0 c i L n,i (x), y(x) can be approximated as

y(x) ≈ y n (x) = n i=0 c i L n,i (x) = C T Φ n (x) ( 25
)
where n is the number of terms of the shifted Legendre polynomial,

C T = [c 0 , c 1 , • • • , c n ]. Then C T Φ n (x), Φ T n (x) = y(x), Φ T n (x) C T = y(x), Φ T n (x) Q -1 Q = Φ n (x), Φ T n (x) = [δ ij ] n i,j=0 (26)
where

δ ij = H 0 L n,i (x)L n,j (x)dx =    0, i = j H i + j + 1 , i = j (i, j = 0, 1, • • • , n). 160
Similarly, a continuous function y(t) in the domain [0, S] can be expanded in terms of shifted Legendre polynomials as y(t) = lim n→∞ n i=0 k i Ln,i (t), y(t) can be approximated as

y(t) ≈ y n (t) = n i=0 k i Ln,i (t) = K T φ n (t) (27)
where

K T = [k 0 , k 1 , • • • , k n ]. Then K T Φ n (t), Φ T n (t) = y(t), Φ T n (t) K T = y(t), Φ T n (t) P -1 P = Φ n (t), Φ T n (t) = [∆ ij ] n i,j=0 (28)
where

∆ ij = S 0 Ln,i (t) Ln,j (t)dt =    0, i = j S i + j + 1 , i = j (i, j = 0, 1, • • • , n).
Two-variable continuous function y(x, t) ∈ L 2 [0, H] × [0, S] can be formulated as

y(x, t) = lim n→∞ n j=0 ( n i=0 c i L n,i (x))k j Ln,j (t) = lim n→∞ n j=0 n i=0 c i k j L n,i (x) Ln,j (t) = lim n→∞ n j=0 n i=0 ω ij L n,i (x) Ln,j (t) = lim n→∞ (Φ n (x)U φ n (t)) (29)
After intercepting the finite term, it can get

y(x, t) ≈ Φ n (x)U φ n (t) (30)
where U = [y ij ] n i,j=0 is a matrix of displacement functions.

Integer order operator matrix

The derivative of Φ n (x) with respect to x is formulated as

dΦ n (x) dx = D x Φ n (x) (31)
Then

D x Φ n (x) = D x ALT n (x) = AL dT n (x) dx = ALET n (x) = ALE(AL) -1 Φ n (x) (32)
where

E = [e ij ] n i,j=0 , e ij = 0, i = j + 1 i, i = j + 1
From Eq. ( 32), D x is obtained as

D x = ALE(AL) -1 (33)
The m exponent th derivative of Φ n (x) with respect to x is formulated as

d m Φ n (x) dx m = D mx Φ n (x), m ∈ N (34)
Then

D mx Φ n (x) = D mx ALT n (x) = AL d m T n (x) dx m = ALE m T n (x) = ALE m (AL) -1 Φ n (x) (35)
Based on Eq. ( 35), D mx can be obtained as

D mx = ALE m (AL) -1 (36)
Similarly, φ n (t) is a series of polynomial matrices with respect to t, the derivative of φ n (t) with respect to t is formulated as

dφ n (t) dt = D t φ n (t) (37)
Then

D t φ n (t) = D t AM T n (t) = AM dT n (t) dt = AM ET n (t) = AM E(AM ) -1 φ n (t) (38)
Therefore, from Eq. ( 38), D t is obtained as

D t = AM E(AM ) -1 (39)
The v exponent th derivative of φ n (t) with respect to t is formulated as

d v φ n (t) dt v = D vt φ n (t), v ∈ N (40)
Then

D vt φ n (t) = D vt AM T n (t) = AM d v T n (t) dt v = AM E v T n (t) = AM E v (AM ) -1 φ n (t) (41
) From Eq. (41), D vt is obtained as

D vt = AM E v (AM ) -1 (42)
Thus,

∂ m+v y(x, t) ∂x m ∂t v ≈ ((D mx Φ n (x)) T U (D vt φ n (t))) (43)
where only U is unknown.

Fractional-order operator matrix 170

The α th derivative of φ n (t) can be expressed as:

d α φ n (t) dt α = D αt φ n (t), α ∈ (0, 1) (44)
D αt φ n (t) = D αt AM T n (t) = AM d α T n (t) dt α = AM F T n (t) = AM F (AM ) -1 φ n (t) (45)
where

F = [f ij ] n i,j=0 , f ij =    Γ(i + 1) Γ(i + 1 -α) t -α , i = j, i ≥ 1 0, otherwise
In view of Eq. ( 45), it gives that

D αt = AM F (AM ) -1 (46)
Thus,

∂ m+α y(x, t) ∂x m ∂t α ≈ ((D mx Φ n (x)) T U (D αt φ n (t))) (47)
where only U is unknown. Basic of the above work, the unknown function is replaced with the approximate function containing differential operator matrix of shifted Legendre polynomials. Afterwards, the considered fractional governing Eq. ( 17) will be converted into the following algebraic equation with unknown U .

I(q α ((D 4x Φ n (x)) T U (D αt φ n (t))) + q 0 ((D 4x Φ n (x)) T U (φ n (t))) +q 1 ((D 4x Φ n (x)) T U (D 2t φ n (t))) + ρA x ((Φ n (x)) T U (D αt φ n (t))) = f (x, t) (48)
The boundary conditions Eq. ( 18) can be transformed into the following form:

Φ n (x)U φ n (0) = 0, Φ n (0)U φ n (t) = 0, Φ n (H)U φ n (t) = 0 D x Φ n (0)U φ n (t) = 0, D x Φ n (H)U φ n (t) = 0 (49)
The variable (x, t) is discretized by the reasonable match points

x i = 2i-1 2(n+1) H, i = 1, 2, • • • , n, t j = 2j-1 2(n+1) S, j = 1, 2, • • • , n
, based on the collection 175 method. Eq. ( 48) is transformed into a set of algebraic equations.

The coefficient y ij (i = 0, 1, 2, • • • n; j = 0, 1, 2, • • • n) can be obtained by using the numerical algorithm in MATLAB platform. The fractional-order equation is transformed to algebraic equations. The least square method is used to calculate the solution of the algebraic equations. The numerical 180 solution of the fractional derivative equation can be obtained.

The proposed algorithm can be summarized as follows:

Algorithm: Numerical solution for the fractional governing equation Input:

α, q α , q 0 , q 1 , ρ, A x , I, f (x, t) Output:

1. Function approximation y(x, t) ≈ Φ n (x)U φ n (t) 2.
Derive calculation of integer and fractional operator matrices 3.

Substitute the operator matrix into the initial equation 4.

The initial equation is transformed into an algebraic equation 5.

Let

x i = 2i-1 2(n+1) H, i = 1, 2, • • • , n, t j = 2j-1 2(n+1) S, j = 1, 2, • • • , n 6.
Solve algebraic equations with MATLAB mathematical software 7.

Obtain the solution of the initial equation y (x, t)

Numerical example

In this section, numerical example is given to demonstrate the effectiveness and accuracy of the algorithm in solving fractional-order equations.

185

The following equation is considered as a numerical example, in which the left term is similar to the governing equation of the viscoelastic beam (defined in Eq. ( 17)) and f (x, t) is derived from algebraic solution.

2x 3 ( ∂ 4.1 y(x, t) ∂x 4 ∂t 0.1 + ∂ 4 y(x, t) ∂x 4 + ∂ 6 y(x, t) ∂x 4 ∂t 2 ) + x 2 ∂ 2 y(x, t) ∂t 2 = f (x, t) f (x, t) = 48 Γ(3) Γ(2.9) x 3 t 1.9 + 48x 3 t 2 + 96x 3 + 2x 4 (2 -x) 2 (50)
where the defined boundary conditions as follows:

y(x, 0) = 0 y(0, t) = 0 y(2, t) = 0 y x (0, t) = 0 y x (2, t) = 0

The algebraic solution of Eq. ( 50) is giving by the following equation:

y(x, t) = x 2 (2 -x) 2 t 2 , x ∈ [0, 2], t ∈ [0, 10] (51)
The numerical solutions of Eq. ( 50) are obtained by the proposed method. The solutions are calculated with different n to validate the efficiency of the method. The values of absolute error between the algebraic and approximate solutions are determined at different x and t, which is summarized in Table 1. The absolute errors are approximately 10 -8 . The absolute error decrease with n. In this research, n is selected as 6 to reduce the absolute error in the calculation.

Eq. ( 50) is in the same fractional form of the fractional governing equation of viscoelastic beam. The right term of Eq. ( 50) is more general and complex than uniformly distributed load. 5. Results and discussions

Determination of material parameters

The values of E and E at different temperatures and frequency were determined by using the dynamic mechanical analysis by Cheng et al. [START_REF] Cheng | Physical modelling, numerical simulation and experimental investigation of microfluidic devices with amorphous thermoplastic polymers using a hot embossing process[END_REF]. The cylindrical compression samples were tested by using the harmonic solicitations with different frequency at various temperatures. The parameters in the viscoelastic constitutive equation are identified based on E and E at T =110 • C, 115 • C, 120 • C and 125 • C.

The identified parameters q α , q 0 , q 1 , α were obtained by inverse method with a least square regression according to Eq. ( 11) and [START_REF] Varghese | Constitutive equations for thrmomechanical deformations of glassy polymers[END_REF]. The identified parameters are summarized in Table 2.

Dynamic analysis of viscoelastic beam under various loading conditions

The dynamic displacement of the PMMA beam of rectangular crosssection is calculated in this section. The input geometric and material parameters of the beam are summarized in Table 3. The viscoelastic parameters are defined in Table 2. Based on the geometric and material parameters of the beam, Eq. ( 17) is developed as follows:

25 6 (q α ∂ 4+α y(x, t) ∂x 4 ∂t α + q 0 ∂ 4 y(x, t) ∂x 4 + q 1 ∂ 6 y(x, t) ∂x 4 ∂t 2 ) +23 ∂ 2 y(x, t) ∂t 2 = f (x, t) (52)
The boundary and initial conditions are as follows:

y(x, 0) = 0 y(0, t) = 0 y(1.6, t) = 0 y x (0, t) = 0

y x (1.6, t) = 0 , x ∈ [0, 1.6], t ∈ [0, 60] (53)

Uniformly distributed load

The uniform transverse distributed load 100 N/m is imposed on the beam. The numerical solution of the displacement of the beam is obtained by using the proposed algorithm.

The displacement of the beam at different t and x are shown in Figure 3. The displacement curves at t = 6, 12 and 60 s, exhibit the maximum value in the middle of the beam and zero value at both ends of the beam.

Further analysis shows that the loading time of uniformly distributed load has an effect on the dynamic displacement of the PMMA beam. As time increases, the dynamic displacement of the PMMA beam gradually increases. In engineering, this work provides a theoretical basis for the protection of load-bearing structures.

Harmonic load

The harmonic load f (x, t) = B sin(100πt) is applied on the viscoelastic beam. The transverse displacement of the viscoelastic beam is obtained with B = 100, 200, 500, and 1000 N/m, as shown in Figure 4 These results show that the temperature rising leads to the internal structure changes of the material, also the PMMA beam is weakened against the external deformation and the dynamic response becomes more obvious with time varying as the temperature rising.

Conclusions

In this paper, the shifted Legendre polynomials is used to solve the governing equation of PMMA viscoelastic beam directly in the time-domain. An improved fractional model is proposed to simulate the constitutive relationship of PMMA above the glass transition temperature. The fractional derivative model is used to describe the viscoelastic behaviour of the beam to predict the behaviour of the beam under uniform and simple harmonic loads. The numerical solution of PMMA viscoelastic beam displacement is obtained. At the same time, the displacements of PMMA beams at different temperatures are analyzed.

1. It is verified that the presented algorithm is accurate and effective through numerical example, and is also found that errors are decreasing with increasing n. Further, the shifted Legendre polynomials algorithm is excepted to solve fractional differential equations with three and more variables. More importantly, the proposed algorithm has shown great potentials for highprecision engineering computing problems.

2. The improved fractional-order model is used to derive the fractionalorder governing differential equations of viscoelastic beams, and the dynamic characteristics of viscoelastic materials are studied. It fills the blank in the field of high precision constitutive model simulation of PMMA materials at glass transition temperature.

3. Analysis shows that the time of uniform load affects the dynamic displacement of PMMA beams. With the time increasing, the dynamic displacement is increasing cumulatively in a fluctuating manner. In engineering, this work provides a theoretical basis for the protection of load-bearing structures above glass transition temperature.

4. The temperature rising leads to the internal structure changes of the material, also the PMMA beam is weakened against with the external deformation and the dynamic response becomes more obvious with time varying as the temperature rising. The results show that as the temperature increases and the time prolongs, the load resistance of PMMA beam gradually decreases.

The numerical algorithm proposed in this paper could be used to solve other engineering problems. The identification method of the viscoelastic properties is suitable for not only solid materials, but also fluid materials. Dynamic mechanical analysis is always used to characterize the viscoelastic properties of material. The parameters in the fractional viscoelastic behaviour law can be identified according to the experimental investigations. The governing equation of the problem, associated with the proper parameters of material behaviour, could be efficiently and accurately solved with the proposed numerical algorithm.

[37] Meng ZJ, Yi traditional constitutive model parameter p αr , q αr traditional fractional constitutive model parameter q α , q 0 , q

Comment 2 :Comment 3 :

 23 In the first equation, the authors have defined one type of the Caputo fractional derivatives, however, there is another type. The authors either should add the full definition or should explain any type. Response: In this paper, only the Caputo fractional derivatives are applied in the calculation. This is the reason why only the Caputo fractional operator is defined in Eq (1). It is not necessary to define and explain the other types of fractional derivative. In the equation-5, equation-6 and after that, the authors should explain what are the type of strain and stress. Are they tensors or less and if they are not tensors what about applying for the polymeric materials or other complex materials.

4 . 5 . 6 .

 456 Page 2, line 11, "Meng et al.[START_REF] Meng | Parameter study of variable order fractional model for the strain hardening behavior of glassy polymers[END_REF] proposed a new variable-order fractional derivative viscoelastic model to describe the strain hardening behavior of amorphous glassy polymers, in which the order function was assumed to be linearly varied with time. The comparison between the model prediction and experimental data confirmed the close relationships of order change and strain hardening" has been added. Page 2, "Henriques et al.[START_REF] Henriques | Viscoelastic behavior of polymeric foams: Experiments and modeling[END_REF] used the fractional derivative model to describe the viscoelastic behaviour of two polymeric foams. Excellent correlation between the experimental data and model predictions was observed for the shear storage modulus." has been added. Page 2, line 20, delete a word "the".

14 .

 14 Page 7, 16, 17, 18, 20, all "Fig." have been revised to "Figure".

18 . 20 .

 1820 Page 7, kinetic energy is redefined as T  19. Page 13, "The  exponent t h derivative of   n x  with respect to t is formulated as" has revised to "The th  derivative of   Page 13, "shifted Legendre polynomials, afterwards the considered" has revised to "shifted Legendre polynomials. Afterwards, the considered" 21. Page 14, line 174, "Based on the collection method, the reasonable match points

27 .

 27 Page 15, line 191, "In this research, n is selected as 6 to reduce the absolute error in the calculation." has been added.[START_REF] Chen | Numerical solution of nonlinear fractional integral differential equations by using the second kind Chebyshev wavelets[END_REF]. Page 15, line 199, "determinated" has revised to "determined" 29. Page 15, line 199, "using the dynamic mechanical analysis by" has been added. 30. Page 15, line 200, "The cylindrical compression samples were tested by using the harmonic solicitations with different frequency at various temperature." has been added. 31. Page 15, 17, 20, all "Tab." have been revised to "Table



  37. Page 17, line 225, "the time of uniform load loading" has revised to "the loading time of uniformly distributed load". 38. Page 18, line 231, " 18, line 233, "A = 100, 200, 500 and 1000" has revised to "B = 100, 200, 500 and 1000 N/m" 40. Page 20, line 250, "As shown in these figure, the higher the temperature, the greater the displacement of the PMMA beam. These show the ability of the PMMA beam to resist loads decreases with increasing temperature and time." has revised to "Based on Figure. 6, the displacement of the beam increases with temperature. It indicates the ability of the PMMA beam to resist load decreases with increasing temperature and time." 41. Page 22, line 289, "The numerical algorithm proposed in this paper could be used to solve other engineering problems. The identification method of the viscoelastic properties is suitable for not only solid materials, but also fluid materials. Dynamic mechanical analysis is always used to characterize the viscoelastic properties of material. The parameters in the fractional viscoelastic behaviour law can be identified according to the experimental investigations. The governing equation of the problem, associated with the proper parameters of material behaviour, could be efficiently and accurately solved with the proposed numerical algorithm." has been added. 42. Acknowledgements has been updated. 43. References have been updated. 44. Appendix have been updated. Highlights: An improved fractional viscoelastic constitutive model is proposed.  The parameters in the proposed viscoelastic model are identified by DMA tests.  An algorithm based on shifted Legendre polynomials is proposed with high accuracy. The fractional governing equation of the viscoelastic beam is solved directly in time domain.  Dynamic analysis of PMMA beam under different load conditions and temperatures.

Figure 1 :

 1 Figure 1: The geometric figure of the viscoelastic beam.

145 3 . 1 .

 31 Shifted Legendre polynomials Legendre polynomials are orthogonal, which used as basic approximate functions. They exhibit simple and convenient form for calculation, compared with other orthogonal polynomials (Chebyshev polynomials...). The interval of Legendre polynomials is -1 to 1. The shifted Legendre polynomials are 150 defined to approximate the unknown function in a larger interval [0, H].

  155

  165

Figure 2 :

 2 Figure 2: Comparison of identified viscoelastic behaviour of PMMA and the experimental date at different temperatures: (a) T = 110 • C, (b) T = 115 • C, (c) T = 120 • C, (d) T = 125 • C

Figure 3 :

 3 Figure 3: The dynamic displacement at different displacements and times.

Figure 4 :

 4 Figure 4: The displacement of PMMA beam under different harmonic loads.

Figure 5 : 5 . 3 .

 553 Figure 5: The displacement of PMMA beam under different harmonic loads when x = 0.8 m.

Figure 6 :

 6 Figure 6: The dynamic displacement y(x, t) at different points under 100 sin(100πt) for (a) T = 110 • C, (b) T = 115 • C, (c) T = 120 • C, and (d) T = 125 • C.

:

  In order to distinguish, We have used k E to represent kinetic energy and T to represent temperature. We keep one thing of A, and the other two things of A are also represented by other symbols.

	results have been added and are marked in cyan in the revised
	manuscript.
	Comment 10: The authors must add some recent studies about using the
	fractional derivatives for the polymeric materials and other complex
	materials in the materials science which are Physica A: Statistical
	Comment 8: Before the algorithm of the numerical method, the authors Mechanics and its Applications Volume 545, 1 May2020,123763 where
	state that they use Matlab software for fitting the coefficients with least the authors studied the viscoelastic with hardening strain with fractional
	square method. The authors should explain if they used the fitting derivatives. Also, Heliyon Volume 6, Issue 7, July2020, e04495 where
	functions found in Matlab software or if they created a program for this the authors studied polymeric materials with fractional derivatives
	purpose. without using Laplace transformation. The authors should add and
	Response: In order to obtain the coefficient matrix of the equation, we discuss these studies. Also, Mechanics of Materials, 103506, 2020 where
	used the numerical algorithm in MATLAB platform. The fractional-order the authors studied the viscoelastic of foam with fractional derivatives.
	equation is transformed to algebraic equations. The least square method
	is used to calculate the solution of the algebraic equations. Thus the
	numerical solution of the fractional order governing equation is obtained.
	Comment 9: In the table-1, the authors should explain why they
	considered the used orders of the derivatives. Also, in the table-2, the
	authors should add more details about the experimental results of the
	elastic modulus of Polymethylmethacrylate and from any study.
	Response: In this paper, n is the number of terms of the shifted Legendre
	polynomial. As n increases, the absolute error gradually decreases,
	indicating the accuracy of the algorithm. Some details about experimental

Table 1 :

 1 Evolution of absolute error in function of n.

	(x, t) Algebraic solution	n = 4	Absolute error value n = 5	n = 6
	(0, 0)	0	0.0140 × 10 -8 0.0913 × 10 -9 0.0236 × 10 -9
	(0, 2)	0	0.0221 × 10 -8 0.1765 × 10 -9 0.0941 × 10 -9
	(0, 4)	0	0.0561 × 10 -8 0.0938 × 10 -9 0.0508 × 10 -9
	(0, 6)	0	0.0488 × 10 -8 0.2287 × 10 -9 0.0765 × 10 -9
	(0, 8)	0	0.0956 × 10 -8 0.5249 × 10 -9 0.0980 × 10 -9
	(1.6, 0)	0	0.0109 × 10 -8 0.0268 × 10 -9 0.0208 × 10 -9
	(1.6, 2)	1.6384	0.0066 × 10 -8 0.0615 × 10 -9 0.0100 × 10 -9
	(1.6, 4)	6.5536	0.0474 × 10 -8 0.0232 × 10 -9 0.0158 × 10 -9
	(1.6, 6)	14.7456	0.0221 × 10 -8 0.1650 × 10 -9 0.0528 × 10 -9
	(1.6, 8)	26.2144	0.0780 × 10 -8 0.2350 × 10 -9 0.0096 × 10 -9

Table 2 :

 2 Identified parameters of the isotropic viscoelastic constitutive behaviour law at different temperatures: T = 110, 115, 120, 125 • C

	T / • C q α /MPa	q 0 /MPa	q 1 /MPa	α
	110		275.8275	582.8400	-0.0002 0.1674
	115		760.7903 -139.9032	0.0001 0.0966
	120 2181.6523 -1939.0033 0.0000 0.0509
	125 1187.7101 -1123.1952 -0.0001 0.0836
	strength/MPa		experimental E experimental E simulative E	′ ′′
	100	200	300	400	500
		f/Hz

′ simulative E ′′

Table 3 :

 3 Input geometric and material parameters of the beam.

	H(m) ρ(kg/m 3 ) A x (m 2 )	I(m 4 )	T ( • C) t(s)
	1.6	1.15 × 10 3	0.02	4.16 × 10 -6	110	60

  MX, Huang J, Song L. Numerical solutions of nonlinear fractional differential equations by alternative Legendre polynomials. Applied Mathematics and Computation, 2018;336:454-464.

	Appendix A. Symbol description
	symbol	explanation
	D α x	Caputo fractional derivative operator
	C, λ	constant
	x	position
	t	time
	p r , q r

  of the shifted Legendre polynomial l n,i (x), L n,i (x), Ln,i (t) shifted Legendre ploynomial A, C, K, U, E, F, L, M, ∆, δ coefficient matrix ϕ n , Φ n , φ n family of shifted Legendre ploynomials T n family of basical ploynomials , Inner product u n n exponent th functional approximation D x first order differential operator matrix for x D mx m exponent th differential operator matrix for x D t first order differential operator matrix for t D vt v exponent th differential operator matrix for t D αt α exponent th differential operator matrix for t

	1 , α	improved fractional constitutive model parameter
	σ(t)	stress
	ε(t)	strain
	ρ	material density
	W	work done by external forces
	y(x, t)	displacement
	E *	complex elastic modulus
	E	storage modulus
	E	loss modulus
	A x	cross-sectional area
	E k	kinetic energy
	V	potential energy
	T	temperature
	I	moment of inertia
	f (x, t)	external load
	H	length
	n	number of terms

Acknowledgements

This work is supported by the Natural Science Foundation of Hebei Province (A2017203100) in China and the LE STUDIUM RESEARCH PROFESSORSHIP award of Centre-Val de Loire region in France. This work has also been supported by the EIPHI Graduate School (contract ANR-17-EURE-0002).

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

Author Contributions: Cao Jiawei wrote some sections of the manuscript and performed analyses; Yiming Chen, Wang Yuanhui, Gang Cheng and Thierry Barrière prepared some other sections of the paper. All authors read and approved the final manuscript.

*Credit Author Statement


Dr Ilya O Aparin

Dr Nina Melnychuk

Dr Andrey S Klymchenko
email: andrey.klymchenko@unistra.fr

Ionic aggregation-induced emission: bulky hydrophobic counterions light up dyes in polymeric nanoparticles

Keywords: dye-loaded polymeric nanoparticles, aggregation-induced emission, aggregationcaused quenching, bulky hydrophobic counterions, single-particle microscopy

Dye-loaded polymeric nanoparticles emerge as powerful bioimaging tools, but their assembly is challenged by aggregation-caused quenching (ACQ) of flat fluorophores. Aggregationinduced emission (AIE) proposes an effective solution against ACQ by exploiting propellershaped neutral fluorophores without polymeric blending. Fighting ACQ of ionic dyes in polymeric nanoparticles can be achieved by bulky hydrophobic counterions. Here, we aim to generate an AIE phenomenon in poorly emissive ionic dyes using bulky counterions. Three cationic hemicyanine dyes of styryl pyridinium family have been synthesized: two planar fluorophores featuring ACQ and one propeller-shaped AIE fluorophore. We found that in water bulky fluorinated tetraphenylborates can light up all three dyes, including planar non-AIE fluorophores. The described "ionic" AIE (iAIE) with bulky counterions enables preparation of 50-nm dye-loaded polymeric NPs showing 40 % quantum yield at 500 mM dye loading, whereas small anions favor poorly emissive large aggregates. Single-particle microscopy reveals that these NPs emit without blinking 50-fold brighter than quantum dots-605 (at 470-nm excitation). iAIE opens the route to assembling charged intrinsically non-AIE fluorophores into bright fluorescent (nano)materials.

Introduction

The surge of fluorescence bioimaging technics encourages the research on new bright dyes and materials with improved properties and fluorescent signal amplification techniques. Fluorescent nanoparticles of small size have been proved to be a promising tool for bioimaging as they largely surpass organic dyes in terms of brightness. [1][2][3][4] Quantum dots, [5][6][7] dye-doped silica nanoparticles, [8,9] nanodiamonds, [10] carbon dots [11] have been intensively studied as bioimaging agents in the last decade. Fluorescent organic nanoparticles assembled from organic dyes is an attractive alternative because they can address the existing challenges in the field, such as poor biodegradability, potential long-term toxicity, limited brightness, etc. [12] In this respect, dye-loaded polymeric nanoparticles are particularly powerful, [13] because polymeric matrix can ensure stability of the nanoparticles preventing organic dyes from bleaching in biological media, while polymer design allows preparation of small NPs with size control in a broad range from 7 till 200 nm. [14,15] However, when the dyes are confined within a nanostructure at very high local concentration, they tend to pi-stack into Haggregates, [16][17][18][19] leading to so-called aggregation-caused quenching (ACQ).

A number of approaches have been developed to fight ACQ in optical materials. [4] The classical approach is to incorporate bulky groups into the fluorophore, which would prevent the dyes form pi-stacking. [20][21][22][23][24] However, this approach requires tedious multi-step synthesis of corresponding dyes and ACQ in these systems is usually not completely prevented. [23] Probably the most popular approach to fight ACQ is aggregation-induced emission (AIE), originally proposed by Tang and co-workers. [25] In this case, the fluorophore generally presents twisted propeller-shaped topology that blocks H-aggregation favoring formation in the solid-state highly emissive aggregates. [26][27][28][29] The effect is further amplified by poor emission of AIE fluorophores (AIEgens) in organic solvents, because their twisted topology favors excitation energy losses by molecular rotation, usually through twisted intramolecular charge transfer (TICT). [30] In the solid state, this rotation is prohibited (so-called restriction of intramolecular rotation, RIR), leading to enhancement of fluorescence emission in the dye aggregates. [26][27][28] AIE expanded to different fields of optical materials, including organic lightemitting diodes, [31] liquid-crystal displays, [32] light-harvesting, [33,34] short-wave nearinfrared, [35] and photo-switching [36] materials, solar concentrators, [37] sensors, [38] theranostics agents, [27] etc, and it significantly impacted development of dye-based nanoparticles. [39][40][41][42][43][44] In particular, AIE enabled preparation of bright fluorescent polymeric nanoparticles, either by direct polymerization of AIE monomers [43,44] or in core-shell systems with the core formed by AIEgen and the shell made of a polymeric biocompatible material. [3,35,[45][46][47] However, preparation of dye-loaded polymeric NPs using AIE dyes meets some challenges, because dispersion of AIEgens within the polymer may lead to dye quenching because its unique crystalline organization is disrupted by the polymer. [48] Moreover, AIEgens are dominated by neutral molecules, whereas the examples of intrinsically charged AIEgens are very rare. [49,50] Charged fluorophores, bearing cationic (or anionic) groups involved in the electronic conjugation, such as cyanines and rhodamines, [51] are particularly attractive because they are characterized by strong electronic delocalization leading to fluorophores of high extinction coefficient and fluorescence quantum yield. However, these dyes are among the most prone to ACQ because they are extended pi-conjugated structures. [16][17][18][19] Moreover, their small inorganic counterions cannot prevent dye-dye pi-stacking, probably because of strong solvation (hydration) and small size of these counterions, as well as dissociative behavior of dye + /anion -pair. These properties do not allow the small inorganic counterions to separate charged fluorophores within the aggregates, and, therefore, the dyes tend to form strongly quenched H-aggregates. [16] In the early studies the groups of Yao [52,53] and Warner [54,55] described the possibility to minimize ACQ in these charged dyes using bulky hydrophobic counterions, based on tetraphenylborate (TPB) and bistriflimide families. These approaches have been applied to variety of fluorophores and counterions. [56,57] However, the NPs of pure dye + /counterion -salts are not stable in biological environment, because they tend to dissolve in hydrophobic compartments of cells showing leakage of the dyes. [57,58] In 2014, we introduced a concept for preparation of dye-loaded polymeric NPs using ion pairs of cationic rhodamine dye (R18) with bulky hydrophobic counterions. [59] Bulky TPBs, especially their perfluorinated analogues (F5-TPB and F12-TPB, Fig. 1), are able to prevent ACQ in polymeric nanoparticles, [57,59,60] because these anions serve as spacers between the cationic dyes and thus prevent their pi-stacking. On the other hand, bulky hydrophobic counterions are essential for encapsulation of the dyes inside the polymer matrix without further leakage into the biological medium. [60] The concept has been extended to cyanines, which enables controlling NPs color allowing multicolor long term tracking of living cells. [61] It has been also found compatible with different types of counterions, including different fluorinated TPBs, tetrakis(trifluoromethyl)borate, aluminium complex and borohydrides, [60,62] showing that the concept of bulky counterion is rather universal and not limited to one type of anion. Importantly, bulky counterions due to ionic forces also brings the dyes together inside the polymer matrix, which leads to fast dye-dye communication by excitation energy migration. [59,63] The latter results in giant light-harvesting phenomenon, where 10,000 dyes can transfer energy to a single energy acceptor. [64] This nanoantenna has already been applied for single molecule detection at ambient sunlight conditions [64] and for the development of ultrabright nanosensors for amplified detection of nucleic acids. [65] So far, it remains unclear whether this counterion approach can be applied to poorly emissive ionic dyes, in order to produce the AIE phenomenon. Recently, we showed that counterions can improve emission properties of AIEgens, especially in the dye-loaded polymeric NPs. [50] In this respect, the current challenge is to apply this bulky counterion approach to produce AIE phenomenon for dyes that intrinsically do not exhibit AIE behavior. In literature, several types of dyes, poorly emissive in organic solvents, have been shown to enhance their emission after association with bulky anions in aqueous media. [52,53,56] However, these systems showed limited fluorescence quantum yields and they have not been really applied for preparation of stable dye-loaded polymeric NPs.

As a promising family of cationic dyes that are poorly emissive in organic solvents without clear AIE characteristics, we selected styryl pyridinium (SP) dyes. Their strongly dipolar nature and flexible conjugation makes them poorly emissive in organic solvents, being highly sensitive to the environmental properties, such as polarity and viscosity. [66][67][68][69] Moreover, lipophilic analogues of SP dyes are practically non-fluorescent in water because of strong ACQ, without showing AIE character. Their fluorogenic (light up) response to viscous environments has been extensively used to design fluorescent probes for biological membranes. [68][69][70][71] A confinement of an SP dye within anionic metal-organic framework have also led to fluorescence light up. [72] However, only one rare report have shown fluorescence enhancement of an SP dye in water due to a bulky anion effect. [73] So far, no attempts have been made to tune AIE properties of SP dyes and they have not been used to fabricate stable polymeric NPs. Thus, SP dyes are promising unexplored dye family for testing our hypothesis.

In the present work, we synthesized three new SP dyes having controlled ACQ/AIE profile: two planar dyes and a propeller shaped molecule, featuring modulated ACQ and AIE characteristics. We found that bulky hydrophobic counterions based on perfluorinated TPB induced drastic fluorescence light up in water. The latter was observed for all SP dyes, including planar derivatives without AIE characteristics. Thus, AIE phenomenon can be observed for non-AIE dyes purely due to counterion effect. This ionic aggregation-induced emission (iAIE) has been applied for preparation of bright dye-loaded polymeric NPs, featuring high fluorescence brightness (50-fold brighter than QDdot-605) and good photostability without detectable blinking. The observed iAIE phenomenon provides a bridge between two complementary concepts of fighting ACQ, AIE and bulky counterions, and opens the route to assembling charged intrinsically non-AIE fluorophores into bright fluorescent (nano)materials.

Results and discussion

Design, synthesis and fluorescent properties of SP dyes

We designed three SP dyes bearing diethylamino, diphenylamino and carbazolyl auxochrome groups: EtSP18, PhSP18 and CzSP18, correspondingly (Fig. 1a). Octadecyl chain was introduced to increase hydrophobicity of the dyes to favor their aggregation and encapsulation into polymeric NPs. [60,74] EtSP18 dye has a planar chromophore structure capable to dissipate excitation energy with rotation of diethylamino moiety and phototoisomerisation. Due to strong polarizability of EtSP18, typical for hemicyanine dyes, it is expected to undergo strong dipole-dipole interaction and aggregation-caused quenching in solid state. To render this dye some AIE character, diethylamino group was replaced with a propeller-shaped diphenylamine moiety (PhSP18), which is extensively used in AIEgens. [27,28] Then to block the AIE property, we designed CzSP18, where we replaced diphenylamine moiety with carbazolyl one rendering the dye perfectly planar and thus bringing back ACQ through pi-stacking. Comparison of CzSP18 and PhSP18 would allow us to evaluate the role of AIE character of the dye in the counterion-controlled solid-state emission. SP dyes PhSP18, EtSP18, CzSP18 were synthesized by adopting previously developed methodology (Fig. 1a). [75] Briefly, 4-picoline was quaternized with octadecyl iodide in acetonitrile under reflux with excellent yield of the picolinium salt. Knoevenagel reaction of N-octadecyl picolinium iodide with 4-(N,N-diphenyl)benzaldehyde, 4-(9carbazolyl)benzaldehyde or 4-(N,N-diethylamino)bezaldehyde was performed in boiling acetonitrile with piperidine as a catalyst, yielding PhSP18, CzSP18, EtSP18, respectively. The obtained dyes were recrystallized either from absolute ethanol or acetonitrile to obtain pure crystalline compounds. All the compounds in form of iodide [I -] salts were characterized with 1H, 13C NMR and HRMS. Then, using simple two-phase ion exchange followed by flash-chromatography purification, we obtained individual ion pairs of the cationic dyes with anions that drastically differ by their diameter: (i) small inorganic anions, perchlorate [ClO4 -] (~0.33 nm) and hexafluorophosphate [PF6 -] (~0.35 nm) and (ii) bulky hydrophobic anions, Next, we studied photophysical properties of EtSP18, PhSP18 and CzSP18 dyes as iodide salts in organic solvents of low and high viscosity: acetonitrile and ethylene glycol, respectively. In the following sequence, CzSP18PhSP18EtSP18, absorption maximum of the dyes in both solvents shifted to the longer wavelengths, which was accompanied by an increase in the extinction coefficient. These bathochromic shifts reflect the increase in the electron donor ability of the aniline nitrogen from CzSP18 to EtSP18. Emission spectra showed quite different tendencies. Narrowest emission band and smallest Stokes shift were observed for EtSP18, which reflect more pronounced cyanine-like electronic configuration with delocalized positive charge in the ground state, typical for styryl pyridinium dyes. [66] All three dyes showed poor fluorescence quantum yield (QY) in acetonitrile, whereas in more viscous glycerol a significant QY increase was observed (Fig. 2b). This result is a clear indication that all three dyes are quenched by intramolecular rotation, similar to molecular rotors [76,77] and AIE dyes. [26][27][28][29] In both solvents, EtSP18 showed higher QY values compared to CzSP18 and PhSP18, probably because of its more pronounced cyanine-like nature. Absorption and emission spectral properties of all three dyes were independent of dye concentration (Figs. S1-S9), indicating that in this studied concentration range the effects of dye aggregation in these solvents could be neglected. All three dyes bear highly hydrophobic groups and therefore are expected to aggregate in water, which should enable us to study AIE properties of these dyes in the presence of different counterions. Dye aggregation was done first by diluting acetonitrile solution of dye iodides in water, followed by addition of 10-mol excess of corresponding salts of [

Aggregation-induced emission of dyes in water: effect of bulky counterions

I -], [ClO4 -], [PF6 -], [TPB], [F5-TPB],
[F12-TPB] in order to achieve an ion exchange. This protocol was successfully used by us and others earlier to prepare fluorescent NPs. [52,53,78] In all studied combinations of three dyes with corresponding counterions, we observed formation of aggregates in the size range 30-90 nm, according to dynamic light scattering (Table S1).

However, the dye salts of iodide, perchlorate and hexafluorophenylphosphate gave aggregates of high polydispersity and they were unstable, showing precipitation within hours/days. By contrast, aggregates of the dyes with [TPB], [F5-TPB] and [F12-TPB] gave stable nanoparticle suspensions of systematically low polydispersity (0.1-0.2) with sizes 40-80 nm (Table S1). Previous studies suggested that bulky hydrophobic anions can provide a stabilizing effect by adsorption at the surface of newly formed nanoparticles. [52,57] Next, we studied optical properties of aggregates of the dyes with the small inorganic anions in water. Low fluorescence signal (QY < 1%) was observed for all dye + /[I -] salts probably due to iodide quenching properties (Fig. 3c). In case of [ClO4 -] and [PF6 -] counterions, EtSP18 and CzSP18 dyes remained poorly emissive, which was not the case of PhSP18 with QY reaching 6.5% and 11% for [ClO4 -] and [PF6 -], correspondingly. Thus, as expected, the diphenylamino propeller group rendered PhSP18 AIE properties, whereas other two SP dyes did not show pronounced AIE character. Nevertheless, fluorescence QY values for PhSP18 with inorganic counterions were quite low compared to classical AIE fluorophores, [26] which suggested a competition between AIE and ACQ phenomena in their aggregates. In the absorption spectra, strong hypochromism was observed for the aggregated dyes in water (Fig. 3a, Figs. S10-S12). In case of EtSP18, this hypochromosim was accompanied by strong blue shift of the absorption band (Fig. 3a), indicating formation of H-aggregates with inorganic anions, in line with the observed low QY values. For all three dyes the absorption and emission spectral properties did not change with dye concentration in the range 0.5-10 µM (Figs. S13-S18), indicating that in these conditions the aggregation behavior was concentration-independent.

Very different emission properties were observed in the presence of bulky hydrophobic counterions [F5-TPB] and [F12-TPB]. Indeed, the latter drastically enhanced fluorescence intensity and QY of all three dyes in water compared to those with inorganic anions (Fig. 3). Parent non-fluorinated TPB showed intermediate behavior (Fig. 3, S10-S12): it enhanced QY of EtSP18 in water compared to inorganic anions, but it was inefficient in case of both PhSP18 and CzSP18. This result shows that in addition to the size of the counterion, fluorination also plays an important role to ensure efficient fluorescence of dyes in the aggregated state. Moreover, for all three dyes with [F5-TPB] and [F12-TPB] anions, fluorescence was much stronger compared to corresponding solutions in acetonitrile (Fig. 3a), suggesting that they behave like typical AIE fluorophores. The AIE phenomenon could be seen by a bare eye: when solutions of CzSP18/[ClO4 -] and CzSP18/[F5-TPB] in acetonitrile were compared to those in water, only CzSP18/[F5-TPB] showed strong emission (Fig. 3e).

In acetonitrile/water solvent mixtures, CzSP18/[F5-TPB] showed strong increase in the fluorescence intensity at higher water fractions, whereas for CzSP18/[ClO4 -] the effect was much less pronounced (Fig. 3d). These observations are remarkable because EtSP18 and CzSP18 do not exhibit pronounced AIE behavior with inorganic anions, but become typical AIEgens in the presence of bulky counterions. On the other hand, PhSP18 has intrinsic AIE properties due to propeller-shaped diphenylamino group, so that all its salts were emissive in water and the differences between small and bulky anions were much less pronounced (Fig. 3c and Fig. S11). Nevertheless, combination of intrinsic AIEgen properties with the effect of bulky counterion made PhSP18/[F5-TPB] and PhSP18/[F12-TPB] salts the most emissive among the studied series (Fig. 3c). This additive effect is very important for design of highly emissive materials and it is supported by our recent work, where emissive properties of an AIEgen could be further improved by using these bulky counterions. [50] In the absorption spectra of the SP dye salts, some hypochromism was still observed compared to the acetonitrile solutions, especially for EtSP18 (Figs. S10-S12). The absorption and emission spectral properties for all three dyes with bulky did not vary with dye concentration in 0.5-10 µM range (as for salts with small inorganic anions), confirming that in the studied conditions the aggregation behavior was concentration-independent (Figs. S13-S18). In addition, absorption bands of all three dyes with fluorinated TPBs were red shifted compared to salts with small anions without significant change in the absorption band shape (Fig. 3a, Figs. S10-S12). The latter could be related to the drop of environment polarity within the aggregates, which is expected to induce bathochromic shift in SP dyes exhibiting negative solvatochromism. [79] Better understanding of the unique effect of bulky counterions on the optical properties of the dye salts could be achieved by X-ray analysis of eventual crystals of the dyes with [F5-TPB] and [F12-TPB]. Unfortunately, these salts were present as amorphous frozen ionic liquid-like materials, which failed to crystallize.

The observed light up of non-AIE dyes in water with assistance of bulky hydrophobic counterions can be called ionic aggregation-induced emission (iAIE). In contrast to classical AIE, this phenomenon does not require AIEgen. On the other hand, it differs from the previously described ion-association approach, which normally use highly fluorescent dyes such as rhodamine [57] and cyanine [52,78] . Some reports showed that poorly fluorescent cationic dyes could enhance their fluorescence with addition of TPB counterions, [50,52] but here we show that this phenomenon can be observed for both AIE and non-AIE fluorophores. However, the limitation of NPs based on neat dye salts is their tendency to dissolve in hydrophobic compartments of cells leading to dye leakage, as it was shown for rhodamine B -R18 salts with [F5-TPB] and [F12-TPB] counterions. [57] Therefore, of particular interest is to encapsulate them into rigid polymer matrix as it has already been realized for conventional bright fluorescent dyes. [15,59,61] Following the encouraging data with the neat dye salts, we further studied behavior of SP dye salts encapsulated into polymeric NPs. Poly(methyl methacrylate-co-methacrylic acid (PMMA-MA) NPs formulated by charged-controlled nanoprecipitation protocol, [14] where polymer solution in acetonitrile with different amounts of the dye was quickly diluted into phosphate buffer with vigorous stirring (Fig. 4a). Commercial PMMA-MA forms small spherical nanoparticles at pH 7.4 stabilized with negative charge of carboxylic groups. [14] Following well-established protocols of ion exchange, [60] we first prepared salts of these three dyes with different counterions. Importantly, dyes with bulky hydrophobic counterions exhibited much higher mobility (Rf) on thin layer chromatography than those with small inorganic anions (Fig. S35, Table S2), indicating that the former dye salts are significantly more hydrophobic. Dyes with fluorinated TPB counterions were loaded into polymer matrix of NPs by simply adding them to acetonitrile solution before the nanoprecipitation. A wide range of dye loadings from 20 mM to 500 mM with respect to PMMA-MA polymer was tested, yielding small NPs with average size of 50-65 nm, which increased slightly with loading, according to DLS measurements (Fig. 4g). By contrast, NPs of dyes with small inorganic anions loaded above 50 mM gave large aggregates (Fig. 4g) that precipitated over time. Transmission electron microscopy (TEM) images of NPs at high loading of EtSP18, PhSP18, and CzSP18/[F5-TPB] (500 mM) showed homogeneous spherical NPs with average size of 52, 50 and 55 nm, respectively. However, NPs formulated with PhSP18/[ClO4 -] showed aggregates already at 50 mM and shapeless precipitate at 500 mM loading (Fig. 5a,b). Bulky TPB counterions are crucial for effective loading of cationic dyes into polymeric NPs due to strong hydrophobic stabilization of ion pairs and shielding of dye positive charge, whereas dyes salts with small inorganic anions can dissociate and undergo ion exchange in aqueous medium. [60] The latter leads to direct interaction of negatively-charged PMMA-MA with cationic dyes and neutralization of PMMA-MA negative surface charge and aggregation (Fig. 5c). [59] Light absorption of NPs increased gradually with dye loading in the range from 20 mM to 500 mM for the ion pairs with bulky counterions (Fig. 4b and Figs. S4-S18). By contrast, for small anions, saturation of absorbance was observed already at 50 mM, followed by further drop of absorbance at 200 mM, which is probably related to particle aggregation and further precipitation (Fig. 4e and Figs. S4-S18). Fluorescence signal also increased with dye loading for the three dyes with bulky counterions, whereas for small anions rapid signal saturation was observed already at 50 mM loading, followed by strong drop at 200 mM (Fig. 4d). Remarkably, fluorescence quantum yields for PhSP18 and CzSP18 NPs with bulky counterions were very high (30-50%), with only a small decrease at the highest (500 mM) loading to 30% and 40%, respectively. The high QY values were observed for both [F5-TPB] and [F12-TPB] anions, whereas for small inorganic anions they were significantly lower even for 100 mM loading (Figs. S19-S21). Thus, bulky counterions seem to play a crucial role in lighting up dyes and preventing their ACQ inside polymeric NPs. Importantly, they enhanced emission for both AIEgen PhSP18 and dyes without AIE character (EtSP18 and CzSP18). In pure salts and in polymer NPs with small anions, PhSP18 show systematically higher QYs values (Figs. 3c and4d), which is probably because of intrinsic AIE property of the dye. However, when bulky counterions are used in polymeric NPs, a non-AIE dye CzSP18 is most emissive (Fig. 4d). This indicates that bulky counterion can light up a poorly emissive dye inside a polymer matrix independently of the AIE character of the dye. In fact, the intrinsic AIE property of the dye does not play important role when it is encapsulated into the hydrophobic matrix, because polymer interferes with dye aggregation and favors uniform distribution in the matrix, as it was shown previously for other AIE dyes. [48] On the other hand, counterion seems to overcome this limitation as by itself it can assemble the dyes inside the polymer minimizing the ACQ. [59,60] To confirm whether the dyes are indeed assembled inside polymeric NPs in a close proximity we made fluorescence anisotropy measurements. This method enables direct access to extent of dye-dye excitation energy migration, which reflects qualitatively the dye-dye proximity. [63,80] It was found that increase in the dye loading for all three dyes with [F5-TPB] counterion led to a drop in the fluorescence anisotropy reaching nearly zero values at 500 mM, indicating a drastic decrease in the dye-dye spacing. It should be noted that the anisotropy loss is usually observed for dye assemblies with small Stokes shift (21 nm / 686 cm -1 for rhodamine B in ethanol), where homo-FRET process can take place easily. [23,63] In the present case, it takes place for dyes with much larger Stokes shift (Table 1) and thus poorer spectral overlap for homo-FRET. Therefore, observation of nearlyzero anisotropy and high QY values at 500 mM loading suggests that the dyes are probably aggregated inside the polymer matrix, but the ACQ is blocked by the counterion, preventing the dyes from direct pi-stacking. Among the three tested dyes, the loss of anisotropy was fastest in case of EtSP18, which is probably because it has the smallest Stokes shift and thus can easier undergo energy hopping through homo-FRET processes than the other dyes. These studies suggest that iAIE phenomenon can be also observed inside polymeric matrixes and thus can be directly applied for preparation of highly fluorescent polymeric NPs. To the best of our knowledge, it is the first report where the counterion approach was successfully used in combination with poorly emissive non-AIE fluorophores to prepare highly emissive dye-loaded polymeric NPs.

Dye-loaded polymeric NPs

Single-particle fluorescence

Next, we examined the developed NPs by single-particle imaging using wide-field fluorescence microscopy. The new NPs were compared with commercial Quantum dots 605 (Qdots 605), presenting emission maximum in similar wavelength range. First, we prepared polymeric NPs loaded with 5 mM, 50 mM and 500 mM loading of PhSP18/[F5-TPB] to quantify single-particle brightness and its distribution (Fig. 6a). After 1000-times dilution of polymeric NPs and Qdots 605 with a phosphate buffer, they were immobilized at the polyethyleneimine-covered glass surface and studied by fluorescence microscopy. Singleparticle brightness increased gradually with the dye loading in polymeric NPs (Fig. 6b).

Images of NPs with lowest PhSP18/[F5-TPB] concentration (5 mM) showed less fluorescent spots, which basically could be explained by presence of empty polymeric NPs without fluorescent dyes. Image analysis revealed that NPs encapsulating 500 mM PhSP18/[F5-TPB] were 50-fold brighter than Qdots 605 at 470 nm excitation wavelength, optimal for both NPs (Fig. 6c). We should note here that this comparison is affected by significantly larger size of our particles compared to Qdots 605 (~50 nm vs ~25 nm). On the other hand, the advantage of our system compared to Qdots is that the size of fluorescent core in our NPs can be tuned to much larger extent within 9-100 nm [14,15] We also compared single-particle fluorescence of NPs loaded with EtSP18, PhSP18 and CzSP18 (Fig. 7a). In line with spectroscopy data in cuvette, PhSP18 NPs demonstrated the highest brightness per particle (Fig. 7b). Video-rate imaging of NPs revealed single-particle traces with good photostability (Fig. 7c) and absence of blinking (Fig. 7d). Interestingly, highest photostability upon irradiation was observed for CzSP18 NPs, which revealed intermediate fluorescence brightness (Fig. 7c). Absence of blinking is an advantage of the newly developed particles compared to Qdots, which showed expected whole particle blinking (Fig. 7e). It is also advantageous compared to our brightest dye-loaded polymeric NPs based on rhodamine (R18) and [F5-TPB], which showed intermittent or complete blinking process. [63] The stable emission of the newly developed NPs is probably related to much larger Stokes shift of this dye, which slows down the energy migration. Even though the energy migration takes place in these NPs, according to fluorescence anisotropy data, it is probably much slower, thus decreasing the dye cooperativity and thus the whole-particle blinking. [4,59,60,63]

Conclusion

We report here an "ionic" aggregation-induced emission (iAIE) promoted by formation of stable ion pairs of poorly fluorescent non-AIE dyes and their application for preparation of bright dye-loaded polymeric nanoparticles. Three hydrophobic derivatives of styryl pyridinium dyes have been developed: two planar fluorophores (EtSP18 and CzSP18), bearing diethylamino and carbazolyl moieties (non-AIE dyes), and propeller shaped dye PhSP18 with diphenylamino group (AIE dye). The dyes showed poor fluorescence in organic solutions, but enhance their emission in viscous media preventing rotational quenching. We found that bulky fluorinated tetraphenylborate counterions [F5-TPB] and [F12-TPB] light up fluorescence of all three dyes in water, whereas with small inorganic anions [PF6 -] and [ClO4 - ] only AIE dye PhSP18 shows some fluorescence enhancement. Thus, we observe "ionic" AIE phenomenon that can operate with intrinsically non-AIE charged dyes. Fluorescence QY of pure PhSP18 (AIE) and CzSP18 (non-AIE) salts with bulky counterions reaches 23% and 14%, respectively, suggesting that the counterion-based iAIE and AIE phenomena are additive. We applied our findings to preparation of dye-loaded polymeric NPs based on biocompatible PMMA-MA. For all three dyes, bulky counterions enable preparation of bright dye-loaded polymeric NPs, featuring ~50 nm size and good fluorescence quantum yields reaching 40 % at 500 mM dye loading. The quantum yield increased in the following order: EtSP18 < PhSP18 < CzSP18, indicating that inside polymeric NPs it is independent of the intrinsic AIE character of dyes, but controlled by the bulky counterions. In sharp contrast, attempts to encapsulate the dye salts with small inorganic anions [I -], [ClO4 -] and [PF6 -] into polymeric NPs gave large aggregates with low quantum yield at high dye loading. Thus, thanks to iAIE with bulky counterions, for the first time polymeric NPs loaded with poorly fluorescent dyes reach 40% quantum yield at 500 mM dye loading, rendering NPs high brightness. Single-particle measurements revealed that polymeric NPs loaded with PhSP18 with [F5-TPB] counterions are 50-fold brighter than quantum dots Qdots 605, excited at 470 nm. The new NPs show stable non-blinking fluorescent signal, which is an additional advantage compared to QDots. Another interesting feature of these NPs is their large Stokes shift (~ 200 nm), which is useful to minimize background noise due to scattering and autofluorescence.

Thus, the described iAIE phenomenon show that weakly emissive dyes presenting or not AIE characteristics can light up in the aggregated form in water or polymeric matrix thanks to bulky hydrophobic counterions. As the studied dyes present molecular rotor properties, they turn on their emission probably because of restriction of intramolecular rotation (RIR) in the aggregated form, as it was reported earlier for AIE dyes. [26][27][28] However, bulky hydrophobic counterions can light up also non-AIE molecular rotor dyes with flat architecture (e.g. CzSP18) probably because they serve as a spacer between the dyes in the aggregated state, as it was observed for highly emissive rhodamine and cyanine dyes. [57,59,60] Comparative studies revealed that counterion fluorination is also very important for efficient emission in the aggregated form, which is probably related to (i) negative charge delocalization among large number of fluorines [81,82] and (ii) super-hydrophobicity of fluorinated compounds. [83] High hydrophobicity and delocalized charge of the anion prevents its solvation and dissociation of the dye salts in water. Therefore, in contrast to highly hydrated inorganic anions that allow flat cationic dyes to pi-stack in water, poorly solvated bulky hydrophobic anions interact strongly with cationic dyes through electrostatic and hydrophobic forces, thus preventing them from pi-stacking. Moreover, high hydrophobicity of these counterions ensures efficient encapsulation of cationic dyes inside hydrophobic matrix of polymeric NPs, preserving their strong emission, in line with our previous works on highly emissive fluorophores based on rhodamine and cyanines. [60][61][62] Overall, the described iAIE phenomenon provides a bridge between AIE and bulky hydrophobic counterions approaches to fight ACQ. On one hand, iAIE extends the counterion concept to large variety of ionic dyes presenting weak fluorescence. On the other hand, it shows the way to generate bright AIE materials based on non-AIE fluorophores using specially selected bulky counterions.

Experimental Section

Materials and methods, synthesis of the styryl pyridinium dyes and their ion-exchange with different anions, the details of nanoparticle preparation, single-particle fluorescence measurements and TEM measurements are described in the Supporting Information.

Supporting Information

Supporting Information is available from the Wiley Online Library or from the author.

  tetraphenylborate [TPB] (~1.08 nm), tetrakis(pentafluorophenyl)borate [F5-TPB] (~1.16 nm) and tetrakis[3,5-bis-(1,1,1,3,3,3-hexafluoro-2-methoxy-2-propyl) phenyl]borate [F12-TPB] (~1.6 nm) (Fig. 1b).

Figure 1 .

 1 Figure 1. Design of SP dyes and nanoparticles. (a) Synthesis of EtSP18, PhSP18 and CzSP18 dyes; (b) Counterions explored in this study; (c) Schematic representation of pure dye NPs and dye-loaded polymeric NPs.

Figure 2 .

 2 Figure 2. Spectroscopic studies of new dyes iodides in organic solvents of different viscosity. Absorption (a) and fluorescence (b) spectra of EtSP18, PhSP18 and CzSP18 in acetonitrile and ethylene glycol (λexc = 400±2 nm for CzSP18 and 460±2 nm for PhSP18 and EtSP18). In the fluorescence spectra, the curves were divided by absorbance values at the excitation wavelength. Concentration of dyes was 10 and 1 µM for absorption and fluorescence spectroscopy, respectively.

Figure 3 .

 3 Figure 3. Neat dye salts in aqueous media. a. Absorption and fluorescence spectra of 5 µM CzSP18, PhSP18 and EtSP18 in acetonitrile and NPs in water with [ClO4 -] and [F5-TPB] counterions (λexc = 400±2 nm for CzSP18 and 460±2 nm for PhSP18 and EtSP18). Fluorescence spectra were divided by absorbance values at the excitation wavelength; b.

Figure 4 .

 4 Figure 4. Dye-loaded polymeric NPs. a. Schematic representation of dye-loaded PMMA-MA NPs precipitation. b. Absorption and fluorescence spectra of PMMA-MA NPs loaded with PhSP18/[F5-TPB] at different loading. c. Integral fluorescence intensity of PMMA-MA NPs with different loading of the dyes with [F5-TPB] and [ClO4 -] counterions. d. Fluorescence quantum yields of dye-loaded PMMA-MA NPs. e. Absorption and fluorescence spectra of PMMA-MA NPs loaded with PhSP18/[ClO4 -] at different loading. f. Anisotropy values of polymeric NPs loaded with 5 mM, 50 mM 500 mM of the dyes per PMMA-MA. g. Average size of polymeric NPs and bulky aggregates by dynamic light scattering.

Figure 5 .

 5 Figure 5. TEM imaging of dye-loaded polymeric NPs. a. TEM images of polymeric NPs samples loaded with 500 mM CzSP18, EtSP18 and PhSP18 dyes with [F5-TPB] counterion; representative NPs with 50 mM loading of PhSP18/[ClO4 -]; bulky shapeless aggregates of PMMA-MA with 500 mM PhSP18/[ClO4 -]. b. Size distribution histograms of NPs corresponding to the TEM images (top). c. Schematic presentation of NPs encapsulated with PhSP18/[F5-TPB] and PhSP18/[ClO4 -].

Figure 6 .

 6 Figure 6. Single-particle characterization of NPs at different loading of PhSP18/[F5-TPB] and reference quantum dots 605 (Qdots 605) by wide-field fluorescence microscopy. a. Single-particle fluorescence images of polymeric NPs loaded with 5 mM, 50 mM and 500 mM PhSP18/[F5-TPB] and Qdots 605. For visibility, images of Qdots 605 and polymeric NPs with 5 mM PhSP18 are presented with 10-fold intensity amplification and NPs with 50 mM loadingwith 2-fold amplification. b. 3D brightness histograms referred to the corresponding fluorescence images (top). c. Corresponding distribution histograms of NPs brightness. Wide-field fluorescence microscopy was performed with excitation at 470 nm (power density of 0.51 W cm -2 ); emission was collected at >500 nm with a long-pass filter.

Figure 7 .

 7 Figure 7. Single-particle characterization of polymeric NPs loaded with 500 mM EtSP18/[F5-TPB], PhSP18/[F5-TPB] and CzSP18/[F5-TPB]. a. Single-particle fluorescence images of dye-loaded NPs. b. Corresponding distribution histograms of NPs brightness. c. Photobleaching profiles of polymeric NPs loaded with 500 mM CzSP18, EtSP18 and PhSP18 dyes with [F5-TPB] counterion. d. Single-particle traces of the NPs. e. Single-particle trace of Qdots 605.

Table 1 .

 1 Spectroscopic properties of EtSP18, PhSP18 and CzSP18 iodides in acetonitrile and ethylene glycol. a

		Solvent ɛ, M -1 *cm -1 λabs, nm λemis, nm	Stokes shift, nm/cm -1	Φf, %
	EtSP18	ACN EG	51800 56600	488 495	620 615	132 / 4363 120 / 3942	2.61 31.9
	PhSP18	ACN EG	42300 46200	464 472	725 680	261 / 7759 208 / 6481	0.61 4.6
	CzSP18	ACN EG	24700 27600	399 406	710 640	311 / 10978 0.47 234 / 9006 2.2
	a Concentration of dyes was 10 and 1 µM for absorption and fluorescence spectroscopy,
	respectively.

Acknowledgements I.O.A. acknowledges Marie Curie post-doctoral research grant (H2020-MSCA-IF-2018, DNANanoProbes 846571) by the European Union. This work was supported by the European Research Council ERC Consolidator grant BrightSens 648525. The authors also thank C. Crucifix from the FRISBI platform for help with electron microscopy.

Conflict of Interest

The authors declare no conflict of interest.

The table of contents entry

Keyword: ionic AIE Ilya O. Aparin, Nina Melnychuk, and Andrey S. Klymchenko*

Ionic aggregation-induced emission: bulky counterions light up dyes in polymeric nanoparticles

A phenomenon of ionic aggregation-induced emission (iAIE) is described, which combines concepts of AIE and bulky hydrophobic counterions to fight aggregation-caused quenching in cationic weakly emissive fluorophores. iAIE has been applied for preparation of bright dyeloaded polymeric nanoparticles exhibiting ~50-fold stronger emission without blinking compared to quantum dots 605 in single-particle microscopy.