R Chakir

J K Hammond

A non-intrusive reduced basis method for elastoplasticity problems in geotechnics

Keywords: Reduced Basis method, Finite Element method, Parametric studies, Elastoplasticity, Soils

This work aims at investigating the use of reduced basis (RB) methods to diminish the cost of the resolution of parameter-dependent partial differential equations (PDEs) present in elastoplasticity problems arising from geotechnics modeling. Computation times for large three-dimensional analysis commonly take tens of hours, making optimization procedures or sensitivity analysis, relying on repeated simulations, hardly feasible. In many cases the analysis requires very specific features such as highly non-linear constitutive laws, involving a complex description of hardening phenomena in soils, which could not be solved using a standard reduced basis method. Given this constraint, an approach making it possible to use the reduced basis framework with any finite element software (without modifying the code) and considering it as a "black box" gives the so-called non-intrusive reduced basis method a versatility of great practical interest. Our approach involves the computation of less expensive (but less accurate) FE approximation during the online stage and improvement of those solutions using a RB-based rectification method. The chosen application belongs to the field of tunnel engineering, the particular problem being the displacement of the soil around a shallow tunnel for varying values of parameters: Young's modulus, friction angle, cohesion coefficient (characterizing the soil) and confinement loss (caused by the excavation of the tunnel). CESAR-LCPC, a FEM-based software, was used as a black-box to compute displacements during the offline and the online stages, whereas Freefem++ was used for the implementation of the RB-based rectification method and analysis of the results. We proposed two rectification methods in the non-intrusive framework, and found that a modified rectification method was more adapted to the problem considered. With this non-intrusive method computation time has been reduced by 85% compared to P 2 finite element method without loss of accuracy.

Introduction

Numerical modeling has met growing success over the last decades, becoming indispensible in the field of geotechnical engineering, leading to the resolution by finite elements of even larger nonlinear problems. This trend stems from the need to account for the influence of constructing new structures, such as deep foundations of high-rise buildings or shallow tunnels for transport infrastructures, on neighboring structures (e.g. sewers, existing buildings, etc.) in dense urban areas. Computation times for large three-dimensional analysis commonly take tens of hours, making sensitivity analysis relying on repeated simulations hardly feasible. A common approach is to develop simplified models, such as metamodels, to approximate the model without significant loss of accuracy. In [START_REF] Bolzon | An effective computational tool for parametric studies and identification problems in materials mechanics[END_REF] a metamodel based on Proper Orthogonal Decomposition (POD) with radial basis functions (RBF) was applied to test problems in material mechanics with the goal of illustrating the capability of these metamodels to reproduce mechanical responses to the loading of complex non-linear material systems. An extended version of the POD-RBF metamodel was proposed in [START_REF] Khaledi | Robust and reliable metamodels for mechanized tunnel simulations[END_REF] to surrogate a 3D finite element simulation of a tunnel using Hardening Soil model.

Another approach to rapidly compute reliable approximations of solutions of complex problems with many parameters is reduced basis (RB) methods [START_REF] Prud'homme | Reliable real-time solution of parametrized partial differential equations: Reduced-basis output bound methods[END_REF]. These methods rely on the parametric structure of the model and that when the parameters vary, the manifold of all possible solutions can be approximated by a low-dimensional space, the reduced basis space. The reduced basis is constructed from solutions of the parametrized problem for a well-chosen set of parameters. Standard reduced basis methods are Galerkin approximations of the full order model within a lower-dimensional reduced basis space. One of the keys of RB techniques is the decomposition of the computational work into offline and online stages. During the offline stage the reduced basis functions are computed, as well as all parameter-independent quantities. This is done only once, whereas parameter-dependent quantities are computed during the online stage. Application of the reduced basis method to linear elastic solid mechanics problems with parameters of different natures (either physical or geometrical) was proposed in [START_REF] Hoang | A fast, certified and "tuning free" two-field reduced basis method for the metamodelling of affinely-parametrised elasticity problems[END_REF][START_REF] Huynh | Reduced basis approximation and a posteriori error estimation for stress intensity factors[END_REF][START_REF] Liu | A novel reduced-basis method with upper and lower bounds for real-time computation of linear elasticity problems[END_REF][START_REF] Milani | Reduced basis method for linear elasticity problems with many parameters[END_REF][START_REF] Veroy | Reduced-basis methods applied to problems in elasticity: Analysis and applications[END_REF]. The efficiency of the reduced basis or POD-based reduction methods relies on liberating online calculation costs from dependency on the discretization. However in elasto-plastic problems with highly nonlinear behavior, not uncommon in the field of soil mechanics, the computational complexity related to the local integration of the nonlinear constitutive laws is not reduced. Several alternative ways to carry out the standard POD-based reduction method for problems with nonlinear behavior were investigated. For example in [START_REF] Corigliano | Model order reduction and domain decomposition strategies for the solution of the dynamic elastic-plastic structural problem[END_REF][START_REF] Radermacher | Model reduction in elastoplasticity: proper orthogonal decomposition combined with adaptive sub-structuring[END_REF] a partial reduction is performed over the region of the domain with elastic behavior, while the plastified region remains unreduced. This selective POD-based model reduction was extented by an adaptive method of sub-structuring POD(A-SPOD) in which the subdomain where model reduction is applied is determined automatically. In [START_REF] Ryckelynck | Hyper Reduction of finite strain elasto-plastic models[END_REF][START_REF] Ryckelynck | Multi-level a priori hyper-reduction of mechanical models involving internal variables[END_REF][START_REF] Ryckelynck | Multidimensional a priori hyperreduction of mechanical models involving internal variables[END_REF] a hyper-reduction approach was proposed by Ryckelynck to treat the problem of local dependency and extended by Zhang [START_REF] Zhang | Efficient hyper reduced-order model (HROM) for parametric studies of the 3d thermo-elasto-plastic calculation[END_REF] to a thermo-elastoplastic model. The hyper reduction method consists in introducing reduced integration domains for internal variables. However these methods require modification of the finite element calculation code leading to an intrusive procedure, which is particularly restrictive in the case of the considered geotechnics modeling applications. Analysis of the displacement field around a tunnel opening using numerical techniques is quite sensitive to the constitutive models of the soil used to described the fundamental behavior of the materials involved. In many cases the analysis requires very specific features which are not available in all finite element softwares, such as highly non-linear constitutive laws, involving a complex description of hardening phenomena in soils. Given this constraint, an approach making it possible to use reduced basis methods with any finite element software (without modifying the code) and considering it as a "black box" gives the so-called non-intrusive reduced basis method a versatility of great practical interest. Our approach involves the computation of less expensive FE approximation (but less accurate) with a black-box FE software and improvement of those solutions using a reduced basis during the online stage. In this work, we aim to demonstrate the feasibility of this approach -a two-grid finite-element/RB method, introduced in [START_REF] Chakir | A two-grid finite-element/reduced basis scheme for the approximation of the solution of parameter dependent PDE[END_REF][START_REF] Chakir | Une méthode combinée d'éléments finis à deux grilles/bases réduites pour l'approximation des solutions d'une EDP paramétrique[END_REF] -to geotechnics modeling.

This paper is organized as follows. In Section 2, we formulate the elastoplastic problem providing a brief description of the physical system, the material behavior laws, the governing equations and boundary conditions. In Section 3, we provide a brief introduction to reduced basis methods and discuss a preliminary analysis of the feasibility and reliability of reduced basis approximations of the elastoplastic problem. In Section 4 the problem is solved with a non-intrusive reduced basis method.

Finally Section 5 presents the conclusions.

The elastoplastic problem

The chosen application belongs to the field of tunnel engineering. In urban areas, it is necessary to consider the possible impact of the construction of a shallow tunnel on existing structures (buildings, foundations, etc.). In many cases, the first step consists of evaluating the settlements that would be induced by the construction of a tunnel in a "greenfield " environment, i.e. with no structure built at the surface. Let us consider a circular tunnel built through a 50-m horizontal ground layer. The tunnel diameter is D=10 m and the axis depth is H=25 m. The analysis is carried out under the plane strain assumption. If the ground is homogeneous and isotropic, only half of the ground layer needs to be considered. For practical reasons, the analysis is limited to a distance of L=100 m from the tunnel axis (see figure 1).

From material behavior laws to the governing equations

In design calculations, materials (soil, concrete, rock, metal, liquid, gas) are considered as continuous mediums (or continua). These materials are thus considered to obey certain general physical and mechanical principles, such as the conservation of energy and momentum. Everyday experience can tell us that different materials do not behave in the same way under the same forces. General physics laws do not allow us to make the distinction between different sorts of materials. We therefore want to characterize the specific behavior of the continuum equivalent to the material under consideration. This is the goal of the constitutive laws associated to a material; the laws must characterize the evolution caused by given exterior forces and be specific to the material in question. When switching from one material to another, the laws must translate the differences in practically observed behavior. The constitutive law associated to a material is necessary to complete the system of equations of any mechanics problem of continua or structural design. The behavior of the soils in our problem is represented by an elastoplastic model used for pulverulent soils (sands) and for long-term coherent soils (clay and silt). Observations show that irreversible deformations appear when the stress exceeds a certain level. Let u be the displacement vector; the deformation is assumed to be infinitesimal so that the strain tensor can be written as (u) = 1 2 (∇u + t ∇u). The framework of plasticity is based on the assumption that the strains can be split into the sum of two terms :

= e + p , ( 1
)
where e is the elastic strain tensor and p is the plastic part of the total strain tensor , which corresponds to the irreversible part of the strain. The elastic part of the behavior of the soil is linear 2 and isotropic and described by Hooke's law (characterized by Young's modulus E and Poisson's coefficient ν). The plastic part of the soil's behavior is considered nonlinear and is obtained via the Mohr-Coulomb model [START_REF] Coquillay | Prise en compte de la non linéarité du comportement des sols soumis à de petites déformations pour le calcul des ouvrages géotechniques[END_REF] (characterized by the cohesion c, the friction angle ϕ, and the dilatancy angle ψ).

2.1.1. Linear elastic behavior : Hooke's Law Hooke's law [START_REF] Khaledi | Robust and reliable metamodels for mechanized tunnel simulations[END_REF] describes the relationship between the stress tensor σ(u) ∈ R d×d and the elastic strain tensor e (u) ∈ R d×d .

σ(u) -σ 0 = E ν (1 + ν)(1 -2ν) tr e (u) I d + E (1 + ν) e (u) (2)
with σ 0 the initial stress tensor, E and ν soil's parameters.

Nonlinear plastic behavior : Mohr Coulomb's model

It is assumed that the plastic strain does not evolve as long as the stress tensor remains in the interior of a region of the stress space, called the elastic domain [START_REF] Coquillay | Prise en compte de la non linéarité du comportement des sols soumis à de petites déformations pour le calcul des ouvrages géotechniques[END_REF][START_REF] Mestat | Lois de comportement des géomatériaux et modélisation par la méthode des éléments finis[END_REF]. The elastic domain is generally defined by a condition of the type f (σ ij ) < 0, where f is called the yield function. The yield surface is the boundary of the elastic domain and thus defined by f (σ) = 0. For sands, the yield function can be expressed as follows.

f

(σ ij ) = (σ -σ s ) -(σ s + σ ) sin ϕ -2c cos ϕ, (3)
where σ and σ s represent the largest and smallest eigenvalues of the stress tensor σ (often call principal stresses in mechanics). The parameters ϕ and c are the friction angle and the cohesion characterizing the soil. In this study, we focus on the case of elastic-perfectly plastic models, in which the yield surface does not evolve with loading. Let us consider the stress tensor σ ij corresponding to a given load. If f (σ ij ) < 0, then σ ij is in the elastic domain, and so we have that the deformation variation is described simply by d = d e .

If f (σ ij ) = 0, then σ ij is on the boundary of the elastic domain. To describe the behavior at this point, we need to know if the material is in loading, in which case the deformation variation is described by

d = d p + d e .
or if the material is unloading and has an elastic behavior. At a regular point σ ij of the elasticity boundary, the plastic deformation can be described by the so-called "plastic flow rule"

d p = d λ ∂g ∂σ ,
where d λ ≥ 0 is a scalar called the plastic multiplier and g is given by

g(σ ij ) = (σ -σ s ) -(σ s + σ ) sin ψ -2c cos ψ. ( 4
)
The problem is closed by the "consistency condition"

df d λ = 0. (5)

Equilibrium equation

We consider a static process, then the equilibrium describing our system is:

div(σ) + ρF = 0, (6)
where ρF = 0 -γ is the external body force and γ is the volumetric weight of the soil.

The elastic deformation is linked to the variation of the stress by a linear relation:

σ -σ 0 = C : e . ( 7
)
From ( 2) one can see that C ijkl (representing the elasticity tensor of the material) is constant, symmetrical

(C ijkl = C jikl = C ijlk = C klij )
and depends only on E and ν.

Boundary conditions

In this paper, we consider a 2D problem on the bounded domain Ω ⊂ R 2 (see figures 2 and 3). We will impose zero horizontal displacement on Γ 1 and Γ 3 and zero horizontal and vertical displacement on Γ 2 . The load consists of a surface density of force I applied on the wall of the tunnel (Γ 6 ) calculated from the initial stress tensor (which we assumed geostatic):

I = λ σ 0 • n, with σ 0 = x 2 K 0 γ 0 0 x 2 γ ,
where λ represents the confinement loss caused by the excavation of the tunnel and K 0 = 1 -sin(ϕ) the coefficient of the earth pressure at rest. Let (u i,i=1,••• ,d ) be the displacement components in the x i -directions in Ω. The boundary conditions of our problem read as follows:

                 σ • - → n = -λ σ 0 • - → n on Γ 6 , u 1 = 0 on Γ 1 and Γ 3 , (σ • - → n ) 2 = 0 on Γ 1 and Γ 3 , u i,i=1,2 = 0, on Γ 2 , σ • - → n = - → 0 on Γ 4 , (8)

Weak formulation and finite element approximation

Associated to the physical domain Ω, we defined the following functional space

X = v ∈ (H 1 (Ω)) 2 v 1 = 0 on Γ 1 , Γ 2 , Γ 3 v 2 = 0 on Γ 2 (9)
We then introduced the weak form of our elastoplastic problem, arising from ( 6) and ( 8):

Find u ∈ X such that, Ω (u) -p (u) : C : (v) dΩ = Ω ρF v dΩ - Ω σ 0 : (v) dΩ - Γ 6 σ 0-→ n • v dΓ, ∀v ∈ X (10) with ij (u) = 1 2 ∂u i ∂x j + ∂u j ∂x i .
In what follows, we fix the Poisson coefficient at ν = 0.3, the volumetric weight of the soil at γ = 20kN/m 3 and assume that the dilatancy angle is equal to the friction angle (ψ = ϕ).

Let us denote by µ = (E, λ, ϕ, c) our parameter set and by D ⊂ R 4 our parameter domain. We will decompose the left hand side of equation ( 10) into a linear term:

a e (u(µ), v; µ) = Ω (u(µ)) : C(µ) : (v) dΩ
and a nonlinear term

a p ( p (u(µ)), v; µ) = Ω p (u(µ)) : C(µ) : (v) dΩ,
and denote by L(v; µ) the right-hand side term

Ω ρF v dΩ - Ω σ 0 : (v) dΩ - Γ 6 σ 0-→ n • v dΓ.
We

(µ) ∈ X such that, ∀v ∈ X, a e (u(µ), v; µ)-a p ( p (u(µ)), v; µ) = L(v; µ). ( 11
)
Let {T h } h be a family of regular triangulations of Ω and denote by X h the following P k finite element space

X h = {v = (v 1 , v 2 ) ∈ X, ∀K ∈ T h , v i | K ∈ P k (K)}.
The finite element discretization of ( 11) is as follows : for a given µ ∈ D, find u h ∈ X h (µ) such that,

a e (u h (µ), v h ; µ) -a p ( p (u h (µ)), v h ; µ) = L(v h ; µ), (12)
In this work, CESAR-LCPC [START_REF] Humbert | CESAR-LCPC: A computation software package dedicated to civil engineering uses[END_REF], a FEM-based software, has been used to solve [START_REF] Ryckelynck | Multi-level a priori hyper-reduction of mechanical models involving internal variables[END_REF], which employs the following iterative procedure (see algorithm 1) to approximate the displacement u h , the stress tensor σ h , the strain tensor h and the plastic strain tensor p h . For more details on the computational procedure see [START_REF] Coquillay | Prise en compte de la non linéarité du comportement des sols soumis à de petites déformations pour le calcul des ouvrages géotechniques[END_REF][START_REF] Mestat | Lois de comportement des géomatériaux et modélisation par la méthode des éléments finis[END_REF].

Algorithm 1 : Resolution in cesar-lcpc 1: Initialization p 0 = 0, σ 0 = σ 0 2: for k = 1 to Nb iter max do 3: Compute displacement u k h such that 4: a e (u k h , v h ) = a p ( p k-1 , v h ) + L(v h ) (13) 5:
Compute

σ * k = C : (u k h ) 6:
for all K h ∈ T h do

7:
for all integration points x q do 8:

if f (σ * k (x q )) > 0 then 9: Compute λk = f (σ * k ) ∂f (σ * k ) ∂σ : C : ∂g(σ * k ) ∂σ 10: Compute k p = λk ∂f ∂σ (σ k-1 h ) 11: Compute σ k h = σ * k -C : p k 12:
else 13:

Set k p = k-1 p and σ k h = σ * k 14:
end if

err = a e (u k h , v h ) -a p ( p k , v h ) -L(v h ) 19:
if err < tol then 20:

Break 21: end if 22: end for 23: Set u h = u k h , σ h = σ k h , p h = k p
We want to apply reduced basis methods within this framework to compute the displacement u h corresponding to different values of E, λ, ϕ, c.

Methodology

The resolution of the problem introduced in the previous section can prove costly, particularly in the many-query context due to its parametric nature, making it an ideal candidate for reduced basis methods. The reduced basis method relies on the fact that when the parameters vary, the set of solutions is often of small Kolmogorov dimension, implying that M h = {u h (µ) ∈ X h | µ ∈ D}, the manifold of all solutions, can be approximated by a finite set of well-chosen FE solutions of the parametrized PDE.

One can identify a set of parameters, S

N = (µ 1 , µ 2 , • • • , µ N ) ∈ D N such that the particular solutions (u h (µ 1 ), • • • , u h (µ N ))
will generate this low dimension space. The idea of reduced basis methods is to compute an inexpensive and accurate approximation, u N h (µ), of the solution to problem [START_REF] Ryckelynck | Hyper Reduction of finite strain elasto-plastic models[END_REF] for any µ ∈ D by seeking a linear combination of the particular solutions (u h (µ 1 ),

• • • , u h (µ N )) u N h (µ) = N i=1 α h i (µ) u h (µ i ). (14)
For a stable implementation of the reduced basis method, it is necessary to build a better basis than the one composed of the {u h (µ i )} 1≤i≤N , usually by a Gram-Schmidt method. In what follows, we denote by {ξ 1 , • • • , ξ N } these L 2 orthonormalized basis functions, and by X N h the approximation space which they engender: the reduced basis space. During the implementation of the reduced basis method, the computational work is separated into two stages: offline and online. This decomposition is a key ingredient of the method. The reduced basis functions, {ξ 1 , • • • , ξ N }, as well as all expensive parameter-independent terms are computed once during the offline stage and stored, whereas during the online stage -for each new value of the parameters -inexpensive parameter-dependent quantities are evaluated, together with the computation of the coefficients α h i (µ) . The usual RB method is a Galerkin method on the space X N h , which is of much smaller dimension than the original approximation space X h , the resolution of the problem (13) in X N h is less expensive than in the true finite element space X h . However, to perform the online stage efficiently, one must isolate the parametric contribution to the corresponding linear system, allowing all parameter-independent matrices and vectors to be built only once and saved during the offline stage. In the case of Mohr-Coulomb's model used in CESAR-LCPC, a parameter-dependent term must be calculated at each integration point of the mesh (see Algorithm 1) during the iterative procedure implemented to solve [START_REF] Ryckelynck | Multi-level a priori hyper-reduction of mechanical models involving internal variables[END_REF]; it is hence impossible to free the online stage of the FE complexity. This entirely nullifies the advantages of the RB method applied to our model. To overcome this flaw, we propose to use an alternative, less intrusive method, introduced in [START_REF] Chakir | Une méthode combinée d'éléments finis à deux grilles/bases réduites pour l'approximation des solutions d'une EDP paramétrique[END_REF][START_REF] Chakir | A two-grid finite-element/reduced basis scheme for the approximation of the solution of parameter dependent PDE[END_REF], where coarse FE approximations are computed during the online stage, then projected into the reduced basis space and improved by a rectification technique.

We will begin by considering an analysis of the feasibility of RB methods for our problem (section 3.1), and will then discuss the non-intrusive method in more detail in section 3.2.

POD analysis

In order to determine if model reduction approaches, such as reduced basis methods or proper orthogonal decomposition (POD), can be applied to this problem, we will try to evaluate the complexity of the manifold M h of all possible solutions induced by varying parameters. This analysis consists in a singular value decomposition method applied to the correlation matrix of solutions of (12) computed for different values of the parameters. Once the rapid decay rate of the singular values is confirmed, one can assume that RB method is worth implementing. Using cesar-lcpc to compute P 1 and P 2 -FE solutions of ( 12) for varying values of µ ∈ Ξ test -a parameter set with sample size of N test = 525 selected over the parameter domain D -a correlation matrix of L 2 -norm scalar products (u h (µ i ), u h (µ j )) L 2 (Ω), 1≤i,j≤Ntest was computed. An L 2 -orthonormalized POD basis was constructed using the following eigenfunctions

w k = 1 √ λ k Ntest =1 v k ( ) u h (µ ) 1 ≤ k ≤ N test , ( 15
)
where v k ( ) represents the th component of the k th eigenvector of the correlation matrix when ordered by decreasing eigenvalues.

Let P P OD k be the L 2 -orthogonal projection operator from X h into the space X k,P OD h , spanned by the k first POD basis functions w k . Each test solution {u h (µ i ), µ i ∈ Ξ test } was projected onto X k,P OD h to analyze the ability of the POD basis to approach the manifold M h , depending on the number of POD modes. In figure 4 we can see the associated errors plotted along with the eigenvalues of the matrix, where the average error is

1 N test Ntest i=1 u h (µ i ) -P P OD k u h (µ i ) L 2
and the maximal error is We observe that the eigenvalues decay rapidly and that the projection errors are quite small. Figure 5 shows the average POD projection errors with respect to so-called reference solutions, computed on a very fine mesh for parameter values not included in the sample space solutions. These errors are plotted with the fine FEM error for comparison. Figure 6 displays the corresponding maximal errors. We can see that with only k = 5 POD modes, the POD projection errors reach the same level of accuracy as the P 1 FEM errors. As for the P 2 FE errors, we only need about k = 10 POD nodes. This suggests that a reduced basis approach is worth implementing.

u h (µ k max ) -P P OD k u h (µ k max ) L 2 with µ k max = argmax µ i ∈Ξtest u h (µ i ) -P P OD k u h (µ i ) L 2 .
u ref (µ) -P POD k u h (µ) L 2 u ref (µ) -u h (µ) L 2
u ref (µ) -P POD k u h (µ) L 2 u ref (µ) -u h (µ) L 2
u ref (µ k max ) -P POD k u h (µ k max ) L 2 u ref (µ k max ) -u h (µ k max ) L 2
u ref (µ k max ) -P POD k u h (µ k max ) L 2 u ref (µ k max ) -u h (µ k max ) L 2

A non-intrusive reduced basis method : two-grid FE/RB method with a rectification approach

A popular strategy for constructing a reduced basis in the case of parameter-dependent problems is to use greedy algorithms, based on the idea of selecting the locally optimal element at each step. This option can be seen as an alternative to the POD strategy of the previous section. If we have an appropriate a priori error estimator to avoid full resolution of the problem to compute the test solutions, the greedy algorithm can be very low-cost. Knowing the Kolmogorov dimension of the solution space is relatively small, we can fix a maximum number N g of basis functions to be computed by the Greedy algorithm (given below, algorithm 2). Additionally, for stable implementation the chosen basis functions are L 2 -orthonormalized with a Gram-Schmidt method.

Algorithm 2 : Greedy's algorithm to build the reduced basis space 1: Initialization: given

Ξ test = (µ 1 , . . . , µ ntest ) ∈ D ntest , n test >> 1 2: Choose randomly µ 1 ∈ D 3: Set S 1 = {µ 1 } and X 1 h = span(u h (µ 1 )). 4: for N = 2 to N g do 5: µ N = argmax µ∈Ξtest u h (µ)-P N -1 u h (µ) L 2 u h (µ) L 2
(where

P N -1 is the L 2 -orthogonal projection operator from X h into X N -1 h ) 6: S N = S N -1 ∪ µ N 7: X N h = X N -1 h + span(u h (µ N )) 8: end for
Standard reduced basis methods are based on a Galerkin approach, whereas the two-grid FE/RB method involves the computation of a less expensive FE solution and improvement of this solution using the reduced basis. The standard reduced basis method aims at evaluating the coefficients α h i (µ) intervening in the decomposition ( 14) of u N h (µ), which can appear as a substitute to the optimal coefficients

β h i (µ) = (u h (µ), ξ i ) L 2 (16)
corresponding to the decomposition of the L 2 -projection of u h (µ) into the space X N h . Let {T H } H be a family of "coarse" regular triangulations of Ω, such that H >> h; we denoted by X H the coarse FE approximation space associated to this mesh, and by u H (µ) the coarse FE approximation of ( 12) on X H . The alternative two-grid FE/RB method consists in proposing another surrogate to the coefficients

β h i (µ) defined by β H i (µ) = (u H (µ), ξ i ) L 2 , (17)
While coarse FE approximations can be computed quickly enough to be used in the online stage, they may not be accurate enough for practical use. As the computation of u H (µ), for H >> h, is significantly less expensive than that of u h (µ), with the mesh size H (chosen adequately) the coefficients β H i (µ) can be used to compute a low-dimensional approximation :

N i=1 β H i (µ) ξ i . ( 18
)
Let P H N be the L 2 -projection operator from X H into the space X N h . Considering that we have used embedded FE spaces, namely X H ⊂ X h , we have

P N u H (µ) = P H N u H (µ)
, and in consequence we will simplify the notation by also denoting by P N the L 2 projection operator P H N . To improve even further the accuracy of this technique we propose to perform a rectification of the P N u H (µ). This is so far an empirical approach, which leads to great improvements in practice. A first explanation of the successful post-processing strategy first presented in [START_REF] Chakir | A two-grid finite-element/reduced basis scheme for the approximation of the solution of parameter dependent PDE[END_REF] and then used in [START_REF] Herrero | RB (Reduced basis) for RB (Rayleigh-Bénard)[END_REF] in the framework of reduced basis simulation of PDE's can be found in [START_REF] Hernandez | Quelques contributions vers la simulation parallèle de la cinétique neutronique et la prise en compte de données observées en temps réel[END_REF]. This treatment will ensure that for the parameters {µ i } 1≤i≤N used in the construction of the reduced basis, the method returns exactly β h i (µ i ). In practice, we want to identify a so-called rectification matrix R N associated to the transformation R N such that :

R N P N u H (µ i ) = P N u h (µ i ) ∀ 1 ≤ i ≤ N.
Since {β h j (µ i )} 1≤j≤N and {β H j (µ i )} 1≤j≤N are the optimal coefficients intervening in the decomposition of P N u h (µ i ) and P N u H (µ i ), the standard matrix, denoted by A N , associated to the transformation R N is equal to

A N = B N h × B N H -1 with A N ∈ R N ×N ,
where

B N h =    β h 1 (µ 1 ) • • • β h 1 (µ N ) . . . . . . . . . β h N (µ 1 ) • • • β h N (µ N )    and B N H =    β H 1 (µ 1 ) • • • β H 1 (µ N ) . . . . . . . . . β H N (µ 1 ) • • • β H N (µ N )    .
Let us note that, contrarily to the u h (µ), which we don't want to compute for a large number of values of µ, the true solutions u h (µ i ) have already been computed to build the reduced basis, making the computation of A N relatively cheap. For each new value of µ, the coefficients β H i (µ) will be replaced by

N k=1 A N ik β H k (µ)
, and an improved two-grid FE/RB approximation to equation ( 12), for R N = A N , can be :

R N P N u H (µ) = N i,j=1 A N ij β H j (µ) ξ N i . ( 19
)
In our problem, we noticed that A N was rather poorly conditioned, and propose here a pre-processing to improve the rectification. Instead of computing the coefficients from the fine and coarse RB solutions, we will consider the previously computed POD basis functions to construct another rectification matrix K N . To do so, in addition to the POD basis function w k introduced in the previous section, we introduced "coarse" POD basis function

w H k = Ntest =1 v k ( ) u H (µ ) 1 ≤ k ≤ N.
We defined a pre-processing matrix

D N = F N h × F N H -1 ,
where

F N h =    (w 1 , ξ 1 ) L 2 • • • (w N , ξ 1 ) L 2 . . . . . . . . . (w 1 , ξ N ) L 2 • • • (w N , ξ N ) L 2    and F N H =    (w H 1 , ξ 1 ) L 2 • • • (w H N , ξ 1 ) L 2 . . . . . . . . . (w H 1 , ξ N ) L 2 • • • (w H N , ξ N ) L 2    .
We then construct the new rectification matrix K N as follows, for a suitable N max . Nmax) . By "cutting off" the chosen rectification before significant increases in the condition number (at N max ), we can prevent associated peaks in error, thus achieving the results of K N . Figure 7 shows condition numbers for the three proposed matrices: A N , D N , and K N . Figure 8 shows rectification errors for the three proposed rectification matrices. We can see that the matrix D N is better conditioned than the matrix A N and that the rectification process is improved. However the most significant improvements are seen with matrix K N .

K N = D Nmax 0 0 T N , with T N = 1 N    1 0 . . . 0 1    ∈ R (N -Nmax)×(N -

Numerical experiments

The above-described method was applied to the problem using CESAR-LCPC for the FE resolution of equation [START_REF] Ryckelynck | Multi-level a priori hyper-reduction of mechanical models involving internal variables[END_REF] and FreeFem++ [START_REF] Hecht | New development in freefem++[END_REF] was used for the implementation of the two-grid FE/RB method and analysis of the results. Three meshes were considered: a coarse mesh T H for the inexpensive computation of coarse solutions u H (µ), a fine mesh T h for the computation of satisfactory solutions used in the construction of the reduced basis, and a reference mesh T ref considered fine enough to provide true solutions used for error calculation. See figure 9 below. A parameter set Ξ trial with sample size of N trial = 16 is selected over D \ Ξ test to test our method with P 2 FEM grids. While in some applications, the simple rectification with R N = A N will achieve the desired results, in this case the significant variation between coarse and fine solutions used to build the rectification matrix caused inadequate rectification results. We thus used matrix R N = K N introduced in the previous section to improve the rectification. Figure 10 shows rectification errors during the offline stage. In figure 11, we can see the two-grid reduced basis method errors using rectification matrix R N = K N , for N = 16; the error reaches the same order of precision as the P 2 -FEM fine solutions. We note that while rectification error in figure 10 does not descend further for N ≥ N max , in contrast to the fine projection errors during the offline stage, figure 11 shows that the rectification approximation online does attain the same precision as the fine FEM solution. Figure 12 shows the actual displacement for a given parameter, µ = µ max = argmax µ∈Ξ trial u h (µ) -R N P N u H (µ) L 2 = (125, 0.35, 23, 0.03). The application of this problem being to evaluate impact on surface structures, we can consider displacement at the surface to be a quantity of interest. Figure 13 shows error maps with respect to the P 2 -FE approximation over the calculation domain various N -values of the two-grid FE/RB method with and without the rectification, where the parameter value µ max = (125, 0.35, 23, 0.03) corresponds to the solution with maximal error. We can see that the errors of the rectified solution with respect to the non rectified solution. Figure 14 shows errors over the calculation domain for N = 15 with respect to the very fine reference solution, again for µ max = (125, 0.35, 23, 0.03). We can see that the rectified solution errors closely resemble the fine FEM errors. 16 show computation times for finite element simulations and the proposed online reduced basis method. We can see that satisfactory results can be obtained in a total of 3.17s over the full domain, a reduction by 85% of computation time compared to a fine finite element approximation. In the case of many-query approximations -such as parametric studies, and possibly optimization procedures which are currently too computationally expensive for practical use -this reduction would prove to be significant.

||u ref (µ) -P N u H (µ) || L 2 ||u ref (µ) -P N u h (µ) || L 2 ||u ref (µ) -R - N P N u H (µ) || L 2 ||u ref (µ) -u h (µ) || L 2 ||u ref (µ) -u H (µ) || L 2

CPU Times

Coarse FEM

2.41s

Two-Grid RB/FE

3.17s

Fine FEM 20.22s

Conclusions

In this paper we proposed a non-intrusive reduced basis method for application to the parametrized PDEs governing an elastoplasticity problem which could not be solved using a standard reduced basis method. We demonstrated the small dimension of the solution space affiliated to the problem using POD analysis. We then proposed two rectification methods in the non-intrusive framework, and found that a modified rectification method was more adapted to the problem considered. The particular problem being the displacement of the soil around a shallow tunnel, the displacement at the surface approximated by the reduced model was considered, showing the successful approximation results held true when considering only the most important area of the domain. The results of this study demonstrate the feasibility of the presented two-grid non-intrusive reduced basis method in geotechnics modeling, a domain for which reduced modeling techniques can provide great benefit. Specifically, this technique is well-adapted to the particular PDE problem studied considering its non-intrusive nature.

Figure 1 :

 1 Figure 1: Geometry of the physical domain

Figure 2 :

 2 Figure 2: 2D representation of the domain Ω

Figure 3 :

 3 Figure 3: Boundaries of the domain Ω

Figure 4 :

 4 Figure 4: Relative errors of the POD projection with P1-FE (left) and P2-FE (right) snapshots.

  k (number of POD modes)

  k (number of POD modes)

Figure 5 :

 5 Figure 5: Average error of the POD projection vs FEM error using P1-FE (left) and P2-FE (right) snapshots

  k (number of POD modes)

  k (number of POD modes)

Figure 6 :

 6 Figure 6: Maximal error of the POD projection vs FEM error using P1 (left) and P2 (right) snapshots

Figure 7 :Figure 8 :

 78 Figure 7: Condition number of the different rectification matrices: AN , DN , and KN during the offline stage (P2 FEM solutions)

P 2

 2 ndof = 1247 P 2 ndof = 4853 P 2 ndof = 19143

Figure 9 :

 9 Figure 9: Coarse (left), fine (middle) and reference (right) embedded meshes used to compute FEM solutions

Figure 10 :

 10 Figure 10: Average (left) and maximal (right) rectified RB projection errors on test space during the offline stage with RN = KN and Nmax = 16

  of the reduced basis space)

  of the reduced basis space)

Figure 11 :

 11 Figure 11: Average (left) and maximal (right) rectified RB projection errors on trial space during the online stage with RN = KN and Nmax = 16

Figure 12 :

 12 Figure 12: Displacement value for µmax = (125, 0.35, 23, 0.03)

Figure 13 :

 13 Figure 13: Relative error maps of the two-grid FE/RB approximation without (left) and with (right) rectification as function of N for µ = µmax = (125, 0.35, 23, 0.03)

Figure 14 :Figure 15 :Figure 1

 14151 Figure 14: Error maps for N=15 and µmax = (125, 0.35, 23, 0.03)

Figure 16 :

 16 Figure 16: Comparaison of calculation times of P2-FE and two-grid FE/RB Methods

  h (µ max ) -P N u H (µ max ) || L 2 ||u h (µ max ) -P N u h (µ max ) || L 2 ||u h (µ max ) -R N P N u H (µ max ) || L 2

							||u h (µ) -P N u H (µ) || L 2
	norm	10 -1					||u h (µ) -P N u h (µ) || L 2 ||u h (µ) -R N P N u H (µ) || L 2	norm	10 -1
	2								2
	Relative error measured in L	10 10 -3 -4 10 -2							Relative error measured in L	10 -4 10 -3 10 -2
		10 -5								10 -5
		5	10	15	20	30	35	40	50	5	10	15	20	30	35	40	50
				N (dimension of the reduced basis space)					N (dimension of the reduced basis space)

||u

Table 1 :

 1 Comparaison of calculation times of P2-FE and two-grid FE/RB Methods (for N = 15)

Let us note that here the term "linear" or "nonlinear" refers to the behavior of the material, not necessarily to a linear or nonlinear equation.

Acknowledgment

The authors would like to thank Emmanuel Bourgeois for fruitful discussions and very valuable advice regarding our numerical simulations with CESAR-LCPC. This research did not receive any specific grant from funding agencies in the public, commercial or not-for-profit sectors.


F Maucher

S Skupin

S A Gardiner

I G Hughes

Creating Complex Optical Longitudinal Polarization Structures

In this paper we show that it is possible to structure the longitudinal polarization component of light. We illustrate our approach by demonstrating linked and knotted longitudinal vortex lines acquired upon non-paraxially propagating a tightly focused sub-wavelength beam. Remaining degrees of freedom in the transverse polarization components can be exploited to generate customized topological vector beams.

The concept of light being a transverse wave represents an approximation that is suitable if the angular spectrum is sufficiently narrow [START_REF] Born | Principles of Optics: Electromagnetic Theory of Propagation, Interference and Diffraction of Light[END_REF]. However, many practical applications ranging from microscopy to data storage require tight focusing. Tight focusing implies a broad angular spectrum and the notion of light being transverse becomes inappropriate. Hence, the longitudinal polarization component can typically not be neglected [START_REF] Richards | [END_REF]3]. To mention a few examples, a "needle beam" with particularly large longitudinal component was proposed in [4], and radial transverse polarization permits the significant decrease of the focal spot size [5,6] while the generated longitudinal component may even dominate the interaction with matter [7]. Last but not least, a Möbius strip in the polarization of light was realized in [8].

In addition, there is current substantial interest in "structured light", that is, generating customized light fields that suit specific needs in applications in a range of fields [9][10][11][12]. Since the proposal of the Gerchberg-Saxton algorithm [13] 1972, advances in light shaping [14][15][16] now permit the realization of complex light patterns in the transverse polarization plane, including light distributions the optical vortex lines of which form knots [17][18][19]. Knotted topological defect lines and their dynamics have been studied in diverse other settings, including for example classical fluid dynamics [20][21][22], excitable media [23][24][25], and nematic colloids [26,27]. To date, the approach has typically been to determine the longitudinal polarization component of the electric field from given transverse components, and attempts to target complex structures in the longitudinal component have not yet been pursued. The reason for this is twofold.

On the one hand, the longitudinal component is not directly accessible by beam shaping techniques. On the other hand, non-paraxial beam configurations are required, and topological light is usually studied in paraxial approximation. It is therefore not immediately evident that the whole range of threedimensional light configurations known for transverse components can be realized in the longitudinal component as well.

In this paper, we will show that complex lightshaping of the longitudinal polarization component is indeed possible. To this end, we firstly identify non-paraxial light patterns that give rise to vortex lines that form knots or links. Secondly, we invert the problem and derive how one must structure the transverse components of a tightly focused beam to give rise to a given complex pattern in the longitudinal component, and thus present the first example of non-transverse non-paraxial knots. Finally, we demonstrate that remaining degrees of freedom in the transverse polarization components allow for simultaneous transverse shaping, which could be interesting for applications, e.g., inscribing vortex lines into Bose-Einstein Condensates.

We begin with the equations describing a monochromatic light beam:

∇ 2 E(r ⊥ , z) + k 2 0 E(r ⊥ , z) = 0, (1)
∇ • E(r ⊥ , z) = ∇ ⊥ • E ⊥ + ∂ z E z = 0. ( 2
)
Here, k 2 0 = ω 2 /c 2 = (2π/λ) 2 , and we have introduced the transverse coordinates r ⊥ = (x, y) and transverse electric field components E ⊥ = (E x , E y ) as we consider propagation in the positive z direction. All three components of E in Eq. ( 1) fulfil the same wave equation, and for a given field con- figuration E f (r ⊥ ) at z = 0 (e.g., at focus) the general solution for propagation in the positive

z di- rection reads Ê(k ⊥ , z) = Êf (k ⊥ ) exp(ik z z), where k z (k ⊥ ) = k 2 0 -k 2 ⊥ , k ⊥ = (k x , k y )
and the symbol ˆdenotes the transverse Fourier domain. The prescribed field configuration Êf must obey certain constraints. Firstly, in order to get a solution, there must be no evanescent fields present, is, Êf = 0 for k 2 ⊥ ≥ k 2 0 . Secondly, Eq. ( 2) implies for solutions propagating in the z direction that Êf z (k ⊥ = 0) = 0.

As preparation for what follows, we first investigate how to obtain a non-paraxial tightly focused knot or link in E z , assuming that we can directly prescribe E f z . For the transverse paraxial case, recipes to generate vortex lines in various shapes are known, and they usually involve linear combinations of Laguerre-Gaussian modes [19,28]. These recipes are not directly applicable to our problem, since there are evanescent fields, due to the nonparaxiality the wavelength cannot be scaled away and it would lead to Eq. ( 2) being violated. Nevertheless, we found that it is possible to adopt those recipes for the non-paraxial case by an educated guess. Starting from a given linear combination of Laguerre-Gaussian modes f , filtering in the transverse Fourier domain [START_REF] Goodman | Introduction to Fourier Optics[END_REF],

H k0 (k ⊥ ) = e - 1 2λ 2 √ k 2 ⊥ -k 0 2 , H 0 (k ⊥ ) = 1-e -(3λk ⊥ ) 2 ,
chops off evanescent amplitudes as well as ampli- tudes close to k ⊥ = 0, and the longitudinal polarization component at z = 0 reads

Êf z (k ⊥ ) = f (k ⊥ )H k0 H 0 , for k 2 ⊥ < k 2 0 0 for k 2 ⊥ ≥ k 2 0 . (3)
Since the higher-order Laguerre-Gaussian modes are broader in Fourier space and thus lose relative weight after attenuation, one must decrease the relative amplitudes of the lower-order modes to a certain extend. We have found that the following field structures produce a Hopf link or trefoil, respectively,

f Hopf = 4LG σ 00 -5LG σ 01 + 11LG σ 02 -8LG σ 20 (4) f Trefoil = 9LG σ 00 -20LG σ 01 + 40LG σ 02 -18LG σ 03 -34LG σ 30 , (5)
for wavelength λ = 780 nm and width σ = 370 nm ≈ λ/2 of the usual Laguerre-Gaussian modes

LG σ ij (r ⊥ ). The resulting amplitudes before and after filtering for f being either f Hopf or f Trefoil defined in Eqs. ( 4) and ( 5) are plotted in Fig. 1 after normalization to unity. We note that individual mode amplitudes can be changed by about 10% without altering topology, demonstrating a degree of robustness and hence experimental feasibility.

Let us now verify that the presented patterns in the focal plane in fact give rise to vortex lines with the desired topology. It is straightforward to propagate the filtered component E f z , as defined in Eq. ( 3), in the z-direction. The vortex lines throughout three-dimensional space are depicted by the black lines in Figs. 2(a,b) together with a slice in the z = 0 plane of the light profile phase. The obtained vortex lines are topologically equivalent to a Hopf link and a trefoil, as drawn in the insets.

We now address the main point of this paper, i.e. how to choose the transverse polarization components to obtain a given longitudinal polarization Amplitude and phase for the linearly polarized transverse components Eqs. ( 6) and ( 7 component. Because only the transverse components E x and E y are accessible to beam shaping, this point is also of great practical relevance. When inspecting Eq. ( 2), at a first glance the problem may seem to be ill-posed, given that only the longitudinal derivative of the longitudinal polarization enters, i.e., ∂ z E z . However, in Fourier space it is easy to see from Eq. ( 2) that a linearly polarized solution to this problem is given by

E f x = -i x -∞ F -1 k z Êf z (x , y)dx , E f y = 0, ( 6
)
where F -1 [ĝ](x, y) = g(x, y) denotes the inverse transverse Fourier transformation. Obviously, an orthogonally solution also exists,

E f x = 0, E f y = -i y -∞ F -1 k z Êf z (x, y )dy . (7)
Both x and y polarized solutions Eqs. ( 6) and ( 7) evaluated for Hopf link and trefoil are depicted in Fig. 3. It is noteworthy that any superposition of real and imaginary parts of the solutions Eqs. ( 6) and ( 7) admissible, as long as the coefficients of this superposition add up to one. Furthermore, an arbitrary solenoidal field can be added without having an effect on the longitudinal component. We will discuss this later in more detail.

Unfortunately, the transverse polarization components computed from Eqs. ( 6) and ( 7) are impractical, since, even though E f z has finite support, the components E f

x or E f y are non-zero on a semi-infinite interval (see Fig. 3). However, simply attenuating these components by multiplying with, e.g., a sufficiently wide super Gaussian profile SG w N (r ⊥ ) = exp(-r 2N ⊥ /w 2N ) allows the resolution of the problem of semi-infinite light distributions without affecting propagation of the longitudinal component close to the optical axis. Evaluating ∇ ⊥ • SG w N (r ⊥ )E f ⊥ (r ⊥ ) reveals that, where ∇ ⊥ SG w N is large and points in the direction of E f ⊥ , additional satellite spots in the longitudinal component will appear. We have checked that using e.g. a super Gaussian with N = 5 and w = 10λ ensures that these additional spots are sufficiently far from the region of interest and both Hopf link and trefoil develop in the propagation of the modified E z component.

So far, we have seen that the answer to the problem of how to choose E f ⊥ (r ⊥ ) for realizing a prescribed E f z is not unique, and there are certain degrees of freedom in the choice of E f ⊥ . The fundamental theorem of vector calculus (Helmholtz decomposition) allows us to decompose a (sufficiently wellbehaved) vector field F into an irrotational (curlfree) and a solenoidal (divergence-free) vector field, and F can be written as F = -∇φ + ∇ × A. We wish to apply this theorem to the transverse plane, that is, we set F = E f ⊥ (r ⊥ ). In this case, the decomposition reduces to

E f ⊥ (r ⊥ ) = -∇φ(r ⊥ ) + ∇ × [A(r ⊥ )e z ] , (8)
with e z the unit vector in z direction. It is straightforward to verify that

φ(k ⊥ ) = -i k z (k ⊥ ) Êf z (k ⊥ ) k 2 ⊥ (9)
gives rise to a valid transverse polarization component E f ⊥ . The scalar function A(r ⊥ ) may be chosen arbitrarily, because the term

E sol ⊥ = ∇ × [A(r ⊥ )e z ]
produces the solenoidal part of E f ⊥ , which does not give rise to any longitudinal polarization component. The irrotational choice for E f ⊥ , that is, evaluating Eqs. ( 8) and ( 9) with A(r ⊥ ) = 0, for Hopf link and trefoil are shown in Fig. 4.

While the transverse polarization components shown in Figs. 3 and4 produce exactly the same longitudinal field, they are completely different, in particular from a topological point of view. Unlike the light distributions in Fig. 3, which do not contain vortices in the transverse polarization components, the irrotational transverse polarization components in Fig. 4 each feature phase singularities. Furthermore, note that the amplitudes required in the irrotational transverse polarizations are only roughly two to three times the peak amplitude in the longitudinal polarization.

We have seen that all possible transverse polarization components producing a certain longitudinal component differ by a solenoidal field E sol ⊥ = ∇ × [A(r ⊥ )e z ], and the function A(r ⊥ ) represents the degrees of freedom one has when shaping the E f ⊥ . As our examples show, it is possible to control the topological structure of longitudinal and transverse electric field components simultaneously. Tightly focused beams containing vortex lines play a role in inscribing vortex lines with specific topology into Bose-Einstein condensates using two-photon Rabitransitions [28,[START_REF] Ruostekoski | [END_REF]. The demonstrated knotted or linked longitudinal vortex lines have an extent of roughly 1 µm 3 and thus match the typical size of a Bose-Einstein Condensate. Being able to exploit the unique features of structured light in all three vector components of the electric field opens new avenues in controlling the interaction of light with matter.

An important practical issue is to actually experimentally detect such a small structure in the longitudinal polarization component. Probing of the longitudinal field using molecules was achieved experimentally roughly 15 years ago [31] and continues to be of interest for light-matter interactions [32]. We propose using a tomographic method using a thermal Rubidium vapour cell that is very thin compared to the wavelength [33] to experimentally access the longitudinal polarization component. Using an additional strong static magnetic field parallel to the optical axis and tuning the light field to reso-nantly drive a π-transition allows the selective coupling of the longitudinal polarization only. To separate the π-transition from the σ ± -transitions beyond Doppler broadening (roughly 0.5GHz at 100 • C) we need a sufficiently large magnetic field (roughly B ∼ 1T ). For such large magnetic fields, isolated pure π-transitions exist e.g. from |5S 1/2 m j m I = |5S 1/2 ± 1/2 ± 3/2 to |5P 3/2 ± 1/2 ± 3/2 . This method of light-matter coupling can however be extended to more general settings, where the angle of the magnetic field can be tuned and thus different components of vectorial topological light can superposed and inscribed into matter.

In conclusion, we have presented a simple algorithm to realize an arbitrary (sufficiently wellbehaved) field in the focal plane in the longitudinal polarization component, and elaborated on how to realize the transverse components for it. We have highlighted the importance of the occurrence of evanescent waves and discussed the important degrees of freedom in the choice of the transverse polarization components. Using this method has the potential to broaden the range of possible vectorial structured light fields extensively and lead to a range of applications in various fields in physics, including nonlinear optics and Bose-Einstein Condensates.

This work is funded by the Leverhulme Trust Research Programme Grant RP2013-K-009, SPOCK: Scientific Properties Of Complex Knots. S.S. acknowledges support by the Qatar National Research Fund through the National Priorities Research Program (Grant No. NPRP 8-246-1-060).

FIG. 1 .

 1 FIG. 1. The profiles f Hopf and f Trefoil of Eqs. (4) and (5) (a,e) at narrow widths contain evanescent waves. This is demonstrated in (c,g), where the profiles are shown in the transverse Fourier domain together with a circle of radius k0. A spectral attenuation (d,h) according to Eq. (3) removes the evanescent amplitudes as well as amplitudes close to k ⊥ = 0 (see text for details), and alters E f z in the focal plane significantly (b,f).

FIG. 2 .

 2 FIG. 2. Propagation of the spectrally attenuated field shown in Figs. 1(b,f) gives rise to the vortex lines (black) in the forms of a Hopf link (a) and a trefoil (b). A phase slice is shown in the xy plane at z 0. As comparison, idealised Hopf link and trefoil are shown as insets.

  FIG. 3.Amplitude and phase for the linearly polarized transverse components Eqs. (6) and (7) that give rise to a longitudinal component forming the Hopf link Fig. 1(b), Fig. 2(a) is shown in (a-d) and the trefoil Fig. 1(f), Fig. 2(b) is shown in (e-h). The colormaps for each figure are on the right of the two rows of plots.

  FIG. 3.Amplitude and phase for the linearly polarized transverse components Eqs. (6) and (7) that give rise to a longitudinal component forming the Hopf link Fig. 1(b), Fig. 2(a) is shown in (a-d) and the trefoil Fig. 1(f), Fig. 2(b) is shown in (e-h). The colormaps for each figure are on the right of the two rows of plots.

FIG. 4 .

 4 FIG. 4. Irrotational transverse amplitude and phase profiles producing a Hopf link (a-d) and a trefoil (e-h) in the longitudinal polarization component.