Frédéric Dufaux
email: frederic.dufaux@l2s.centralesupelec.fr

Grand Challenges in Image Processing

Introduction

The field of image processing has been the subject of intensive research and development activities for several decades. This broad area encompasses topics such as image/video processing, image/video analysis, image/video communications, image/video sensing, modeling and representation, computational imaging, electronic imaging, information forensics and security, 3D imaging, medical imaging, and machine learning applied to these respective topics. Hereafter, we will consider both image and video content (i.e. sequence of images), and more generally all forms of visual information.

Rapid technological advances, especially in terms of computing power and network transmission bandwidth, have resulted in many remarkable and successful applications. Nowadays, images are ubiquitous in our daily life. Entertainment is one class of applications that has greatly benefited, including digital TV (e.g. broadcast, cable, and satellite TV), Internet video streaming, digital cinema, and video games. Beyond entertainment, imaging technologies are central in many other applications, including digital photography, video conferencing, video monitoring and surveillance, satellite imaging, but also in more distant domains such as healthcare and medicine, distance learning, digital archiving, cultural heritage or the automotive industry.

In this paper, we highlight a few research grand challenges for future imaging and video systems, in order to achieve breakthroughs to meet the growing expectations of end users. Given the vastness of the field, this list is by no means exhaustive.

A brief historical perspective

We first briefly discuss a few key milestones in the field of image processing. Key inventions in the development of photography and motion pictures can be traced to the 19 th century. The earliest surviving photograph of a real-world scene was made by Nicéphore Niépce in 1827 [START_REF] Hirsch | Seizing the light: a history of photography[END_REF]. The Lumière brothers made the first cinematographic film in 1895, with a public screening the same year [Lumière, 1996]. After decades of remarkable developments, the second half of the 20 th century saw the emergence of new technologies launching the digital revolution. While the first prototype digital camera using a Charge-Coupled Device (CCD) was demonstrated in 1975, the first commercial consumer digital cameras started appearing in the early 1990s. These digital cameras quickly surpassed cameras using films and the digital revolution in the field of imaging was underway. As a key consequence, the digital process enabled computational imaging, in other words the use of sophisticated processing algorithms in order to produce high quality images.

In 1992, the Joint Photographic Experts Group (JPEG) released the JPEG standard for still image coding [START_REF] Wallace | The JPEG still picture compression standard[END_REF]. In parallel, in 1993, the Moving Picture Experts Group (MPEG) published its first standard for coding of moving pictures and associated audio, MPEG-1 [ [START_REF] Gall | MPEG: a video compression standard for multimedia applications[END_REF], and a few years later MPEG-2 [START_REF] Haskell | Digital video: an introduction to MPEG-2[END_REF]. By guaranteeing interoperability, these standards have been essential in many successful applications and services, for both the consumer and business markets. In particular, it is remarkable that, almost 30 years later, JPEG remains the dominant format for still images and photographs.

In the late 2000s and early 2010s, we could observe a paradigm shift with the appearance of smartphones integrating a camera. Thanks to advances in computational photography, these new smartphones soon became capable of rivaling the quality of consumer digital cameras at the time. Moreover, these smartphones were also capable of acquiring video sequences. Almost concurrently, another key evolution was the development of high bandwidth networks. In particular, the launch of 4G wireless services circa 2010 enabled users to quickly and efficiently exchange multimedia content. From this point, most of us are carrying a camera, anywhere and anytime, allowing to capture images and videos at will and to seamlessly exchange them with our contacts.

As a direct consequence of the above developments, we are currently observing a boom in the usage of multimedia content. It is estimated that today 3.2 billion images are shared each day on social media platforms, and 300 hours of video are uploaded every minute on YouTube1 . In a 2019 report, Cisco estimated that video content represented 75% of all Internet traffic in 2017, and this share is forecasted to grow to 82% in 2022 [Cisco, 2019]. While Internet video streaming and Over-The-Top (OTT) media services account for a significant bulk of this traffic, other applications are also expected to see significant increases, including video surveillance and Virtual Reality (VR) / Augmented Reality (AR).

Hyper-realistic and immersive imaging

A major direction and key driver to research and development activities over the years has been the objective to deliver an ever-improving image quality and user experience.

For instance, in the realm of video, we have observed constantly increasing spatial and temporal resolutions, with the emergence nowadays of Ultra High Definition (UHD). Another aim has been to provide a sense of the depth in the scene. For this purpose, various 3D video representations have been explored, including stereoscopic 3D and multi-view [START_REF] Dufaux | Emerging Technologies for 3D Video: Creation, Coding, Transmission and Rendering[END_REF].

In this context, the ultimate goal is to be able to faithfully represent the physical world and to deliver an immersive and perceptually hyperrealist experience. For this purpose, we discuss hereafter some emerging innovations. These developments are also very relevant in VR and AR applications [START_REF] Slater | Grand challenges in virtual environments[END_REF]. Finally, while this paper is only focusing on the visual information processing aspects, it is obvious that emerging display technologies [START_REF] Masia | A survey on computational displays: Pushing the boundaries of optics, computation, and perception[END_REF] and audio also plays key roles in many application scenarios.

Light fields, point clouds, volumetric imaging

In order to wholly represent a scene, the light information coming from all the directions has to be represented. For this purpose, the 7D plenoptic function is a key concept [START_REF] Adelson | The plenoptic function and the elements of early vision[END_REF], although it is unmanageable in practice.

By introducing additional constraints, the light field representation collects radiance from rays in all directions. Therefore, it contains a much richer information, when compared to traditional 2D imaging that captures a 2D projection of the light in the scene integrating the angular domain. For instance, this allows post-capture processing such as refocusing and changing the viewpoint. However, it also entails several technical challenges, in terms of acquisition and calibration, as well as computational image processing steps including depth estimation, super-resolution, compression and image synthesis [START_REF] Ihrke | Principles of Light Field Imaging: Briefly revisiting 25 years of research[END_REF][START_REF] Wu | Light Field Image Processing: An Overview[END_REF]. The resolution trade-off between spatial and angular resolutions is a fundamental issue. With a significant fraction of the earlier work focusing on static light fields, it is also expected that dynamic light field videos will stimulate more interest in the future. In particular, dense multi-camera arrays are becoming more tractable. Finally, the development of efficient light field compression and streaming techniques is a key enabler in many applications [START_REF] Conti | Dense Light Field Coding: A Survey[END_REF].

Another promising direction is to consider a point cloud representation. A point cloud is a set of points in the 3D space represented by their spatial coordinates and additional attributes, including color pixel values, normals, or reflectance. They are often very large, easily ranging in the millions of points, and are typically sparse. One major distinguishing feature of point clouds is that, unlike images, they do not have a regular structure, calling for new algorithms. To remove the noise often present in acquired data, while preserving the intrinsic characteristics, effective 3D point cloud filtering approaches are needed [Han et al., 2017]. It is also important to develop efficient techniques for Point Cloud Compression (PCC). For this purpose, MPEG is developing two standards: Geometry-based PCC (G-PCC) and Video-based PCC (V-PCC) [START_REF] Graziosi | An overview of ongoing point cloud compression standardization activities: Video-based (V-PCC) and geometrybased (G-PCC)[END_REF]. G-PCC considers the point cloud in its native form and compress it using 3D data structures such as octrees. Conversely, V-PCC projects the point cloud onto 2D planes and then applies existing video coding schemes. More recently, deep learning-based approaches for PCC have been shown to be effective [START_REF] Guarda | Adaptive Deep Learning-based Point Cloud Geometry Coding[END_REF]. Another challenge is to develop generic and robust solutions able to handle potentially widely varying characteristics of point clouds, e.g. in terms of size and non-uniform density. Efficient solutions for dynamic point clouds are also needed. Finally, while many techniques focus on the geometric information or the attributes independently, it is paramount to process them jointly.

High Dynamic Range and Wide Color Gamut

The human visual system is able to perceive, using various adaptation mechanisms, a broad range of luminous intensities, from very bright to very dark, as experienced every day in the real world. Nonetheless, current imaging technologies are still limited in terms of capturing or rendering such a wide range of conditions. High Dynamic Range (HDR) imaging aims at addressing this issue. Wide Color Gamut (WCG) is also often associated with HDR in order to provide a wider colorimetry.

HDR has reached some levels of maturity in the context of photography. However, extending HDR to video sequences raises scientific challenges in order to provide high quality and cost-effective solutions, impacting the whole imaging processing pipeline, including content acquisition, tone reproduction, color management, coding, and display [START_REF] Dufaux | High Dynamic Range Video -From Acquisition, to Display and Applications[END_REF][START_REF] Chalmers | HDR video past, present and future: A perspective[END_REF]. Backward compatibility with legacy content and traditional systems is another issue. Despite recent progress, the potential of HDR has not been fully exploited yet.

Coding and transmission

Three decades of standardization activities have continuously improved the hybrid video coding scheme based on the principles of transform coding and predictive coding. The Versatile Video Coding (VVC) standard has been finalized in 2020 [START_REF] Bross | Developments in international video coding standardization after AVC, with an overview of Versatile Video Coding (VVC)[END_REF], achieving approximately 50% bit rate reduction for the same subjective quality when compared to its predecessor, High Efficiency Video Coding (HEVC). While substantially outperforming VVC in the short term may be difficult, one encouraging direction is to rely on improved perceptual models to further optimize compression in terms of visual quality. Another direction, which has already shown promising results, is to apply deep learning-based approaches [START_REF] Ding | Advances In Video Compression System Using Deep Neural Network: A Review And Case Studies[END_REF]. Here, one key issue is the ability to generalize these deep models to a wide diversity of video content. The second key issue is the implementation complexity, both in terms of computation and memory requirements, which is a significant obstacle to a widespread deployment. Besides, the emergence of new video formats targeting immersive communications is also calling for new coding schemes [START_REF] Wien | Standardization Status of Immersive Video Coding[END_REF].

Considering that in many application scenarios, videos are processed by intelligent analytic algorithms rather than viewed by users, another interesting track is the development of video coding for machines [START_REF] Duan | Video coding for machines: A paradigm of collaborative compression and intelligent analytics[END_REF]. In this context, the compression is optimized taking into account the performance of video analysis tasks.

The push towards hyper-realistic and immersive visual communications entails most often an increasing raw data rate. Despite improved compression schemes, more transmission bandwidth is needed. Moreover, some emerging applications, such as VR/AR, autonomous driving, and Industry 4.0, bring a strong requirement for low latency transmission, with implications on both the imaging processing pipeline and the transmission channel. In this context, the emergence of 5G wireless networks will positively contribute to the deployment of new multimedia applications, and the development of future wireless communication technologies points towards promising advances [START_REF] Da Costa | Grand Challenges in Wireless Communications[END_REF].

Human perception and visual quality assessment

It is important to develop effective models of human perception. On the one hand, it can contribute to the development of perceptually inspired algorithms. On the other hand, perceptual quality assessment methods are needed in order to optimize and validate new imaging solutions.

The notion of Quality of Experience (QoE) relates to the degree of delight or annoyance of the user of an application or service [START_REF] Callet | Qualinet white paper on definitions of quality of experience. European network on quality of experience in multimedia systems and services[END_REF]. QoE is strongly linked to subjective and objective quality assessment methods. Many years of research have resulted in the successful development of perceptual visual quality metrics based on models of human perception [START_REF] Lin | Perceptual visual quality metrics: A survey[END_REF][START_REF] Bovik | Automatic prediction of perceptual image and video quality[END_REF]. More recently, deep learning-based approaches have also been successfully applied to this problem [START_REF] Bosse | Deep neural networks for noreference and full-reference image quality assessment[END_REF]. While these perceptual quality metrics have achieved good performances, several significant challenges remain. First, when applied to video sequences, most current perceptual metrics are applied on individual images, neglecting temporal modeling. Second, whereas color is a key attribute, there are currently no widely accepted perceptual quality metrics explicitly considering color. Finally, new modalities, such as 360° videos, light fields, point clouds, and HDR, require new approaches.

Another closely related topic is image aesthetic assessment [START_REF] Deng | Image aesthetic assessment: An experimental survey[END_REF]. The aesthetic quality of an image is affected by numerous factors, such as lighting, color, contrast, and composition. It is useful in different application scenarios such as image retrieval and ranking, recommendation, and photos enhancement. While earlier attempts have used handcrafted features, most recent techniques to predict aesthetic quality are data driven and based on deep learning approaches, leveraging the availability of large annotated datasets for training [START_REF] Murray | AVA: A large-scale database for aesthetic visual analysis[END_REF]. One key challenge is the inherently subjective nature of aesthetics assessment, resulting in ambiguity in the ground-truth labels. Another important issue is to explain the behavior of deep aesthetic prediction models.

Analysis, interpretation and understanding

Another major research direction has been the objective to efficiently analyze, interpret and understand visual data. This goal is challenging, due to the high diversity and complexity of visual data. This has led to many research activities, involving both low-level and high-level analysis, addressing topics such as image classification and segmentation, optical flow, image indexing and retrieval, object detection and tracking, and scene interpretation and understanding. Hereafter, we discuss some trends and challenges.

Keypoints detection and local descriptors

Local imaging matching has been the cornerstone of many analysis tasks. It involves the detection of keypoints, i.e. salient visual points that can be robustly and repeatedly detected, and descriptors, i.e. a compact signature locally describing the visual features at each keypoint. It allows to subsequently compute pairwise matching between the features to reveal local correspondences. In this context, several frameworks have been proposed, including Scale Invariant Feature Transform (SIFT) [START_REF] Lowe | Distinctive image features from scale-invariant keypoints[END_REF] and Speeded Up Robust Features (SURF) [START_REF] Bay | Speeded-up robust features (SURF)[END_REF], and later binary variants including Binary Robust Independent Elementary Feature (BRIEF) [START_REF] Calonder | Brief: Binary robust independent elementary features[END_REF], Oriented FAST and Rotated BRIEF (ORB) [START_REF] Rublee | ORB: An efficient alternative to SIFT or SURF[END_REF] and Binary Robust Invariant Scalable Keypoints (BRISK) [START_REF] Leutenegger | BRISK: Binary robust invariant scalable keypoints[END_REF]. Although these approaches exhibit scale and rotation invariance, they are less suited to deal with large 3D distortions such as perspective deformations, out-of-plane rotations, and significant viewpoint changes. Besides, they tend to fail under significantly varying and challenging illumination conditions. These traditional approaches based on handcrafted features have been successfully applied to problems such as image and video retrieval, object detection, visual Simultaneous Localization And Mapping (SLAM), and visual odometry. Besides, the emergence of new imaging modalities as introduced above can also be beneficial for image analysis tasks, including light fields [START_REF] Galdi | Light fields for face analysis[END_REF], point clouds [START_REF] Guo | A review of algorithms for filtering the 3D point cloud[END_REF], and HDR [START_REF] Rana | Learning-based tone mapping operator for efficient image matching[END_REF]. However, when applied to highdimensional visual data for semantic analysis and understanding, these approaches based on handcrafted features have been supplanted in recent years by approaches based on deep learning.

Deep learning-based methods

Data-driven deep learning-based approaches [START_REF] Lecun | Deep learning[END_REF], and in particular the Convolutional Neural Network (CNN) architecture, represent nowadays the state-of-the-art in terms of performances for complex pattern recognition tasks in scene analysis and understanding. By combining multiple processing layers, deep models are able to learn data representations with different levels of abstraction.

Supervised learning is the most common form of deep learning. It requires a large and fully labeled training dataset, a typically time-consuming and expensive process needed whenever tackling a new application scenario. Moreover, in some specialized domains, e.g. medical data, it can be very difficult to obtain annotations. To alleviate this major burden, methods such as transfer learning and weakly supervised learning have been proposed.

In another direction, deep models have been shown to be vulnerable to adversarial attacks [START_REF] Akhtar | Threat of adversarial attacks on deep learning in computer vision: A survey[END_REF]. Those attacks consist in introducing subtle perturbations to the input, such that the model predicts an incorrect output. For instance, in the case of images, imperceptible pixel differences are able to fool deep learning models. Such adversarial attacks are definitively an important obstacle to the successful deployment of deep learning, especially in applications where safety and security are critical. While some early solutions have been proposed, a significant challenge is to develop effective defense mechanisms against those attacks.

Finally, another challenge is to enable low complexity and efficient implementations. This is especially important for mobile or embedded applications. For this purpose, further interactions between signal processing and machine learning can potentially bring additional benefits. For instance, one direction is to compress deep neural networks in order to enable their more efficient handling. Moreover, by combining traditional processing techniques with deep learning models, it is possible to develop low complexity solutions while preserving high performance.

Explainability in deep learning

While data-driven deep learning models often achieve impressive performances on many visual analysis tasks, their black-box nature often makes it inherently very difficult to understand how they reach a predicted output and how it relates to particular characteristics of the input data. However, this is a major impediment in many decision-critical application scenarios. Moreover, it is important not only to have confidence in the proposed solution, but also to gain further insights from it. Based on these considerations, some deep learning systems aim at promoting explainability [START_REF] Adadi | Peeking inside the black-box: a survey on explainable artificial intelligence (XAI)[END_REF][START_REF] Xie | Explainable deep learning: A field guide for the uninitiated[END_REF]. This can be achieved by exhibiting traits related to confidence, trust, safety, and ethics. However, explainable deep learning is still in its early phase. More developments are needed, in particular to develop a systematic theory of model explanation. Important aspects include the need to understand and quantify risk, to comprehend how the model makes predictions for transparency and trustworthiness, and to quantify the uncertainty in the model prediction. This challenge is key in order to deploy and use deep learning-based solutions in an accountable way, for instance in application domains such as healthcare or autonomous driving.

Self-supervised learning

Self-supervised learning refers to methods that learn general visual features from large-scale unlabeled data, without the need for manual annotations. Self-supervised learning is therefore very appealing, as it allows exploiting the vast amount of unlabeled images and videos available. Moreover, it is widely believed that it is closer to how humans actually learn. One common approach is to use the data to provide the supervision, leveraging its structure. More generally, a pretext task can be defined, e.g. image inpainting, colorizing grayscale images, predicting future frames in videos, by withholding some parts of the data and by training the neural network to predict it [START_REF] Jing | Self-supervised visual feature learning with deep neural networks: A survey[END_REF]. By learning an objective function corresponding to the pretext task, the network is forced to learn relevant visual features in order to solve the problem. Self-supervised learning has also been successfully applied to autonomous vehicles perception. More specifically, the complementarity between analytical and learning methods can be exploited to address various autonomous driving perception tasks, without the prerequisite of an annotated data set [START_REF] Chiaroni | Self-supervised learning for autonomous vehicles perception: A conciliation between analytical and learning methods[END_REF].

While good performances have already been obtained using self-supervised learning, further work is still needed. A few promising directions are outlined hereafter. Combining self-supervised learning with other learning methods is a first interesting path. For instance, semi-supervised learning [Van Engelen and Hoos, 2020] and few-short learning [START_REF] Fei-Fei | One-shot learning of object categories[END_REF] methods have been proposed for scenarios where limited labeled data is available. The performance of these methods can potentially be boosted by incorporating a self-supervised pre-training. The pretext task can also serve to add regularization. Another interesting trend in self-supervised learning is to train neural networks with synthetic data. The challenge here is to bridge the domain gap between the synthetic and real data. Finally, another compelling direction is to exploit data from different modalities. A simple example is to consider both the video and audio signals in a video sequence. In another example in the context of autonomous driving, vehicles are typically equipped with multiple sensors, including cameras, LIght Detection And Ranging (LIDAR), Global Positioning System (GPS), and Inertial Measurement Units (IMU). In such cases, it is easy to acquire large unlabeled multimodal datasets, where the different modalities can be effectively exploited in self-supervised learning methods.

Reproducible research and large public datasets

The reproducible research initiative is another way to further ensure high-quality research for the benefit of our community [START_REF] Vandewalle | Reproducible research in signal processing[END_REF]. Reproducibility, referring to the ability by someone else working independently to accurately reproduce the results of an experiment, is a key principle of the scientific method. In the context of image and video processing, it is usually not sufficient to provide a detailed description of the proposed algorithm. Most often, it is essential to also provide access to the code and data. This is even more imperative in the case of deep learning-based models.

In parallel, the availability of large public datasets is also highly desirable in order to support research activities. This is especially critical for new emerging modalities or specific application scenarios, where it is difficult to get access to relevant data. Moreover, with the emergence of deep learning, large datasets, along with labels, are often needed for training, which can be another burden.

Conclusion and perspectives

The field of image processing is very broad and rich, with many successful applications in both the consumer and business markets. However, many technical challenges remain in order to further push the limits in imaging technologies. Two main trends are on the one hand to always improve the quality and realism of image and video content, and on the other hand to be able to effectively interpret and understand this vast and complex amount of visual data. However, the list is certainly not exhaustive and there are many other interesting problems, e.g. related to computational imaging, information security and forensics, or medical imaging. Key innovations will be found at the crossroad of image processing, optics, psychophysics, communication, computer vision, artificial intelligence, and computer graphics. Multi-disciplinary collaborations are therefore critical moving forward, involving actors from both academia and the industry, in order to drive these breakthroughs.

The "Image Processing" section of Frontier in Signal Processing aims at giving to the research community a forum to exchange, discuss and improve new ideas, with the goal to contribute to the further advancement of the field of image processing and to bring exciting innovations in the foreseeable future.

https://www.brandwatch.com/blog/amazing-social-media-statistics-and-facts/ (accessed on Feb.

23, 2021)


Paolo Piersanti

ON THE IMPROVED INTERIOR REGULARITY OF THE SOLUTION OF A SECOND ORDER ELLIPTIC BOUNDARY VALUE PROBLEM MODELLING THE DISPLACEMENT OF A LINEARLY ELASTIC ELLIPTIC MEMBRANE SHELL SUBJECT TO AN OBSTACLE

Keywords: 2020 Mathematics Subject Classification. Primary: 74B05, 74G40; Secondary: 39A06 Elliptic membrane shell, Obstacle problems, Elliptic variational inequalities, Improved regularity, Second order problems

In this paper we show that the solution of an obstacle problem for linearly elastic elliptic membrane shells enjoys higher differentiability properties in the interior of the domain where it is defined.

Introduction.

In this paper we study the regularity of the solution of an obstacle problem for linearly elastic elliptic membrane shells, which was obtained as a result of a rigorous asymptotic analysis in the papers [START_REF] Ciarlet | Un problème de confinement pour une coque membranaire linéairement élastique de type elliptique[END_REF][START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF]. Since the solution of this second order problem is uniquely determined, the problem in object is formulated as a set of variational inequalities posed over a non-empty, closed, and convex subset of a Sobolev space.

The study of the augmentation of regularity of solutions for boundary value problems modelled via elliptic equations began between the end of the Fifties and the early Sixties, when Agmon, Douglis & Nirenberg published the two pioneering papers [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF] and [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF] about the regularity properties of solutions of elliptic systems near the boundary of the integration domain.

The augmentation of regularity for solutions of variational inequalities was first addressed by Frehse in the early Seventies [START_REF] Frehse | Zum Differenzierbarkeitsproblem bei Variationsungleichungen höherer Ordnung[END_REF][START_REF] Frehse | On the regularity of the solution of the biharmonic variational inequality[END_REF]. In the late Seventies and early Eighties, Caffarelli and his collaborators published the two papers [START_REF] Caffarelli | The obstacle problem for the biharmonic operator[END_REF][START_REF] Caffarelli | The two-obstacle problem for the biharmonic operator[END_REF], where they proved that the solution of an obstacle problem for the biharmonic operator (cf., e.g., Section 6.7 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]) could not be too regular. It was recently established in [START_REF] Piersanti | On the improved interior regularity of the solution of a fourth order elliptic problem modelling the displacement of a linearly elastic shallow shell subject to an obstacle[END_REF] that the solution of an obstacle problem for linearly elastic shallow shells enjoys higher regularity properties in the interior of the domain where it is defined.

To our best knowledge, there is no record in the literature treating the augmentation of regularity of the solution of second order variational inequalities in the case where one such solution is a vector field and the constraint defining the nonempty, closed, and convex subset of the Sobolev space where the solution is sought is expressed in terms of all of the three components of the displacement vector field. The purpose of this paper is exactly to remedy this situation.

This paper is divided into five sections (including this one). In section 2 we present some background and notation. In section 3 we recall the formulation and the properties of a three-dimensional obstacle problem for a "general" linearly elastic shell. It is worth mentioning that this three-dimensional problem is the starting point for deriving the variational formulation of the two-dimensional problem, whose solution regularity is the object of interest of this paper. In section 4 we rigorously state the two-dimensional problem we are interested in studying. Finally, in section 5, we show that the solution of the problem announced in section 4 enjoys higher differentiability properties in the interior of the domain where it is defined.

We will see that, differently from [START_REF] Frehse | Zum Differenzierbarkeitsproblem bei Variationsungleichungen höherer Ordnung[END_REF] and [START_REF] Piersanti | On the improved interior regularity of the solution of a fourth order elliptic problem modelling the displacement of a linearly elastic shallow shell subject to an obstacle[END_REF], the fact that the constraint is defined in terms of all of the three components of the displacement vector field introduces new challenges for defining a suitable test function; a test function that will be used to recover uniform estimates for the first order finite difference quotient. We will see that the choice of the test function we will use to recover the desired conclusion is justified by the density property established in [START_REF] Ciarlet | Un problème de confinement pour une coque membranaire linéairement élastique de type elliptique[END_REF][START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF] and also used in [START_REF] Ciarlet | A confinement problem for a linearly elastic Koiter's shell[END_REF][START_REF] Ciarlet | Obstacle problems for Koiter's shells[END_REF] in the context of a the justification of Koiter's model for a linearly elastic elliptic membrane shell subject to an obstacle.

2. Background and notation. For a complete overview about the classical notions of differential geometry used in this paper, see, e.g. [START_REF] Ciarlet | Mathematical Elasticity[END_REF] or [START_REF] Ciarlet | An Introduction to Differential Geometry with Applications to Elasticity[END_REF].

Greek indices, except ε, take their values in the set {1, 2}, while Latin indices, except when they are used for ordering sequences, take their values in the set {1, 2, 3}, and the summation convention with respect to repeated indices is used jointly with these two rules. As a model of the three-dimensional "physical" space R 3 , we take a real three-dimensional affine Euclidean space, i.e., a set in which a point O ∈ R 3 has been chosen as the origin and with which a real three-dimensional Euclidean space, denoted E 3 , is associated. We equip E 3 with an orthonormal basis consisting of three vectors e i , with components e i j = δ i j . The definition of R 3 as an affine Euclidean space means that with any point x ∈ R 3 is associated an uniquely determined vector Ox ∈ E 3 . The origin O ∈ R 3 and the orthonormal vectors e i ∈ E 3 together constitute a Cartesian frame in R 3 and the three components x i of the vector Ox over the basis formed by e i are called the Cartesian coordinates of x ∈ R 3 , or the Cartesian components of Ox ∈ E 3 . Once a Cartesian frame has been chosen, any point x ∈ R 3 may be thus identified with the vector Ox = x i e i ∈ E 3 . As a result, a set in R 3 can be identified with a "physical" body in the Euclidean space E 3 . The Euclidean inner product and the vector product of u, v ∈ E 3 are respectively denoted by u • v and u ∧ v; the Euclidean norm of u ∈ E 3 is denoted by |u|. The notation δ j i designates the Kronecker symbol.

Given an open subset Ω of R n , where n ≥ 1, we denote the usual Lebesgue and Sobolev spaces by L 2 (Ω), L 1 loc (Ω), H 1 (Ω), H 1 0 (Ω), H 1 loc (Ω), and the notation D(Ω) designates the space of all functions that are infinitely differentiable over Ω and have compact supports in Ω. We denote • X the norm in a normed vector space X. Spaces of vector-valued functions are written in boldface. The Euclidean norm of any point x ∈ Ω is denoted by |x|.

The boundary Γ of an open subset Ω in R n is said to be Lipschitz-continuous if the following conditions are satisfied (cf., e.g., Section 1.18 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]): Given an integer s ≥ 1, there exist constants α 1 > 0 and L > 0, a finite number of local coordinate systems, with coordinates

φ r = (φ r 1 , . . . , φ r n-1 ) ∈ R n-1 and φ r = φ r n , 1 ≤ r ≤ s, sets ωr := {φ r ∈ R n-1 ; |φ r | < α 1 }, 1 ≤ r ≤ s,
and corresponding functions

θr : ωr → R, 1 ≤ r ≤ s, such that Γ = s r=1
{(φ r , φ r ); φ r ∈ ωr and φ r = θr (φ r )}, and

| θr (φ r ) -θr (υ r )| ≤ L|φ r -υ r |
, for all φ r , υ r ∈ ωr , and all 1 ≤ r ≤ s.

We observe that the second last formula takes into account overlapping local charts, while the last set of inequalities expresses the Lipschitz continuity of the mappings θr .

An open set Ω is said to be locally on the same side of its boundary Γ if, in addition, there exists a constant α 2 > 0 such that {(φ r , φ r ); φ r ∈ ωr and θr (φ r ) < φ r < θr (φ r ) + α 2 } ⊂ Ω, for all 1 ≤ r ≤ s, {(φ r , φ r ); φ r ∈ ωr and θr (φ r ) -α 2 < φ r < θr (φ r )} ⊂ R n \ Ω, for all 1 ≤ r ≤ s.

A domain in R n is a bounded and connected open subset Ω of R n , whose boundary ∂Ω is Lipschitz-continuous, the set Ω being locally on a single side of ∂Ω.

Let ω be a domain in R 2 with boundary γ := ∂ω, and let ω 1 ⊂ ω. The special notation ω 1 ⊂⊂ ω means that ω 1 ⊂ ω and dist(γ, ∂ω 1 ) > 0. Let y = (y α ) denote a generic point in ω, and let ∂ α := ∂/∂y α . A mapping θ ∈ C 1 (ω; E 3 ) is said to be an immersion if the two vectors

a α (y) := ∂ α θ(y)
are linearly independent at each point y ∈ ω. Then the set θ(ω) is a surface in E 3 , equipped with y 1 , y 2 as its curvilinear coordinates. Given any point y ∈ ω, the linear combinations of the vectors a α (y) span the tangent plane to the surface θ(ω) at the point θ(y), the unit vector

a 3 (y) := a 1 (y) ∧ a 2 (y) |a 1 (y) ∧ a 2 (y)|
is orthogonal to θ(ω) at the point θ(y), the three vectors a i (y) form the covariant basis at the point θ(y), and the three vectors a j (y) defined by the relations a j (y) • a i (y) = δ j i form the contravariant basis at θ(y); note that the vectors a β (y) also span the tangent plane to θ(ω) at θ(y) and that a 3 (y) = a 3 (y).

The first fundamental form of the surface θ(ω) is then defined by means of its covariant components a αβ := a α • a β = a βα ∈ C 0 (ω), or by means of its contravariant components

a αβ := a α • a β = a βα ∈ C 0 (ω).
Note that the symmetric matrix field (a αβ ) is then the inverse of the matrix field (a αβ ), that a β = a αβ a α and a α = a αβ a β , and that the area element along θ(ω) is given at each point θ(y), y ∈ ω, by a(y) dy, where

a := det(a αβ ) ∈ C 0 (ω).
Given an immersion θ ∈ C 2 (ω; E 3 ), the second fundamental form of the surface θ(ω) is defined by means of its covariant components

b αβ := ∂ α a β • a 3 = -a β • ∂ α a 3 = b βα ∈ C 0 (ω),
or by means of its mixed components b β α := a βσ b ασ ∈ C 0 (ω), and the Christoffel symbols associated with the immersion θ are defined by

Γ σ αβ := ∂ α a β • a σ = Γ σ βα ∈ C 0 (ω). The Gaussian curvature at each point θ(y), y ∈ ω, of the surface θ(ω) is defined by κ(y) := det(b αβ (y)) det(a αβ (y)) = det b β α (y) .
Observe that the denominator in the above relation does not vanish since θ is assumed to be an immersion. Note that the Gaussian curvature κ(y) at the point θ(y) is also equal to the product of the two principal curvatures at this point.

Given an immersion θ ∈ C 2 (ω; E 3 ) and a vector field η = (η i ) ∈ C 1 (ω; R 3 ), the vector field η := η i a i may be viewed as the displacement field of the surface θ(ω), thus defined by means of its covariant components η i over the vectors a i of the contravariant bases along the surface. If the norms η i C 1 (ω) are small enough, the mapping (θ+η i a i ) ∈ C 1 (ω; E 3 ) is also an immersion, so that the set (θ + η i a i )(ω) is again a surface in E 3 , equipped with the same curvilinear coordinates as those of the surface θ(ω) and is called the deformed surface corresponding to the displacement field η = η i a i .

It is thus possible to define the first fundamental form of the deformed surface in terms of its covariant components by

a αβ (η) :=(a α + ∂ α η) • (a β + ∂ β η) =a αβ + a α • ∂ β η + ∂ α η • a β + ∂ α η • ∂ β η.
The linear part with respect to η in the difference (a αβ (η) -a αβ )/2 is called the linearised change of metric, or strain, tensor associated with the displacement field η i a i , the covariant components of which are thus defined by

γ αβ (η) := 1 2 (a α • ∂ β η + ∂ α η • a β ) = 1 2 (∂ β η α + ∂ α η β ) -Γ σ αβ η σ -b αβ η 3 = γ βα (η).
In this paper, we shall consider a specific class of surfaces, according to the following definition: Let ω be a domain in R 2 . Then a surface θ(ω) defined by means of an immersion θ ∈ C 2 (ω; E 3 ) is said to be elliptic if its Gaussian curvature is everywhere > 0 in ω, or equivalently, if there exists a constant κ 0 such that 0 < κ 0 ≤ κ(y), for all y ∈ ω.

It turns out that, when an elliptic surface is subjected to a displacement field η i a i whose tangential covariant components η α vanish on the entire boundary of the domain ω, the following inequality holds. Note that the components of the displacement fields and linearised change of metric tensors appearing in the next theorem are no longer assumed to be continuously differentiable functions; they are instead to be understood in a generalised sense, since they now belong to ad hoc Lebesgue or Sobolev spaces.

Theorem 1. Let ω be a domain in R 2 and let an immersion θ ∈ C 3 (ω; E 3 ) be given such that the surface θ(ω) is elliptic. Define the space

V M (ω) := H 1 0 (ω) × H 1 0 (ω) × L 2 (ω). Then, there exists a constant c 0 = c 0 (ω, θ) > 0 such that α η α 2 H 1 (ω) + η 3 2 L 2 (ω) 1/2 ≤ c 0 α,β γ αβ (η) L 2 (ω) 1/2 for all η = (η i ) ∈ V M (ω).
The above inequality, which is due to [START_REF] Ciarlet | On the ellipticity of linear membrane shell equations[END_REF] and [START_REF] Ciarlet | An existence and uniqueness theorem for the two-dimensional linear membrane shell equations[END_REF] (see also Theorem 2.7-3 of [START_REF] Ciarlet | Mathematical Elasticity[END_REF]), constitutes an example of a Korn's inequality on a surface, in the sense that it provides an estimate of an appropriate norm of a displacement field defined on a surface in terms of an appropriate norm of a specific "measure of strain" (here, the linearised change of metric tensor) corresponding to the displacement field considered.

3. The three-dimensional obstacle problem for a "general" linearly elastic shell. Let ω be a domain in R 2 , let γ := ∂ω, and let γ 0 be a non-empty relatively open subset of γ. For each ε > 0, we define the sets

Ω ε = ω × ]-ε, ε[ and Γ ε 0 := γ 0 × ]-ε, ε[ , we let x ε = (x ε
i ) designate a generic point in the set Ω ε , and we let ∂ ε i := ∂/∂x ε i . Hence we also have x ε α = y α and ∂ ε α = ∂ α . Given an immersion θ ∈ C 3 (ω; E 3 ) and ε > 0, consider a shell with middle surface θ(ω) and with constant thickness 2ε. This means that the reference configuration of the shell is the set Θ(Ω ε ), where the mapping Θ : Ω ε → E 3 is defined by

Θ(x ε ) := θ(y) + x ε 3 a 3 (y) at each point x ε = (y, x ε 3 ) ∈ Ω ε .
One can then show (cf., e.g., Theorem 3.1-1 of [START_REF] Ciarlet | Mathematical Elasticity[END_REF]) that, if ε > 0 is small enough, such a mapping Θ ∈ C 2 (Ω ε ; E 3 ) is an immersion, in the sense that the three vectors

g ε i (x ε ) := ∂ ε i Θ(x ε
) are linearly independent at each point x ε ∈ Ω ε ; these vectors then constitute the covariant basis at the point Θ(x ε ), while the three vectors g j,ε (x ε ) defined by the relations g j,ε (x ε ) • g ε i (x ε ) = δ j i constitute the contravariant basis at the same point. It will be implicitly assumed in the sequel that ε > 0 is small enough so that Θ : Ω ε → E 3 is an immersion.

One then defines the metric tensor associated with the immersion Θ by means of its covariant components

g ε ij := g ε i • g ε j ∈ C 1 (Ω ε
), or by means of its contravariant components

g ij,ε := g i,ε • g j,ε ∈ C 1 (Ω ε ).
Note that the symmetric matrix field (g ij,ε ) is then the inverse of the matrix field (g ε ij ), that g j,ε = g ij,ε g ε i and g ε i = g ε ij g j,ε , and that the volume element in Θ(Ω ε ) is given at each point Θ(x ε ), x ε ∈ Ω ε , by g ε (x ε ) dx ε , where

g ε := det(g ε ij ) ∈ C 1 (Ω ε
). One also defines the Christoffel symbols associated with the immersion Θ by

Γ p,ε ij := ∂ i g ε j • g p,ε = Γ p,ε ji ∈ C 0 (Ω ε ). Note that Γ 3,ε α3 = Γ p,ε 33 = 0. Given a vector field v ε = (v ε i ) ∈ C 1 (Ω ε ; R 3
), the associated vector field ṽε := v ε i g i,ε can be viewed as a displacement field of the reference configuration Θ(Ω ε ) of the shell, thus defined by means of its covariant components v ε i over the vectors g i,ε of the contravariant bases in the reference configuration.

If the norms

v ε i C 1 (Ω ε ) are small enough, the mapping (Θ + v ε i g i,ε
) is also an immersion, so that one can also define the metric tensor of the deformed configuration

(Θ + v ε i g i,ε )(Ω ε ) by means of its covariant components g ε ij (v ε ) :=(g ε i + ∂ ε i ṽε ) • (g ε j + ∂ ε j ṽε ) =g ε ij + g ε i • ∂ j ṽε + ∂ ε i ṽε • g ε j + ∂ i ṽε • ∂ j ṽε .
The linear part with respect to ṽε in the difference (g ε ij (v ε ) -g ε ij )/2 is then called the linearised strain tensor associated with the displacement field v ε i g i,ε , the covariant components of which are thus defined by

e ε i j (v ε ) := 1 2 g ε i • ∂ j ṽε + ∂ ε i ṽε • g ε j = 1 2 (∂ ε j v ε i + ∂ ε i v ε j ) -Γ p,ε ij v ε p = e ε j i (v ε ).
The functions e ε i j (v ε ) are called the linearised strains in curvilinear coordinates associated with the displacement field v ε i g i,ε . We assume throughout this paper that, for each ε > 0, the reference configuration Θ(Ω ε ) of the shell is a natural state (i.e., stress-free) and that the material constituting the shell is homogeneous, isotropic, and linearly elastic. The behaviour of such an elastic material is thus entirely governed by its two Lamé constants λ ≥ 0 and µ > 0 (for details, see, e.g., Section 3.8 of [START_REF] Ciarlet | Mathematical Elasticity. Vol. I: Three-Dimensional Elasticity[END_REF]).

We will also assume that the shell is subjected to applied body forces whose density per unit volume is defined by means of its covariant components f i,ε ∈ L 2 (Ω ε ), and to a homogeneous boundary condition of place along the portion Γ ε 0 of its lateral face (i.e., the displacement vanishes on Γ ε 0 ). In this paper, we consider a specific obstacle problem for such a shell, in the sense that the shell is also subjected to a confinement condition, expressing that any admissible deformed configuration remains in a half-space of the form

H := {Ox ∈ E 3 ; Ox • q ≥ 0},
where q ∈ E 3 is a non-zero vector given once and for all. In other words, any admissible displacement field must satisfy

Θ(x ε ) + v ε i (x ε )g i,ε (x ε ) • q ≥ 0 for all x ε ∈ Ω ε ,
or possibly only for almost all (a.a. in what follows) x ε ∈ Ω ε when the covariant components v ε i are required to belong to the Sobolev space H 1 (Ω ε ) as in Theorem 2 below.

We will of course assume that the reference configuration satisfies the confinement condition, i.e., that

Θ(Ω ε ) ⊂ H.
It is to be emphasised that the above confinement condition considerably departs from the usual Signorini condition favoured by most authors, who usually require that only the points of the undeformed and deformed "lower face" ω × {-ε} of the reference configuration satisfy the confinement condition (see, e.g., [START_REF] Léger | Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF], [START_REF] Léger | A linearly elastic shell over an obstacle: The flexural case[END_REF], [START_REF] Mezabia | Modelling of frictionless Signorini problem for a linear elastic membrane shell, Applicable Analysis[END_REF], [START_REF] Rodríguez-Arós | Mathematical justification of the obstacle problem for elastic elliptic membrane shells[END_REF]). Clearly, the confinement condition considered in the present paper is more physically realistic, since a Signorini condition imposed only on the lower face of the reference configuration does not prevent -at least "mathematically" -other points of the deformed reference configuration to "cross" the plane {Ox ∈ E 3 ; Ox • q = 0} and then to end up on the "other side" of this plane. The mathematical models characterised by the confinement condition introduced beforehand, confinement condition which is also considered in the seminal paper [START_REF] Léger | Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF] in a different geometrical framework, do not take any traction forces into account. Indeed, there could be no traction forces applied to the portion of the three-dimensional shell boundary that engages contact with the obstacle. Besides, friction is not considered in the context of this analysis.

Unlike the classical Signorni condition, the confinement condition here considered is more suitable in the context of multi-scales multi-bodies problems like, for instance, the study of the motion of the human heart valves, conducted by Quarteroni and his associates in [START_REF] Piersanti | Modeling cardiac muscle fibers in ventricular and atrial electrophysiology simulations[END_REF][START_REF] Regazzoni | Active force generation in cardiac muscle cells: mathematical modeling and numerical simulation of the actinmyosin interaction[END_REF][START_REF] Zingaro | Hemodynamics of the heart's left atrium based on a Variational Multiscale-LES numerical method[END_REF] and the references therein.

Such a confinement condition renders the study of this problem considerably more difficult, however, as the constraint now bears on a vector field, the displacement vector field of the reference configuration, instead of on only a single component of this field.

The mathematical modelling of such an obstacle problem for a linearly elastic shell is then clear; since, apart from the confinement condition, the rest, i.e., the function space and the expression of the quadratic energy J ε , is classical (see, e.g. [START_REF] Ciarlet | Mathematical Elasticity[END_REF]). More specifically, let

A ijk ,ε := λg ij,ε g k ,ε + µ g ik,ε g j ,ε + g i ,ε g jk,ε = A jik ,ε = A k ij,ε
denote the contravariant components of the elasticity tensor of the elastic material constituting the shell. Then the unknown of the problem, which is the vector field u ε = (u ε i ) where the functions u ε i : Ω ε → R are the three covariant components of the unknown "three-dimensional" displacement vector field u ε i g i,ε of the reference configuration of the shell, should minimise the energy

J ε : H 1 (Ω ε ) → R defined by J ε (v ε ) := 1 2 Ω ε A ijk ,ε e ε k (v ε )e ε i j (v ε ) √ g ε dx ε - Ω ε f i,ε v ε i √ g ε dx ε for each v ε = (v ε i ) ∈ H 1 (Ω ε
) over the set of admissible displacements defined by:

U (Ω ε ) :={v ε = (v ε i ) ∈ H 1 (Ω ε ); v ε = 0 on Γ ε 0 , (Θ(x ε ) + v ε i (x ε )g i,ε (x ε )) • q ≥ 0 for a.a. x ε ∈ Ω ε }.
The solution to this minimisation problem exists and is unique, and it can be also characterised as the unique solution of a set of appropriate variational inequalities (cf., Theorem 2.1 of [START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF]).

Theorem 2. The quadratic minimisation problem: Find a vector field

u ε ∈ U (Ω ε ) such that J ε (u ε ) = inf v ε ∈U (Ω ε ) J ε (v ε )
has one and only one solution. Besides, the vector field u ε is also the unique solution of the variational problem P(Ω ε ): Find u ε ∈ U (Ω ε ) that satisfies the following variational inequalities:

Ω ε A ijk ,ε e ε k (u ε ) e ε i j (v ε ) -e ε i j (u ε ) √ g ε dx ε ≥ Ω ε f i,ε (v ε i -u ε i ) √ g ε dx ε for all v ε = (v ε i ) ∈ U (Ω ε ). Since θ(ω) ⊂ Θ(Ω ε
), it evidently follows that θ(y) • q ≥ 0 for all y ∈ ω. But in fact, a stronger property holds (cf., Lemma 2.1 of [START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF]):

Lemma 1. Let ω be a domain in R 2 , let θ ∈ C 1 (ω; E 3
) be an immersion, let q ∈ E 3 be a non-zero vector, and let ε > 0. Then the inclusion

Θ(Ω ε ) ⊂ H = {x ∈ E 3 ; Ox • q ≥ 0} implies that min y∈ω (θ(y) • q) > 0.
Clearly, the assumed inclusion Θ(Ω ε ) ⊂ H implies that min y∈ω (θ(y) • q) depends on ε > 0. In order to identify the two-dimensional problem that constitutes the object of interest of this paper via a rigorous asymptotic analysis, we need to assume that min y∈ω (θ(y) • q) is strictly positive and independent of ε > 0.

4. The scaled three-dimensional problem for a family of linearly elastic elliptic membrane shells. In section 3, we considered an obstacle problem for "general" linearly elastic shells. From now on, we will restrict ourselves to a specific class of shells, according to the following definition (proposed in [START_REF] Ciarlet | Asymptotic analysis of linearly elastic shells. I. Justification of membrane shell equations[END_REF]; see also [START_REF] Ciarlet | Mathematical Elasticity[END_REF]).

Consider a linearly elastic shell, subjected to the various assumptions set forth in section 3. Such a shell is said to be a linearly elastic elliptic membrane shell (from now on simply membrane shell ) if the following two additional assumptions are satisfied: first, γ 0 = γ, i.e., the homogeneous boundary condition of place is imposed over the entire lateral face γ × ]-ε, ε[ of the shell, and second, its middle surface θ(ω) is elliptic, according to the definition given in section 2.

In this paper, we consider the obstacle problem (as defined in section 3) for a family of membrane shells, all sharing the same middle surface and whose thickness 2ε > 0 is considered as a "small" parameter approaching zero. In order to conduct an asymptotic analysis on the three-dimensional model as the thickness ε → 0, we resorted to a (by now standard) methodology first proposed in [START_REF] Ciarlet | A justification of the two-dimensional linear plate model[END_REF]: To begin with, we "scale" each problem P(Ω ε ), ε > 0, over a fixed domain Ω, using appropriate scalings on the unknowns and assumptions on the data.

More specifically, let Ω := ω × ]-1, 1[ , let x = (x i ) denote a generic point in the set Ω, and let

∂ i := ∂/∂x i . With each point x = (x i ) ∈ Ω, we associate the point x ε = (x ε i ) defined by x ε α := x α = y α and x ε 3 := εx 3 , so that ∂ ε α = ∂ α and ∂ ε 3 = ε -1 ∂ 3 .
To the unknown u ε = (u ε i ) and to the vector fields v ε = (v ε i ) appearing in the formulation of the problem P(Ω ε ) corresponding to a membrane shell, we then associate the scaled unknown u(ε) = (u i (ε)) and the scaled vector fields v = (v i ) by letting

u i (ε)(x) := u ε i (x ε ) and v i (x) := v ε i (x ε ) at each x ∈ Ω.
Finally, we assume that there exist functions f i ∈ L 2 (Ω) independent on ε such that the following assumptions on the data hold

f i,ε (x ε ) = f i (x) at each x ∈ Ω.
Note that the independence on ε of the Lamé constants assumed in Section 3 in the formulation of problem P(Ω ε ) implicitly constituted another assumption on the data.

The variational problem P(ε; Ω) defined in the next theorem will constitute the point of departure of the asymptotic analysis performed in [START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF].

Theorem 3. For each ε > 0, define the set

U (ε; Ω) := {v = (v i ) ∈ H 1 (Ω); v = 0 on γ × ]-1, 1[ , θ(y) + εx 3 a 3 (y) + v i (x)g i (ε)(x) • q ≥ 0 for a.a. x = (y, x 3 ) ∈ Ω}, where g i (ε)(x) := g i,ε (x ε ) at each x ∈ Ω.
Then the scaled unknown of the variational problem P(Ω ε ) is the unique solution of the variational problem P(ε; Ω): Find u(ε) ∈ U (ε; Ω) that satisfies the following variational inequalities:

Ω A ijk (ε)e k (ε; u(ε)) e i j (ε; v) -e i j (ε; u(ε)) g(ε) dx ≥ Ω f i (v i -u i (ε)) g(ε) dx for all v ∈ U (ε; Ω),
where

g(ε)(x) := g ε (x ε ) and A ijk (ε)(x) := A ijk ,ε (x ε ) at each x ∈ Ω, e α β (ε; v) := 1 2 (∂ β v α + ∂ α v β ) -Γ k αβ (ε)v k = e β α (ε; v), e 3 α (ε; v) := 1 2 1 ε ∂ 3 v α + ∂ α v 3 -Γ σ α3 (ε)v σ = e α 3 (ε; v), e 3 3 (ε; v) := 1 ε ∂ 3 v 3 , where Γ p ij (ε)(x) := Γ p,ε ij (x ε ) at each x ∈ Ω.
The problem we are interested in is derived as a result of the rigorous asymptotic analysis conducted in Theorem 4.1 of [START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF].

Theorem 4. Let ω be a domain in R 2 , let θ ∈ C 3 (ω; E 3 ) be an immersion such that the surface θ(ω) is elliptic (cf. section 2), and let q ∈ E 3 be a non-zero vector. Define the space and sets

V M (ω) := H 1 0 (ω) × H 1 0 (ω) × L 2 (ω), U M (ω) := {η = (η i ) ∈ H 1 0 (ω) × H 1 0 (ω) × L 2 (ω); θ(y) + η i (y)a i (y) • q ≥ 0 for a.a. y ∈ ω}, ŨM (ω) := {η = (η i ) ∈ H 1 0 (ω) × H 1 0 (ω) × H 1 0 (ω); θ(y) + η i (y)a i (y) • q ≥ 0 for a.a. y ∈ ω},
and assume that the immersion θ is such that

d := min y∈ω (θ(y) • q) > 0,
is independent of ε, and assume that the following "density property" holds:

ŨM (ω) is dense in U M (ω) with respect to the norm of • H 1 (ω)×H 1 (ω)×L 2 (ω) .
Let there be given a family of membrane shells with the same middle surface θ(ω) and thickness 2ε > 0, and let

u(ε) = (u i (ε)) ∈ U (ε; Ω) := {v = (v i ) ∈ H 1 (Ω); v = 0 on γ × ]-1, 1[ , θ(y) + εx 3 a 3 (y) + v i (x)g i (ε)(x) • q ≥ 0 for a.a. x = (y, x 3 ) ∈ Ω}
denote for each ε > 0 the unique solution of the corresponding problem P(ε; Ω) introduced in Theorem 3.

Then there exist functions u α ∈ H 1 (Ω) independent of the variable x 3 and satisfying

u α = 0 on γ × ]-1, 1[ ,
and there exists a function u 3 ∈ L 2 (Ω) independent of the variable x 3 , such that

u α (ε) → u α in H 1 (Ω) and u 3 (ε) → u 3 in L 2 (Ω).

Define the average

u = (u i ) := 1 2 1 -1 u dx 3 ∈ V M (ω). Then u = ζ,
where ζ is the unique solution of the two-dimensional variational problem P M (ω):

Find ζ ∈ U M (ω) that satisfies the following variational inequalities ω a αβστ γ στ (ζ)γ αβ (η -ζ) √ a dy ≥ ω p i (η i -ζ i ) √ a dy for all η = (η i ) ∈ U M (ω),
where

a αβστ := 4λµ λ + 2µ
a αβ a στ + 2µ a ασ a βτ + a ατ a βσ and p i :=

1 -1 f i dx 3 .
Critical to establish the convergence of the family (u(ε)) ε>0 is the "density property" assumed there, which asserts that the set ŨM (ω) is dense in the set U M (ω) with respect to the norm • H 1 (ω)×H 1 (ω)×L 2 (ω) . The same "density property" is used to provide a justification, via a rigorous asymptotic analysis, of Koiter's model for membrane shells subject to an obstacle (cf. [START_REF] Ciarlet | A confinement problem for a linearly elastic Koiter's shell[END_REF], [START_REF] Ciarlet | Obstacle problems for Koiter's shells[END_REF]). We hereby recall a sufficient geometric condition ensuring the assumed "density property" (cf. Theorem 5.1 of [START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF]).

Theorem 5. Let θ ∈ C 2 (ω; E 3 ) be an immersion with the following property: There exists a non-zero vector q ∈ E 3 such that min y∈ω (θ(y) • q) > 0 and min y∈ω (a 3 (y) • q) > 0.

Define the sets

U M (ω) := {η =(η i ) ∈ H 1 0 (ω) × H 1 0 (ω) × L 2 (ω); θ(y) + η i (y)a i (y) • q ≥ 0 for a.a. y ∈ ω}, U M (ω) ∩ D(ω) := {η =(η i ) ∈ D(ω) × D(ω) × D(ω); θ(y) + η i (y)a i (y) • q ≥ 0 for a.a. y ∈ ω}. Then the set U M (ω) ∩ D(ω) is dense in the set U M (ω) with respect to the norm • H 1 (ω)×H 1 (ω)×L 2 (ω) .
Examples of elliptic membrane shells satisfying the "density property" thus include those whose middle surface is a portion of an ellipsoid that is strictly contained in one of the open half-spaces that contain two of its main axes (see, e.g., that drawn in Figure 4.1-1 in [START_REF] Ciarlet | Mathematical Elasticity[END_REF]), the boundary of the half-space coinciding with the obstacle in this case.

As a final step, we de-scale Problem P M (ω) and we get the following variational formulation (cf. Theorem 4.2 of [START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF]).

Problem P ε M (ω). Find ζ ε = (ζ ε i ) ∈ U M (ω)
satisfying the following variational inequalities:

ε ω a αβστ γ στ (ζ ε )γ αβ (η -ζ ε ) √ a dy ≥ ω p i,ε (η i -ζ ε i ) √ a dy, for all η = (η i ) ∈ U M (ω), where p i,ε := ε -ε f i,ε dx ε 3 .
By virtue of the Korn inequality recalled in Theorem 1, it results that Problem P ε M (ω) admits a unique solution. The main objective of this paper is to show that the solution of the Problem P ε M (ω), derived from the asymptotic analysis conducted in Theorem 4, is also of class H 2 loc (ω) × H 2 loc (ω) × H 1 loc (ω). We will see that a crucial role for establishing the desired higher regularity property will be played by the assumed "density property".

Augmentation of interior regularity of the solution of Problem

P ε M (ω). Let ω 0 ⊂ ω and ω 1 ⊂ ω be such that ω 1 ⊂⊂ ω 0 ⊂⊂ ω. (1)
By the definition of the symbol ⊂⊂ in (1), we obtain that the quantity

d := 1 2 min{dist(∂ω 1 , ∂ω 0 ), dist(∂ω 0 , γ)} (2)
is, actually, strictly greater than zero. Let ϕ 1 ∈ D(ω) be such that supp ϕ 1 ⊂⊂ ω 1 and 0 ≤ ϕ 1 ≤ 1.

In the rest of the paper, the immersion θ modelling the middle surface of the membrane shell under consideration is assumed to be of class C 3 (ω; E 3 ).

Denote by D ρh the first order (forward) finite difference quotient of either a function or a vector field in the canonical direction e ρ of R 2 and with increment size 0 < h < d sufficiently small. We can regard the first order (forward) finite difference quotient of a function as a linear operator defined as follows:

D ρh : L 2 (ω) → L 2 (ω 0 ).
The first order finite difference quotient of a function ξ in the canonical direction e ρ of R 2 and with increment size 0 < h < d is defined by:

D ρh ξ(y) := ξ(y + he ρ ) -ξ(y) h ,
for all (or, possibly, a.a.) y ∈ ω such that (y + he ρ ) ∈ ω.

The first order finite difference quotient of a vector field ξ = (ξ α ) in the canonical direction e ρ of R 2 and with increment size 0 < h < d is defined by

D ρh ξ(y) := ξ(y + he ρ ) -ξ(y) h ,
or, equivalently, D ρh ξ(y) = (D ρh ξ α (y)). Similarly, we can show that the first order (forward) finite difference quotient of a vector field is a linear operator from L 2 (ω) to L 2 (ω 0 ).

We define the second order finite difference quotient of a function ξ in the canonical direction e ρ of R 2 and with increment size 0 < h < d by

δ ρh ξ(y) := ξ(y + he ρ ) -2ξ(y) + ξ(y -he ρ ) h 2 ,
for all (or, possibly, a.a.) y ∈ ω such that (y ± he ρ ) ∈ ω.

The second order finite difference quotient of a vector field ξ = (ξ α ) in the canonical direction e ρ of R 2 and with increment size 0 < h < d is defined by

δ ρh ξ(y) := ξ α (y + he ρ ) -2ξ α (y) + ξ α (y -he ρ ) h 2 ,
for all (or, possibly, a.a.) y ∈ ω such that (y ± he ρ ) ∈ ω. Note in passing that the second order finite difference quotient of a function ξ can be expressed in terms of the first order finite difference quotient via the following identity:

δ ρh ξ = D -ρh D ρh ξ.

(3) Similarly, the second order finite difference quotient of a vector field ξ can be expressed in terms of the first order finite difference quotient via the following identity:

δ ρh ξ = D -ρh D ρh ξ.

(4) Before proving the main result of this section, we briefly outline the main steps we are going to implement for showing that the solution

ζ ε of Problem P ε M (ω) is also of class H 2 loc (ω) × H 2 loc (ω) × H 1 loc (ω).
The main difficulty in doing so amounts to constructing an admissible test vector field from which it is possible to infer the desired augmented regularity property.

Our strategy consists in locally perturbing the solution ζ ε by a vector field with compact support in a neighbourhood of ω, compact support which is sufficiently far away from the boundary γ.

In order to implement this local perturbation strategy, we have to first fix a neighbourhood U 1 ⊂ ω sufficiently far away from the boundary γ, and consider a smooth function ϕ whose support is compact in U 1 .

Second, we multiply the smooth function ϕ by a suitable coefficient > 0, = O(h 2 ) and a suitable second order finite difference quotient. The purpose of the positive coefficient is to ensure that the geometrical constraint appearing in the set U M (ω) is satisfied in the selected neighbourhood U 1 . The argument of the second order difference quotient is the product between the smooth function ϕ and a smooth approximation of the solution. The approximation by smooth functions requirement cannot be dropped since the transverse component of the solution, namely ζ ε 3 , is a priori only of class L 2 (ω). This regularity is in general too low for working with a second order difference quotient.

Thanks to the assumed "density property", we are in a position to consider an approximation by smooth functions of the solution ζ ε , approximations which all belong to the non-empty, closed and convex subset U M (ω).

Finally, using Theorem 1 and classical integration-by-parts formulas for finite difference quotients (cf., e.g., [START_REF] Evans | Partial Differential Equations[END_REF]), we aim to show that the finite difference quotient D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) is uniformly bounded with respect to h, where the presence of the smooth function is needed to insure the local validity of the result.

To begin with, we establish two abstract preparatory lemmas.

Lemma 2. Let ω, ω 0 and ω 1 be subsets of R 2 as in (1). For each y ∈ ω, let a i (y) denote a vector of the contravariant basis at the point θ(y) of the surface θ(ω).

Assume that ϑ ∈ C 1 (ω; E 3 ) is such that ϑ • q concave on ω 0 . Let η = (η i ) ∈ H 1 (ω) × H 1 (ω) × L 2 (ω) be such that (ϑ(y) + η i (y)a i (y)) • q ≥ 0 for a.a. y ∈ ω 0 .

Then, for each 0 < h < d and all 0 < < h 2 /2, the vector field η = (η ,i ) ∈ L 2 (ω 1 ) defined in a way such that η ,j a j := η i a i + ϕ 1 δ ρh (η i a i ), is such that (ϑ(y) + η ,i (y)a i (y)) • q ≥ 0 for a.a. y ∈ ω 1 .

Proof. For a.a. y ∈ ω 1 we have that

η ,j (y)a j (y) = η i (y)a i (y) + ϕ 1 (y)δ ρh (η i a i )(y) = η i (y)a i (y) + ϕ 1 (y) η i (y + he ρ )a i (y + he ρ ) -2η i (y)a i (y) + η i (y -he ρ )a i (y -he ρ ) h 2 = h 2 ϕ 1 (y)(η i a i )(y + he ρ ) + 1 - 2 h 2 ϕ 1 (y) (η i a i )(y) + h 2 ϕ 1 (y)(η i a i )(y -he ρ ).
By virtue of the properties of ϕ 1 and , the functions c 1 (y) = c -1 (y) := h -2 ϕ 1 (y) and c 0 (y) := 1 -2 h -2 ϕ 1 (y) are non-negative in ω. Combining these properties with the assumption that (ϑ + η i a i ) • q ≥ 0 almost everywhere in ω 0 gives:

η ,j (y)a j (y) • q ≥ -c 1 (y)ϑ(y + he ρ ) • q -c 0 (y)ϑ(y) • q -c -1 (y)ϑ(y -he ρ ) • q = -ϑ(y) + ϕ 1 (y) ϑ(y + he ρ ) -2ϑ(y) + ϑ(y -he ρ ) h 2 • q,
for a.a. y ∈ ω 1 .

The assumed concavity of ϑ • q in ω 0 gives

δ ρh ϑ(y) • q = ϑ(y + he ρ ) -2ϑ(y) + ϑ(y -he ρ ) h 2 • q < 0,
for all 0 < h < d and all y ∈ ω 1 . Recalling that supp ϕ 1 ⊂⊂ ω 1 , we derive (ϑ(y) + η ,i (y)a i (y)) • q ≥ 0 for a.a. y ∈ ω 1 , as it was to be proved.

We observe that the concavity of the function ϑ • q has to be assumed over ω 0 since the points of the form (y ± he ρ ), with y ∈ ω 1 , may lie outside ω 1 .

The next result is a direct application of Lemma 2 and of the "density property" (cf. ( 6) below).

Corollary 1. Let ω, ω 0 and ω 1 be subsets of R 2 as in [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF]. For each y ∈ ω, let a i (y) denote a vector of the contravariant basis at the point θ(y) of the surface θ(ω).

Assume that ϑ ∈ C 1 (ω;

E 3 ) is such that ϑ • q concave on ω 0 . Let η = (η i ) ∈ H 1 (ω) × H 1 (ω) × L 2 (
ω) be such that (ϑ(y) + η i (y)a i (y)) • q ≥ 0 for a.a. y ∈ ω 0 . Let 0 < h < d be given, and let 0 < < h 2 /2. In correspondence of one such number h, define the number = (h) > 0 by

0 < ≤ h 2   inf h>0 ρ∈{1,2} min y∈ω1 (δ ρh (-ϑ • q)) 2 max{ a i • q C 0 (ω) ; 1 ≤ i ≤ 3}   . (5)
Assume that in correspondence of such a number there exists an element ξ = (ξ i ) ∈ C 1 (ω) such that (ϑ(y) + ξ i (y)a i (y)) • q ≥ 0 for all y ∈ ω 0 and such that |ξ i (y) -η i (y)| ≤ 3 for a.a. y ∈ ω 1 and all 1 ≤ i ≤ 3. (6)

Then, the vector field η defined in a way such that η ,j a j := η i a i + ϕ 1 δ ρh (ξ i a i ),

is such that η = (η ,i ) ∈ H 1 (ω) × H 1 (ω) × L 2 (ω)
and (ϑ(y) + η ,i (y)a i (y)) • q ≥ 0 for a.a. y ∈ ω 1 .

Proof. In the same spirit as Lemma 2, we have that for a.a. y ∈ ω 1 η ,j (y)a j (y)

• q ≥ η i (y)a i (y) • q - ϕ 1 (y) h 2 (ϑ(y + he ρ ) + ϑ(y -he ρ )) • q - 2 ϕ 1 (y) h 2 ξ i (y)a i (y) • q ≥ η i (y)a i (y) • q - ϕ 1 (y) h 2 (ϑ(y + he ρ ) + ϑ(y -he ρ )) • q - 2 ϕ 1 (y) h 2 max{ a j • q C 0 (ω) ; 1 ≤ j ≤ 3} 3 i=1 |ξ i (y) -η i (y)| - 2 ϕ 1 (y) h 2 η i (y)a i (y) • q = 1 - 2 ϕ 1 (y) h 2 η i (y)a i (y) • q - ϕ 1 (y) h 2 (ϑ(y + he ρ ) + ϑ(y -he ρ )) • q - 2 ϕ 1 (y) h 2 max{ a j • q C 0 (ω) ; 1 ≤ j ≤ 3} 3 i=1 |ξ i (y) -η i (y)| ≥ -1 - 2 ϕ 1 (y) h 2 ϑ(y) • q - ϕ 1 (y) h 2 (ϑ(y + he ρ ) + ϑ(y -he ρ )) • q - 2 ϕ 1 (y) h 2 max{ a j • q C 0 (ω) ; 1 ≤ j ≤ 3} 3 i=1 |ξ i (y) -η i (y)| ≥ -ϑ(y) • q + ϕ 1 (y)δ ρh (-ϑ(y) • q) - 2 ϕ 1 (y) h 2 max{ a j • q C 0 (ω) ; 1 ≤ j ≤ 3} ≥ -ϑ(y) • q + ϕ 1 (y)δ ρh (-ϑ(y) • q) -ϕ 1 (y)   inf h>0 ρ∈{1,2} min y∈ω1 (δ ρh (-ϑ(y) • q))   ≥ -ϑ(y) • q,
where the second last inequality derives from the definition of in ( 5) and the assumed concavity of ϑ • q over the set ω 0 .

We observe that the number defined in (5) depends on h, and is such that = O(h 2 ) -big "Oh" in Landau notation.

For treating the case where the concavity assumption fails, we need the following auxiliary result.

Lemma 3. Let the function ϑ = (ϑ 1 , ϑ 2 , ϑ 3 ) ∈ C 1 (ω; E 3 ) be such that ϑ • q > 0 in ω.
Then, for every y 0 = (y 0,1 , y 0,2 ) ∈ ω, there exists a neighbourhood U of y 0 and numbers B ∈ R, B 0 > 0, r > 0 such that the function

(-ϑ • q + B)g is convex in U, where g(y) := 1 - 1 2 2 ρ=1
exp(ry ρ -ry 0,ρ ) for all y = (y 1 , y 2 ) ∈ ω.

Moreover, it results g(y) ≥ B 0 , for all y ∈ U.

Proof. Fix y 0 ∈ ω. Owing to the fact that ϑ ∈ C 1 (ω; E 3 ) and ϑ • q > 0 in ω, we can find numbers B ∈ R and T > 0, and a neighbourhood U 0 of y 0 such that -ϑ(y) • q + B ≤ -T < 0 for all y ∈ U 0 .

For all y = (y 1 , y 2 ) ∈ ω, define the function:

Π(y) := 2 ρ=1
exp(ry ρ -ry 0,ρ ).

We recall that a function is convex if and only if it is convex along any lines that intersect the function domain (cf., e.g., page 67 of [START_REF] Boyd | Convex optimization[END_REF]). In other words, checking the convexity of the function (-ϑ • q + B)g in U amounts to checking that, for all v = (v 1 , v 2 ) ∈ R 2 and all y = (y 1 , y 2 ) ∈ U, the function

H(t) := (-ϑ(y + tv) • q + B)g(y + tv), t ∈ R, (y + tv) = (y 1 + tv 1 , y 2 + tv 2 ) ∈ U, is convex. Let us fix an arbitrary point y ∈ U, a vector v = (v 1 , v 2 ) ∈ R 2
and a scalar t with the aforementioned properties. By direct computation, we have that

dH dt (t) = [-(∇ϑ 1 (y + tv)) • vq 1 -(∇ϑ 2 (y + tv)) • vq 2 -(∇ϑ 3 (y + tv)) • vq 3 ] g(y + tv) + [-ϑ(y + tv) • q + B] d dt (g(y + tv)) ,
and that

d 2 H dt 2 (t) = -v T (∇ 2 ϑ 1 (y + tv) + ∇ 2 ϑ 2 (y + tv) + ∇ 2 ϑ 3 (y + tv))v g(y + tv) + 2 [-∇ϑ 1 (y + tv) • vq 1 -∇ϑ 2 (y + tv) • vq 2 -∇ϑ 3 (y + tv) • vq 3 ] d dt (g(y + tv)) + [-ϑ(y + tv) • q + B] d 2 dt 2 (g(y + tv))
. By the definition of g, we have that dg dt (y + tv) = ∇g(y + tv)

• v = - 1 2 dΠ dt (y + tv) = - r 2 (v 1 + v 2 )Π(y + tv),
and that

d 2 g dt 2 (y + tv) = v T ∇ 2 g(y + tv)v = - r 2 2 (v 1 + v 2 ) 2 Π(y + tv).
Thanks to the global uniform boundedness of the first and second derivative of ϑ • q in ω, the properties of the numbers B and T , and the positiveness of the function Π we derive that there exists a positive number M such that:

1 Π d 2 H dt 2 (t) ≥ -M 1 Π - 1 2 -rM + T 2 (v 1 + v 2 ) 2 r 2 .
We observe that, for r sufficiently large,

p 1 (r) := T 2 (v 1 + v 2 ) 2 r 2 -M r -M > 0.
Let us choose a neighbourhood U of y 0 such that U ⊂⊂ U 0 and

|y ρ -y 0,ρ | ≤ 1 2r ln 3 2
for all y ∈ U and all ρ ∈ {1, 2}.

On the one hand, it is immediate to see that, by the monotonicity of the exponential and ( 7

On the other hand, again by [START_REF] Ciarlet | Mathematical Elasticity[END_REF], we have that

1 Π(y) - 1 2 = 2 ρ=1 exp(-r(y ρ -y 0,ρ )) - 1 2 ≤ 1 for all y ∈ U, (9)
and, by virtue of (8), we have that 1 Π(y)

- 1 2 ≥ 2 3 - 1 2 > 0 for all y ∈ U. (10)
In conclusion, putting ( 9) and ( 10) together, we have that

1 Π(y) - 1 2 ≤ 1 for all y ∈ U. (11)
An application of (11) immediately gives

1 Π d 2 H dt 2 (t) ≥ -M 1 Π - 1 2 -M r + T 2 (v 1 + v 2 ) 2 r 2 ≥ p 1 (r) ≥ 0.
By virtue of (8), we also obtain g(y) ≥ 1 -Π/2 ≥ 1/4, for all y ∈ U. Letting B 0 := 1/4 completes the proof.

Let us observe that the assumption ϑ • q > 0 in ω in Lemma 3 is related to the result stated in Lemma 1.

In the next preparatory lemma, we prove some convergence properties of the operator D ρh for functions. As a consequence of Theorem 6.6-4 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF] we have that

C 1 (ω) • H 1 (ω) = H 1 (ω). Lemma 4. Let {v k } k≥1 be a sequence in C 1 (ω) that converges to a certain element v ∈ H 1 (ω) with respect to the norm • H 1 (ω)
. Then, we have that for all 0 < h < d and all ρ ∈ {1, 2},

D ρh v ∈ H 1 (ω 0 ) and D ρh v k → D ρh v in H 1 (ω 0 ) as k → ∞. (12)
Proof. Fix 0 < h < d and fix ρ ∈ {1, 2}. Clearly, by definition of D ρh we have that

D ρh v ∈ L 2 (ω 0 ). Let us show that D ρh is continuous. A direct computation gives ω0 |D ρh v k (y) -D ρh v(y)| 2 dy = ω0 v k (y + he ρ ) -v k (y) h - v(y + he ρ ) -v(y) h 2 dy ≤ 4 h 2 v k -v 2 L 2 (ω) → 0 as k → ∞.
We now show that D ρh v ∈ H 1 (ω 0 ). For each β ∈ {1, 2} and all ψ ∈ D(ω 0 ), we have that

ω0 D ρh v(y)∂ β ψ(y) dy = 1 h ω0 v(y + he ρ )∂ β ψ(y) dy - 1 h ω0 v(y)∂ β ψ(y) dy = 1 h ω0+heρ v(z)∂ β ψ(z -he ρ ) dz - 1 h ω0 v(y)∂ β ψ(y) dy = - ω0 v(y)∂ β (D -ρh ψ)(y) dy = ω0 ∂ β v(y)D -ρh ψ(y) dy = - ω0 (D ρh (∂ β v))(y)ψ(y) dy,
where the second equality takes into account the change of variables z = y + he ρ , the third equality holds for 0 < h < d, being ψ ∈ D(ω 0 ), the fourth equality is owing to the fact that v ∈ H 1 (ω 0 ) and, finally, the fifth equality is owing to the properties of D ρh . Since, by definition of D ρh , we have that D ρh (∂ β v) ∈ L 2 (ω 0 ), we then obtain that

D ρh v ∈ H 1 (ω 0 ) with ∂ β (D ρh v) = D ρh (∂ β v) ∈ L 2 (ω 0 ).
Finally, we observe that the assumed convergence

v k -v H 1 (ω) → 0 as k → ∞ gives: ω0 |∂ β (D ρh v k ) -∂ β (D ρh v)| 2 dy = ω0 |D ρh (∂ β (v k -v))| 2 dy = ω0 ∂ β v k (y + he ρ ) -∂ β v k (y) h - ∂ β v(y + he ρ ) -∂ β v(y) h 2 dy ≤ 4 h 2 ∂ β v k -∂ β v 2 L 2 (ω) → 0 as k → ∞, thus showing that ∂ β (D ρh v k ) → ∂ β (D ρh v) in L 2 (ω 0 ) as k → ∞.
In conclusion, we have obtained that for a given 0 < h < d and a given ρ ∈ {1, 2},

D ρh v k -D ρh v H 1 (ω0) ≤ 2 h v k -v H 1 (ω) → 0 as k → ∞,
and this completes the proof.

As a remark, we observe that the convergence in ( 12) is to be read as follows: Let h > 0 be given. For each ˜ > 0 (in general, independent of h) there exists an integer K = K(h, ρ, ˜ ) ≥ 1 such that for all k ≥ K it results

D ρh v k -D ρh v H 1 (ω0) < ˜ .
We are now ready to prove the main result of this paper. Theorem 6. Assume that the "density property" stated in Theorem 4, namely,

ŨM (ω) is dense in U M (ω) with respect to the norm of • H 1 (ω)×H 1 (ω)×L 2 (ω)
holds. Assume also that the vector f ε = (f i,ε ) ∈ L 2 (Ω ε ), defining the applied body force density, is such that f i,ε ∈ H 1 (Ω ε ).

Then, the solution

ζ ε of Problem P ε M (ω) is of class V M (ω) ∩ H 2 loc (ω) × H 2 loc (ω) × H 1 loc (ω). Proof.
For sake of clarity we break the proof into two steps, numbered (i) and (ii). Fix ϕ ∈ D(U) such that supp ϕ ⊂⊂ U 1 ⊂⊂ U, for some neighbourhood U 1 of y 0 , such that 0 ≤ ϕ ≤ 1, and such that ϕ ≡ 1 in a compact strict subset of its support. Following the notation of Corollary 1, here we let ω 0 := U, and ω 1 := U 1 . In what follows, the number 0 < h < d is fixed, where the number d has been defined in [START_REF] Agmon | Estimates near the boundary for solutions of elliptic partial differential equations satisfying general boundary conditions[END_REF].

(i) Construction of a candidate test field. Let 0 < < h 2 /2. Let y 0 ∈ ω and let U be the neighbourhood of y 0 constructed in Lemma 3 in correspondence of the immersion θ which is, by assumption, of class C 3 (ω; E 3 ). Let B ∈ R and the function g be defined as in Lemma 3. The vector field

θ := --θ + Bq |q| 2 g, is of class C 3 (ω; E 3
) and, by Lemma 3, it is such that θ • q is concave in U. Without loss of generality, we can assume that U 1 is a domain; otherwise, we take a polygon Q such that supp ϕ ⊂⊂ Q ⊂⊂ U 1 and we rename it U 1 without any loss of rigour. Fix = (h) > 0 in the same spirit as (5) of Corollary 1, i.e.,

0 < ≤ h 2    inf h>0 ρ∈{1,2} min y∈U1 (δ ρh (-θ • q)) 2 max{ a i • q C 0 (ω) ; 1 ≤ i ≤ 3}   , (13)
and observe that, in the notation of Corollary 1, here we have let ϑ := θ.

Since, by Lemma 3, we have g ≥ B 0 > 0 in U, and since ζ ε ∈ U M (ω), writing q in terms of the contravariant basis gives

ζ ε i (y)a i (y) + Bq |q| 2 g(y) • q = ζ ε i (y) + Bq • a i (y) |(q • a (y))a (y)| 2 a i (y) • (q • a j (y))a j (y) g(y)
≥ (-θ(y) • q + B)g(y) = -θ(y) • q, for a.a. y ∈ ω. [START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF] Let {ζ ε, } ≥1 ⊂ U M (ω)∩D(ω) be one of the sequences generated by the "density property". As a result of the assumed "density property", there exists a positive integer k 0 = k 0 (ε, , U) for which Besides, in correspondence of any integer k ≥ k 0 , the assumed "density property" and the same computations yielding [START_REF] Ciarlet | An obstacle problem for elliptic membrane shells[END_REF] give

ζ ε -ζ ε,k H 1 (ω)×H 1 (ω)×L 2 (ω) ≤ 3 max{g(y); y ∈ U} , |ζ ε,k i (y) -ζ ε i (y)| ≤ 3
-θ(y) • q ≤ ζ ε,k i (y) + B q • a i (y) |q| 2 a i (y)g(y) • q, for all y ∈ ω. (15)
Again, with respect to the notation of Corollary 1, we hereby let

ξ i := ζ ε,k i + B q • a i |q| 2 g, so that ξ ∈ C 1 (ω), η i := ζ ε i + B q • a i |q| 2 g, so that η ∈ H 1 (ω) × H 1 (ω) × L 2 (ω),
and we observe that a straightforward computation gives

|ξ i (y)-η i (y)| ≤ (max{g(y); y ∈ U})|ζ ε,k i (y)-ζ ε i (y)| ≤ 3
, for all 1 ≤ i ≤ 3 and a.a. y ∈ ω.

(16) Thanks to the concavity of θ•q in U, the definition in [START_REF] Ciarlet | Un problème de confinement pour une coque membranaire linéairement élastique de type elliptique[END_REF], and ( 14)-( 16), letting ϕ 1 = ϕ 2 in the notation of Corollary 1, we obtain (once again up to re-writing q in terms of the contravariant basis) by Corollary 1,

-θ(y) • q ≤ ζ ε i (y)a i (y) + Bq |q| 2 g(y) • q + ϕ 2 (y)δ ρh ζ ε,k i (y)a i (y) + Bq |q| 2 g(y) • q, (17)
for a.a. y ∈ U 1 . Dividing (17) by g and then subtracting B from each member of [START_REF] Ciarlet | An existence and uniqueness theorem for the two-dimensional linear membrane shell equations[END_REF] gives

-θ(y) • q ≤ ζ ε i (y)a i (y) • q + (g(y)) -1 ϕ 2 (y)δ ρh ζ ε,k i (y)a i (y) + Bq |q| 2 g(y) • q for a.a. y ∈ U 1 . (18)
We then define the vector field ζε,k in a way such that ζε,k ,j (y)a j (y) = ζ ε i (y)a i (y)

+ (g(y)) -1 ϕ 2 (y)δ ρh ζ ε,k i (y)a i (y) + Bq |q| 2 g(y)
for a.a. y ∈ U 1 .

(19) A direct application of ( 18) and [START_REF] Frehse | Zum Differenzierbarkeitsproblem bei Variationsungleichungen höherer Ordnung[END_REF] gives that θ(y) • q + ζε,k ,i (y)a i (y) • q ≥ 0, for a.a. y ∈ U 1 . By virtue of the definition of the covariant and contravariant bases, we have that ζε,k

∈ H 1 (U 1 ) × H 1 (U 1 ) × L 2 (U 1 ). (ii) The solution ζ ε is of class H 2 loc (ω) × H 2 loc (ω) × H 1 loc (ω)
. Since, by definition, we have that ζε (y) = ζ ε (y) for a.a. y ∈ U 1 \ supp ϕ, we are in a position to extend ζε,k by ζ ε outside U 1 ; we denote this extension by ζε,k . Besides, by virtue of (i), such a vector field satisfies the geometrical constraint. As a result, we can infer that ζε,k ∈ U M (ω).

Observe that the assumption that the f 3,ε ∈ H 1 (Ω ε ) implies, by Theorem 4.2-1(b) of [START_REF] Ciarlet | Mathematical Elasticity[END_REF], that the transverse component p 3,ε of the vector p ε = (p j,ε ) appearing in the statement of Problem

P ε M (ω) is of class H 1 (ω). Specialising η = ζε,k in Problem P ε M (ω) gives: -ε U1 a αβστ γ στ (ζ ε )γ αβ g-1 ϕ 2 δ ρh ζ ε,k i a i + Bq |q| 2 g • a j 3 j=1 √ a dy ≤ - U1 p j,ε g-1 ϕ 2 δ ρh ζ ε,k i a i + Bq |q| 2 g • a j √ a dy. (20)
Let us define the translation operator E in the canonical direction e ρ of R 2 and with increment size 0 < h < d for a smooth enough function v : U 1 → R by

E ρh v(y) := v(y + he ρ ), E -ρh v(y) := v(y -he ρ ).
Moreover, the following identities can be easily checked out (cf. [START_REF] Frehse | Zum Differenzierbarkeitsproblem bei Variationsungleichungen höherer Ordnung[END_REF] and [START_REF] Piersanti | On the improved interior regularity of the solution of a fourth order elliptic problem modelling the displacement of a linearly elastic shallow shell subject to an obstacle[END_REF]):

D ρh (vw) = (E ρh w)(D ρh v) + vD ρh w, (21)
D -ρh (vw) = (E -ρh w)(D -ρh v) + vD -ρh w, (22)
δ ρh (vw) = wδ ρh v + (D ρh w)(D ρh v) + (D -ρh w)(D -ρh v) + vδ ρh w. (23)
Given two functions f 1 , f 2 ∈ L 2 (U 1 ), we say that f 1 f 2 if there exists a constant C > 0 independent of h and eventually dependent on

p ε H 1 (ω) , g, its reciprocal g-1 , ζ ε H 1 (ω)×H 1 (ω)×L 2 (ω)
, ϕ and its derivatives such that:

U1 f 1 √ a dy ≤ U1 f 2 √ a dy + C(1 + D ρh (ϕζ ε,k ) H 1 (U1)×H 1 (U1)×L 2 (U1) ). ( 24
)
We first show that there exists a constant C > 0 independent of h but eventually dependent on p ε H 1 (ω) , g, its reciprocal g-1 , ζ ε H 1 (ω)×H 1 (ω)×L 2 (ω) , ϕ and its derivatives such that:

- U1 p j,ε g-1 ϕ 2 δ ρh ζ ε,k i a i + Bq |q| 2 g • a j √ a dy ≤ C(1 + D ρh (ϕζ ε,k ) H 1 (U1)×H 1 (U1)×L 2 (U1)
).

An application of formulas ( 3), ( 21)-( 23), Hölder's inequality, the assumed boundedness of ϕ and its derivatives, the fact that p i,ε ∈ H 1 (ω), the smoothness of the immersion θ ∈ C 3 (ω; E 3 ) (and so of the corresponding fundamental forms, viz. section 2), the integration-by-parts formula for the first order finite difference quotient D -ρh (cf., e.g., formula [START_REF] Ciarlet | Obstacle problems for Koiter's shells[END_REF] in Section 6.3.1 of [START_REF] Evans | Partial Differential Equations[END_REF]) and the assumed "density property" (cf. Theorem 5) give:

- U1 p j,ε g-1 ϕ 2 δ ρh ζ ε,k i a i + Bq |q| 2 g • a j √ a dy = - U1 (p j,ε g-1 ϕa j ) • δ ρh ζ ε,k i a i + Bq |q| 2 gϕ -(D ρh ϕ)D ρh ζ ε,k i a i + Bq |q| 2 g -(D -ρh ϕ)D -ρh ζ ε,k i a i + Bq |q| 2 g -(δ ρh ϕ) ζ ε,k i a i + Bq |q| 2 g √ a dy = - U1 D ρh (p j,ε g-1 ϕa j ) • D ρh ζ ε,k 3 a 3 gϕ dy - U1 (p j,ε g-1 ϕa j ) • δ ρh ζ ε,k α a α gϕ dy - U1 (p j,ε g-1 ϕa j ) • δ ρh gϕ Bq |q| 2 dy + U1 (p j,ε g-1 ϕa j D ρh ϕ) • D ρh ζ ε,k i a i + Bq |q| 2 g √ a dy + U1 (p j,ε g-1 ϕa j D -ρh ϕ) • D -ρh ζ ε,k i a i + Bq |q| 2 g √ a dy + U1 (p j,ε ϕa j δ ρh ϕ) ζ ε,k i a i + Bq |q| 2 √ a dy ≤ C p j,ε g-1 ϕa j H 1 (ω) D ρh (ζ ε,k 3 a 3 gϕ) L 2 (U1) + p j,ε g-1 ϕa j L 2 (ω) D ρh (ζ ε,k α a α gϕ) H 1 (U1) + 2 p j,ε g-1 ϕa j L 2 (ω) D ρh ζ ε,k i a i + Bq |q| 2 g L 2 (U1) + p j,ε g-1 ϕa j L 2 (ω) ζ ε,k i a i + Bq |q| 2 g L 2 (U1) ≤ C p j,ε g-1 ϕa j H 1 (ω) ga 3 C 1 (U1;E 3 ) D ρh (ϕζ ε,k 3 ) L 2 (U1) + ϕζ ε,k 3 L 2 (U1) + p j,ε g-1 ϕa j H 1 (ω) ga α C 1 (U1;E 3 ) D ρh (ϕζ ε,k α ) H 1 (U1) + ϕζ ε,k α H 1 (U1) + p j,ε g-1 ϕa j H 1 (ω) a i C 1 (ω;E 3 ) D ρh ζ ε,k i L 2 (U1) + ζ ε,k i L 2 (U1) + (meas ω) g Bq |q| 2 C 1 (U1) ≤ C D ρh (ϕζ ε,k ) H 1 (U1)×H 1 (U1)×L 2 (U1) + ζ ε,k H 1 (U1)×H 1 (U1)×L 2 (U1) + 1 + D ρh ζ ε,k L 2 (U1) ≤ C D ρh (ϕζ ε,k ) H 1 (U1)×H 1 (U1)×L 2 (U1) + ζ ε,k H 1 (U1)×H 1 (U1)×L 2 (U1) + 1 + (E ρh ϕ)(D ρh ζ ε,k ) L 2 (U1) = C D ρh (ϕζ ε,k ) H 1 (U1)×H 1 (U1)×L 2 (U1) + ζ ε,k H 1 (U1)×H 1 (U1)×L 2 (U1) + 1 + D ρh (ϕζ ε,k ) -ζ ε,k D ρh ϕ L 2 (U1) ≤ C 1 + D ρh (ϕζ ε,k ) H 1 (U1)×H 1 (U1)×L 2 (U1) ,
for some constant C > 0 independent of h but eventually dependent on p ε H 1 (ω) , g, its reciprocal g-1 , ζ ε H 1 (ω)×H 1 (ω)×L 2 (ω) , ϕ and its derivatives as in [START_REF] Piersanti | On the improved interior regularity of the solution of a fourth order elliptic problem modelling the displacement of a linearly elastic shallow shell subject to an obstacle[END_REF]. In the previous calculations, we used Hölder's inequality to derive the first estimate, and we used [START_REF] Léger | Mathematical justification of the obstacle problem in the case of a shallow shell[END_REF] to derive the second estimate.

By Lemma 3 we know that the derivatives of g and g-1 are uniformly bounded in U. As a result, we have that (20) and the computations above imply:

-a αβστ γ στ (ζ ε )γ αβ g-1 ϕ 2 δ ρh ζ ε,k i a i + Bq |q| 2 g • a j 3 j=1 0. ( 25
)
Let us now use ( 21)-( 23) for computing

δ ρh (ζ ε,k ϕ). Define A := -ϕ(D ρh ϕ)(D ρh ζ ε,k )- ϕ(D -ρh ϕ)(D -ρh ζ ε,k
) and observe that [START_REF] Mezabia | Modelling of frictionless Signorini problem for a linear elastic membrane shell, Applicable Analysis[END_REF] gives:

ϕδ ρh (ζ ε,k ϕ) + A = ϕ[ϕδ ρh ζ ε,k + ζ ε,k δ ρh ϕ]. (26)
The next step consists in showing that:

-a αβστ γ στ (ϕζ ε )γ αβ (δ ρh (ϕζ ε,k )) -a αβστ γ στ (ζ ε )γ αβ (ϕδ ρh (ϕζ ε,k )). (27)
Recalling the definition of the change of metric tensor components γ αβ (cf. section 2) and recalling that θ ∈ C 3 (ω; E 3 ), we have that the integral U1 -a αβστ γ στ (ϕζ ε )γ αβ (δ ρh (ϕζ ε,k )) √ a dy can be estimated by estimating the following main nine addenda of its. In the evaluation of the following nine terms, the indices are assumed to be fixed, i.e., the summation rule with respect to repeated indices is not enforced in ( 28)-(36) below. First, thanks to an application of Green's formula (cf., e.g., Theorem 6.6-7 of [START_REF] Ciarlet | Linear and Nonlinear Functional Analysis with Applications[END_REF]), we estimate:

U1 -a αβστ ∂ σ (ϕζ ε τ )∂ β (δ ρh (ϕζ ε,k α )) √ a dy = U1 -a αβστ [(∂ σ ϕ)ζ ε τ + ϕ∂ σ ζ ε τ ]∂ β (δ ρh (ϕζ ε,k α )) √ a dy = U1 ∂ β (a αβστ (∂ σ ϕ)ζ ε τ )δ ρh (ϕζ ε,k α ) √ a dy + U1 -a αβστ ∂ σ ζ ε τ [ϕ∂ β (δ ρh (ϕζ ε,k α ))] √ a dy ≤ C ζ ε τ H 1 (U1)×H 1 (U1)×L 2 (U1) D ρh (ϕζ ε,k α ) H 1 (U1) + U1 -a αβστ ∂ σ ζ ε τ ∂ β (ϕδ ρh (ϕζ ε,k α )) √ a dy + U1 a αβστ (∂ σ ζ ε τ )(∂ β ϕ)δ ρh (ϕζ ε,k α ) √ a dy ≤ U1 -a αβστ (∂ σ ζ ε τ )∂ β (ϕδ ρh (ϕζ ε,k α )) √ a dy + C D ρh (ϕζ ε,k α ) H 1 (U1) . (28)
Second, we estimate:

U1 -a αβστ (-Γ ς στ ϕζ ε ς )∂ β (δ ρh (ϕζ ε,k α )) √ a dy = - U1 ∂ β (a αβστ Γ ς αβ ϕζ ε ς )δ ρh (ϕζ ε,k α ) √ a dy ≤ C ζ ε ς H 1 (U1) D ρh (ϕζ ε,k α ) H 1 (U1) , (29)
where the equality holds as a consequence of Green's formula. Third, we estimate:

U1 -a αβστ (-b αβ ϕζ ε 3 )∂ β (δ ρh (ϕζ ε,k α )) √ a dy = U1 a αβστ (b αβ ζ ε 3 )∂ β (ϕδ ρh (ϕζ ε,k α )) √ a dy - U1 a αβστ b αβ ζ ε 3 (∂ β ϕ)δ ρh (ϕζ ε,k α ) √ a dy ≤ U1 a αβστ (b αβ ζ ε 3 )∂ β (ϕδ ρh (ϕζ ε,k α )) √ a dy + C D ρh (ϕζ ε,k α ) H 1 (U1) . (30)
Fourth, we estimate:

U1 -a αβστ ∂ σ (ϕζ ε τ )Γ υ αβ δ ρh (ϕζ ε,k υ ) √ a dy = U1 -a αβστ (∂ σ ϕ)ζ ε τ Γ υ αβ δ ρh (ϕζ ε,k υ ) √ a dy + U1 -a αβστ ϕ(∂ σ ζ ε τ )Γ υ αβ δ ρh (ϕζ ε,k υ ) √ a dy ≤ C D ρh (ϕζ ε,k α ) H 1 (U1) + U1 -a αβστ (∂ σ ζ ε τ )Γ υ αβ [ϕδ ρh (ϕζ ε,k υ )] √ a dy. (31)
Fifth, we straightforwardly observe that:

U1 -a αβστ (-Γ ς στ ζ ε ς ϕ)Γ υ αβ δ ρh (ζ ε,k υ ϕ) √ a dy ≤ C(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) ) + U1 -a αβστ (-Γ ς στ ζ ε ς )Γ υ αβ [ϕδ ρh (ζ ε,k υ ϕ)] √ a dy. (32)
Sixth, we straightforwardly observe that:

U1 -a αβστ (b στ ζ ε 3 ϕ)Γ υ αβ δ ρh (ϕζ ε,k υ ) √ a dy ≤ C(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) ) + U1 -a αβστ b στ ζ ε 3 Γ υ αβ [ϕδ ρh (ϕζ ε,k υ )] √ a dy. (33)
Seventh, we straightforwardly observe that:

U1 -a αβστ (-Γ ς στ ϕζ ε ς )b αβ δ ρh (ζ ε,k 3 ϕ) √ a dy ≤ C(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) ) + U1 -a αβστ (-Γ ς στ ζ ε ς )b αβ [ϕδ ρh (ζ ε,k 3 ϕ)]
where in the last equality we used the integration-by-parts formula for finite difference quotients. Ninth, and last, we straightforwardly observe that

In conclusion, combining ( 28)-(36) together gives:

-a αβστ γ στ (ϕζ ε )γ αβ (δ ρh (ϕζ ε,k )) -a αβστ γ στ (ζ ε )γ αβ (ϕδ ρh (ϕζ ε,k )).

(37) Using ( 23), [START_REF] Regazzoni | Active force generation in cardiac muscle cells: mathematical modeling and numerical simulation of the actinmyosin interaction[END_REF] in the left hand side of (25) and applying (37) gives , where the latter term corresponds to the left hand side of [START_REF] Piersanti | Modeling cardiac muscle fibers in ventricular and atrial electrophysiology simulations[END_REF].

To sum up, we have that -a αβστ γ στ (ϕζ ε )γ αβ (δ ρh (ϕζ ε,k )) 0, which means that there exists a constant C > 0 independent of h and eventually dependent on g, g-1 , ζ ε H 1 (ω)×H 1 (ω)×L 2 (ω) , ϕ and its derivatives such that -ε U1 a αβστ γ στ (ϕζ ε )γ αβ (δ ρh (ϕζ ε,k )) √ a dy ≤ C(1+ D ρh (ϕζ ε,k ) H 1 (U1)×H 1 (U1)×L 2 (U1) ).

Writing the second order finite difference quotient δ ρh in the form (3) and recalling the properties of supp ϕ, we are then in a position to apply the integration-byparts formula for the first order finite difference quotient D -ρh (cf., e.g., formula ( 16) in section 6.3.1 of [START_REF] Evans | Partial Differential Equations[END_REF]), thus getting: where the constant C 1 > 0 in second inequality is associated with the continuity of the bilinear form appearing in Problem P ε M (ω), the third inequality holds by the continuity of D ρh shown in Lemma 4 (up to taking a suitably large k in step (i)), and the last inequality holds by the assumed "density property".

Since, = O(h 2 ), we have that the term 2h -1 (C + εC 1 ) introduced in the above set of inequalities is uniformly bounded by a positive constant.

Therefore, up to renaming the bounding constant Ĉ, the following estimate holds: Since θ ∈ C 3 (ω; E 3 ), we obtain that the estimates in Theorem 1 and the uniform positive definiteness of the two-dimensional fourth order elasticity tensor (cf., e.g., In conclusion, we have obtained that:

D ρh (ϕζ ε ) 2 H 1 (U1)×H 1 (U1)×L 2 (U1) ≤ c Ĉ(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1)
). An application of Theorem 3 of Section 5.8.2 of [START_REF] Evans | Partial Differential Equations[END_REF] together with the fact that ϕ ≡ 1 in a compact strict subset of its support shows that ζ ε ∈ H 2 loc (ω) × H 2 loc (ω) × H 1 loc (ω), as it was to be proved. As a final remark, we point out that the higher regularity for the applied body force density f ε has been used in part (ii) of Theorem 6 in order to apply formulas ( 16) and ( 17) in section 6.3.1 of [START_REF] Evans | Partial Differential Equations[END_REF], and estimate the right hand side of (20) in terms of D ρh (ϕζ ε,k

3 ) L 2 (U1) . This, in general, would not possible if one considered a second order finite difference quotient as in [START_REF] Frehse | Zum Differenzierbarkeitsproblem bei Variationsungleichungen höherer Ordnung[END_REF] and [START_REF] Piersanti | On the improved interior regularity of the solution of a fourth order elliptic problem modelling the displacement of a linearly elastic shallow shell subject to an obstacle[END_REF].

  max{g(y); y ∈ U} , for a.a. y ∈ ω and all 1 ≤ i ≤ 3, for all k ≥ k 0 .

D

  ∂ σ (ϕζ ε τ )b αβ δ ρh (ζ ε,k 3 ϕ) √ a dy = U1 -a αβστ (∂ σ ζ ε τ )b αβ [ϕδ ρh (ζ ε,k 3 ϕ)] √ a dy + U1 -a αβστ (∂ σ ϕ)ζ ε τ b αβ δ ρh (ζ ε,k 3 ϕ) √ a dy = U1 -a αβστ (∂ σ ζ ε τ )b αβ [ϕδ ρh (ζ ε,k 3 ϕ)] ρh (-a αβστ (∂ σ ϕ)ζ ε τ b αβ )D ρh (ϕζ ε,k 3 ) (∂ σ ζ ε τ )b αβ [ϕδ ρh (ζ ε,k 3 ϕ)] √ a dy + C(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) ),

U1-

  a αβστ (b στ ζ ε 3 ϕ)b αβ δ ρh (ζ ε,k 3 ϕ) √ a dy = U1 -a αβστ (b στ ζ ε 3 )b αβ [ϕδ ρh (ζ ε,k 3 ϕ)] √ a dy ≤ C(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) ) + U1 -a αβστ (b στ ζ ε 3 )b αβ [ϕδ ρh (ζ ε,k 3 ϕ)]√ a dy.

-

  a αβστ γ στ (ϕζ ε )γ αβ (δ ρh (ϕζ ε,k ))-a αβστ γ στ (ζ ε )γ αβ (ϕδ ρh (ϕζ ε,k )) -a αβστ γ στ (ζ ε )γ αβ (g -1 ϕ 2 [gδ ρh ζ ε,k + (D ρh g)(D ρh ζ ε,k ) + (D -ρh g)(D -ρh ζ ε,k )]) -a αβστ γ στ (ζ ε )γ αβ g-1 ϕ 2 gδ ρh ζ ε,k i a i + Bq |q| 2 • a j + (D ρh g) D ρh ζ ε,k i a i + Bq |q| 2 • a j + (D -ρh g) D -ρh ζ ε,k i a i + Bq |q| 2 • a j + ζ ε,k i a i + Bq |q| 2 • a j δ ρh g -a αβστ γ στ (ζ ε )γ αβ g-1 ϕ 2 δ ρh ζ ε,k i a i + Bq |q| 2 g • a j 3 j=1

εa

  U1 αβστ γ στ (D ρh (ϕζ ε ))γ αβ (D ρh (ϕζ ε,k )) √ a dy ≤ C(1+ D ρh (ϕζ ε,k ) H 1 (U1)×H 1 (U1)×L 2 (U1) ).A simple computation thus givesε U1 a αβστ γ στ (D ρh (ϕζ ε ))γ αβ (D ρh (ϕζ ε )) √ a dy ≤ C(1 + D ρh (ϕζ ε,k ) H 1 (U1)×H 1 (U1)×L 2 (U1) ) + ε U1 a αβστ γ στ (D ρh (ϕζ ε ))γ αβ (D ρh (ϕ(ζ εζ ε,k ))) √ a dy ≤ C(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) ) + C D ρh (ϕ(ζ εζ ε,k )) H 1 (U1)×H 1 (U1)×L 2 (U1) + εC 1 D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) D ρh (ϕ(ζ εζ ε,k )) H 1 (U1)×H 1 (U1)×L 2 (U1) ≤ C(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) ) + 2 h (C + εC 1 ) ζ ε,kζ ε H 1 (ω)×H 1 (ω)×L 2 (ω) ≤ C(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) ) + 2 h (C + εC 1 ) ,

εa

  U1 αβστ γ στ (D ρh (ϕζ ε ))γ αβ (D ρh (ϕζ ε )) √ a dy ≤ Ĉ(1 + D ρh (ϕζ ε ) H 1 (U1)×H 1 (U1)×L 2 (U1) ).

Theorem 3 . 3 . 2 H 1 ((ϕζ ε ) 2 H 1 (

 332121 -2 of[START_REF] Ciarlet | Mathematical Elasticity[END_REF]) give the existence of a constant c = c(ω, θ, µ) > 0 such thatε U1 a αβστ γ στ (D ρh (ϕζ ε ))γ αβ (D ρh (ϕζ ε )) U1) + D ρh (ϕζ ε 3 ) U1)×H 1 (U1)×L 2 (U1) .

E-mail address: ppiersan@iu.edu


Xingyuan Xu

Jiayang Wu

Sai T Chu

Brent E Little

Roberto Morandotti

Arnan Mitchell

David J Moss
email: dmoss@swin.edu.au

-Combs Xingyuan Xu

Microwave and RF Photonic Applications of Integrated Kerr Micro-Combs

Keywords:

come

Photonic integrated circuits that exploit nonlinear optics for all-optical signal processing have been demonstrated, particularly in silicon, including all-optical logic [1], demultiplexing from 160Gb/s [2] to over 1Tb/s [3], to optical performance monitoring using slow light at speeds of up to 640Gb/s [4][5], all-optical regeneration [6], and many others. The efficiency of all-optical devices depends on the waveguide nonlinear parameter, γ =  n2 / c Aeff. Although silicon can achieve extremely high nonlinearities, it suffers from high nonlinear losses due to twophoton absorption (TPA) and the resulting free carriers. Even if the free carriers are eliminated by p-i-n junctions, silicon's poor nonlinear figure of merit (FOM = n2 / (β λ), where β is the TPA) of 0.3 in the telecom band is too low to achieve high performance. While TPA can be turned to advantage for all-optical functions [7][8][9], for the most part silicon's low FOM in the telecom band is a limitation.

In 2008-2010, new CMOS compatible platforms for nonlinear optics were introduced, including Hydex [10 -18] and silicon nitride [19]. These platforms exhibit negligible nonlinear absorption in the telecom band, and have revolutionized micro-resonator optical frequency combs. The first integrated CMOS compatible integrated optical parametric oscillators were reported in 2010 [11,19], showing that Kerr frequency comb sources could be realised in chip form by using ring resonators with relatively modest Q-factors. Following this, a stable modelocked laser with pulse repetition rates at 200GHz was reported [12]. The success of this platform arises from its very low linear loss, its moderate nonlinearity (γ  233W -1 km -1 ) and negligible TPA [18].

Here, we review radio-frequency (RF) and microwave photonic applications of micro-combs. Optical frequency combs, and micro-combs in particular, are ideal sources for microwave and RF transversal filters [20]. Radar mapping, measurement, imaging as well as the realization of advanced modulation formats for digital communications, require the generation, analysis and processing of analogue RF signals where both the amplitude and phase of the signals are important. In-phase and quadrature-phase components of a signal can be obtained via a photonic based Hilbert transform [21]. Figure 1 shows a schematic of our ring resonator, with a Q factor of 1.2 million and FSR of 200GHz [11]. The comb spectrum has a spacing corresponding to the FSR of the microring, had reasonably high coherence, not operating in the chaotic regime [22]. For the microwave demonstrated here this was more than adequate. Figure 1 shows the ring resonator device and the output micro-comb. Figure 2 shows the experimental setup and Figure 3 shows the input and output to a device configured to perform high order RF differentiation. We review a wide range of other novel RF functions based on microcombs [23 -26]. These devices have significant potential for applications in telecommunications and on-chip WDM optical interconnects for computing and many other applications.

Fig. 1 .

 1 Fig. 1. Integrated ring resonator (top left), output comb spectrum (top right).

Fig. 2 .

 2 Fig. 2. Experimental setup. TLS = tunable laser source, BPF=bandpass filter, VNA=vector network analyzer, TCS=thermally controlled stage, MZM=machzenhder modulator.

Fig. 3 .

 3 Fig. 3. Experimental measurements of high order differentiation of a microwave pulse a) input pulse b) 1 st order derivative, c) second order derivative and d) 3 rd order derivative.