Martin Petek
email: martin.petek@polito.it

Javier Rivero
email: javier.rivero@polito.it

Jorge Alberto Tobón Vásquez

Senior Member, IEEE Guido Valerio
email: guido.valerio@sorbonne-universite.fr

Fellow, IEEE Oscar Quevedo-Teruel
email: oscarqt@kth.se

Senior Member, IEEE Francesca Vipiana
email: francesca.vipiana@polito.it

Jorge Alberto

Member, IEEE Tobón Vásquez

Tobon Vasquez

Method of Moments for the Dispersion Modelling of Glide-Symmetric Periodic Structures

Keywords: Glide symmetry, dispersion diagram, computational electromagnetics, Green's function, periodic structures, metasurfaces

come

L'archive ouverte

I. INTRODUCTION

M ODERN telecommunication systems have increasingly higher demands in terms of data rates, due to increase of users and proliferation of video traffic [START_REF]IMT traffic estimates for the years 2020 to 2030[END_REF]. With 5G and the emerging 6G technologies, the trend is to increase bandwidth by moving to higher frequencies [START_REF] Wang | 5G Mobile: Spectrum Broadening to Higher-Frequency Bands to Support High Data Rates[END_REF], [START_REF] Jiang | The Road Towards 6G: A Comprehensive Survey[END_REF]. This results in increased attenuation of electromagnetic waves in materials and in free-space. One way to mitigate increased losses in materials is to utilize fully-metallic structures, thus completely removing wave attenuation in dielectric materials. On the other hand, free-space losses can be mitigated with high gain antennas.

A particular topic of interest, in the recent years, to design highly efficient structures, has been the use of glide-symmetric metasurfaces [START_REF] Quevedo-Teruel | On the Benefits of Glide Symmetries for Microwave Devices[END_REF]. Glide-symmetric structures are periodic and invariant with respect to a translation by half a period and a mirroring. The effects of introducing such higher symmetry on propagation of waves in periodic structures have first been studied in the 60s and 70s [START_REF] Crepeau | Consequences of symmetry in periodic structures[END_REF], [START_REF] Hessel | Propagation in periodically loaded waveguides with higher symmetries[END_REF]. Recently, their properties have been exploited to design highly efficient fully-metallic graded index lens antennas [START_REF] Quevedo-Teruel | Glide-symmetric fully metallic Luneburg lens for 5G communications at Ka-band[END_REF], [START_REF] Fan | Wideband Low-Profile Luneburg Lens Based on a Glide-Symmetric Metasurface[END_REF], leaky wave antennas [START_REF] Zetterstrom | Low-Dispersive Leaky-Wave Antennas for mmWave Pointto-Point High-Throughput Communications[END_REF]- [START_REF] Zeng | Fully Metallic Glide-Symmetric Leaky-Wave Antenna at Ka-Band With Lens-Augmented Scanning[END_REF], and filters [START_REF] Monje-Real | Holey Glide-Symmetric Filters for 5G at Millimeter-Wave Frequencies[END_REF]. Furthermore, glide symmetries have been used for suppression of leakage between two connecting parts of a waveguide [START_REF] Sipus | Glidesymmetric holey structures applied to waveguide technology: design considerations[END_REF] and flanges [START_REF] Ebrahimpouri | Using Glide-Symmetric Holes to Reduce Leakage Between Waveguide Flanges[END_REF].

As a result of increased interest, methods to accurately compute the dispersion diagram (i.e., the propagation constants of Bloch modes) of glide-symmetric structures has also garnered attention from the antenna and microwave communities. A commonly applied technique is mode-matching [START_REF] Valerio | Glide-Symmetric All-Metal Holey Metasurfaces for Low-Dispersive Artificial Materials: Modeling and Properties[END_REF], [START_REF] Alex-Amor | Glide-Symmetric Metallic Structures With Elliptical Holes for Lens Compression[END_REF], which is particularly suited to holey structures. At discontinuities, the coupling of modes in different parts of the structure is done by enforcing boundary conditions. If the structure has known analytical field expansions, mode-matching is a computationally efficient and a straightforward method to implement. Furthermore, at lower frequencies, quasi-static approximation can be used with mode matching to faster obtain the effective refractive index [START_REF] Fischer | Quasi-Static Homogenization of Glide-Symmetric Holey Parallel-Plate Waveguides With Ultra-Wideband Validity[END_REF]. Recently, a hybrid method named multimodal transfer matrix method (MMTMM) has been proposed for analysis of glide-symmetric structures [START_REF] Mesa | Simulation-Assisted Efficient Computation of the Dispersion Diagram of Periodic Structures: A comprehensive overview with applications to filters, leaky-wave antennas and metasurfaces[END_REF]. In this method, waveguide ports are placed at the edge of the unit cell and a full-wave solver is used to compute the scattering parameters. In this step, the definition of multiple modes on each port is necessary for each waveguide in order to obtain correct results. A post-processing procedure is used to obtain the propagation characteristics of Bloch modes in the structure. The method can find both real and complex solutions and has recently been applied to study the stopband attenuation of glide-symmetric unit cells with circular holes [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF].

In this work, we propose a method of moments (MoM) modelling approach tailored to fully-metallic glide-symmetric structures. Method of moments has been previously used for modal analysis of printed structures [START_REF] Bozzi | Efficient calculation of the dispersion diagram of planar electromagnetic band-gap structures by the MoM/BI-RME method[END_REF]- [START_REF] Valerio | Input Impedance of Nonperiodic Sources Exciting 1-D Periodic Shielded Microstrip Structures[END_REF], electromagnetic bandgap structures [START_REF] Barlevy | Characterization of electromagnetic band-gaps composed of multiple periodic tripods with interconnecting vias: concept, analysis, and design[END_REF], and 1-D periodic glide-symmetric parallel-plate waveguides [START_REF] Amari | Accurate analysis of periodic structures with an additional symmetry in the unit cell from classical matrix eigenvalues[END_REF]. In this paper, we present a novel integration kernel for 2-D glide-symmetric structures, which reduces the computational domain to the bottom half of the unit cell only. The new kernel is applied to analyze a rectangular holey structure [START_REF] Valerio | Glide-Symmetric All-Metal Holey Metasurfaces for Low-Dispersive Artificial Materials: Modeling and Properties[END_REF] and a circular holey structure [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF]. The obtained results are compared to a commercial eigenvalue solver and to the MMTMM analysis [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF], finding, in both cases, a very good agreement. Several benefits are offered by the method presented here over other approaches. Most importantly, the MoM is able to compute the attenuation constant in the stopbands and for complex modes, both of which are currently unavailable from commercial eigenvalue solvers. Furthermore, the method provides additional information on the mode parity and can be easily applied to arbitrary geometries, an advantage over mode-matching. As opposed to MMTMM, it does not require a truncation of waveguide modes at the edges of the unit cell [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF]. Hence, less physical insight, necessary to understand which modes to include in the computation, is required to obtain accurate results. Furthermore, since commercial software include waveguideport modes by following a strict order of cutoff frequencies and cannot arbitrarily select the modes to retain, using MMTMM can result in an increased computational effort, since many higher-order modes must be included in the simulation in order to keep a few modes of very high order, as in [START_REF] Castillo-Tapia | Study of Forward and Backward Modes in Double-Sided Dielectric-Filled Corrugated Waveguides[END_REF].

The paper is organized as follows. In Sect. II, the integral equation formulation and the novel Green's function are presented. Then, in Sect. III, the equation discretization, the computation of singular integrals and the zero-search algorithm are described. In Sect. IV, numerical results are compared to a commercial eigenvalue solver and MMTMM analysis. The conclusion and perspectives are in Sect. V. Preliminary results of the proposed modelling of 2-D periodic glide-symmetric structures were presented in [START_REF] Vásquez | Periodic Integral Equation Formulation for the Numerical Analysis of Glide Structures[END_REF].

II. INTEGRAL EQUATION FORMULATION

A. Integral equation and periodic Green's function

Let us consider the task of finding modes in a fullymetallic periodic structure whose surface is described with the geometry S. If losses due to finite conductivity of metal are negligible, the surface can be described as a perfect electric conductor (PEC) in free space without significantly impacting the accuracy of the solution. On S, the tangential component of the electric field is zero:

E tan = -jωA tan -∇ S Φ = 0, (1)
where ω is the angular frequency, ω = 2πf , and the field is decomposed into the vector and the scalar potentials, A and Φ respectively. The vector potential represents the contribution to the E tan due to the electrical surface current density J

A(r) = µ S G(r, r ′ ) • J(r ′ ) dS ′ , (2)
and the scalar potential the contribution of the surface charge density

Φ(r) = - 1 jωε S ∇ ′ • J(r ′ )G(r, r ′ ) dS ′ . (3)
Both integrals are computed over the entire surface S. In ( 2) and (3), ε and µ are the permittivity and the permeability of the surrounding medium, J(r ′ ) is the unknown electrical surface current density at the source point r ′ and G(r, r ′ ) is the scalar Green's function at the observation point r. The dyadic Green's function G is defined as a product of the identity dyad I = x x + ŷ ŷ + ẑ ẑ and the scalar Green's function as

G(r, r ′ ) = I G(r, r ′ ). (4)
Since the structure is periodic and thus spans to infinity in both dimensions, the computation domain can be restricted to the finite size of one single unit cell (the surface S in (2) and (3) becomes the metallic surface of a single cell) if the periodicity is captured in the Green's function. This is achieved by introducing the free space periodic Green's function (FSPGF), that can be written as a double infinite sum of free space Green's functions [START_REF] Valerio | Comparative Analysis of Acceleration Techniques for 2-D and 3-D Green's Functions in Periodic Structures Along One and Two Directions[END_REF]

G(r, r ′ ) = +∞ m=-∞ +∞ n=-∞ G mn , (5)
where

G mn = e -jkt00
•ρmn e -jkRmn 4πR mn .

Here, k t00 is the wave vector controlling the (possibly complex) propagation between different points in the lattice. The distance R mn is given by

R mn = (z -z ′ ) 2 + |ρ -ρ ′ -ρ mn | 2 , (7)
where

ρ mn = m s 1 + n s 2 . (8)
Here, ρ mn describes the translation from the unit cell indexed by m = n = 0 to cell indexed by m and n of the lattice, s 1 and s 2 being the lattice periodicity vectors.

B. Glide-symmetric periodic Green's function

We will now extend the concept of FSPGF to introduce a novel Green's function which exploits the structure's internal symmetry. For mirror-and glide-symmetric structures, we introduce the scalar higher-symmetric periodic Green's function (HSPGF)

G(r, r ′ ) = G B ± G T (9)
where G B represents the bottom array of sources as

G B (r, r ′ ) = +∞ m=-∞ +∞ n=-∞ G mn , (10)
and G T the top array of sources as

G T (r, r ′ ) = e -jkt00•ρg +∞ m=-∞ +∞ n=-∞ G mn,t . (11)
The Green's function is represented in Fig. 1 as a double array of point sources. To complement the scalar HSPGF, its dyadic counterpart is defined as:

G(r, r ′ ) = I G B ± R z G T , (12)
where the mirroring dyad R z is defined as

R z = x x + ŷ ŷ -ẑ ẑ. (13)
In the computation of G T , the additional factor, e -jkt00•ρg , is added due to the translation of the source points of the top part by glide-symmetric structures b 1 = b 2 = 0.5. Only a combination of the values of 0.0 and 0.5 are possible, as in other cases the higher-symmetry is broken and the current of the top part of the structure cannot be expressed with a phase shift and attenuation of the current in the bottom part [START_REF] Hessel | Propagation in periodically loaded waveguides with higher symmetries[END_REF]. The terms G mn,t are defined as

ρ g = b 1 s 1 + b 2 s 2 , (14)
G mn,t = e -jkt00•ρmn e -jkRmn,t 4πR mn,t , (15)
similarly to ( 5), but with the distance R mn changed to

R mn,t = (z + z ′ ) 2 + |ρ -ρ ′ -ρ mn -ρ g | 2 . ( 16
)
The choice of summing or subtracting in [START_REF] Zetterstrom | Low-Dispersive Leaky-Wave Antennas for mmWave Pointto-Point High-Throughput Communications[END_REF] represents two different types of symmetries. For mirror-symmetric structures, these correspond to perfect magnetic conductor (PMC) (summation) and PEC (subtraction) mirrored modes, but for glide-symmetric structures we will refer to them as the plus (+) and minus (-) symmetries.

C. Ewald acceleration

The expressions ( 9) and ( 12) could be used in simulations as they are, but, in practice, their use is prohibited by a very slow convergence [START_REF] Frutos | Application of Kummer's Transformation to the Efficient Computation of the 3-D Green's Function With 1-D Periodicity[END_REF]- [START_REF] Guerin | Combined Method for the Computation of the Doubly Periodic Green's Function[END_REF]. In this work, we extend the Ewald method [START_REF] Ewald | Die berechnung optischer und elektrostatischer gitterpotentiale[END_REF], [START_REF] Oroskar | Efficient computation of the 2D periodic Green's function using the Ewald method[END_REF] to the concept of HSPGF to allow for efficient computation of the Green's function.

First, we consider the evaluation of ( 10) and [START_REF] Zeng | Fully Metallic Glide-Symmetric Leaky-Wave Antenna at Ka-Band With Lens-Augmented Scanning[END_REF] as two separate sums. Then, the individual sums are separately split to spatial and spectral sums [START_REF] Oroskar | Efficient computation of the 2D periodic Green's function using the Ewald method[END_REF], [START_REF] Linton | The Green's function for the two-dimensional Helmholtz equation in periodic domains[END_REF] as

+∞ m,n=-∞ G mn = +∞ m,n=-∞ G E mn + +∞ p,q=-∞ G E pq , (17)
where all sums are double sums with summation indices m, n for the spatial sums and p, q for the spectral sums. When numerically computing the Green's function, the series in [START_REF] Fischer | Quasi-Static Homogenization of Glide-Symmetric Holey Parallel-Plate Waveguides With Ultra-Wideband Validity[END_REF] are truncated. In this work, the truncation limits are determined by the procedure described in [START_REF] Celepcikay | Choosing splitting parameters and summation limits in the numerical evaluation of 1-D and 2-D periodic Green's functions using the Ewald method[END_REF]. For producing the plots in this paper, we have not observed values larger than ±3. The Ewald method introduces the splitting parameter E, which controls the convergence properties. For frequencies where the wavelength is larger than the periodicity of the unit cell, the value of E opt = π/A, which provides fastest convergence, can be used [START_REF] Jordan | An efficient numerical evaluation of the Green's function for the Helmholtz operator on periodic structures[END_REF]. Otherwise, when the wavelength is smaller than the periodicity of the unit cell, the values of the spectral and spatial parts become similar in magnitude but different in sign. Due to subtraction of two large numbers and the finite precision arithmetic, this leads to a loss of significant digits. The issue can be alleviated by changing the value of E at the cost of slower convergence. The procedure from [START_REF] Celepcikay | Choosing splitting parameters and summation limits in the numerical evaluation of 1-D and 2-D periodic Green's functions using the Ewald method[END_REF] is used in this work to choose E which ensures automatically, without the need of setting a threshold value of frequency, a set number of significant digits even at higher frequencies. The values of the summands in [START_REF] Fischer | Quasi-Static Homogenization of Glide-Symmetric Holey Parallel-Plate Waveguides With Ultra-Wideband Validity[END_REF] are:

G E mn = e -jkt00•ρmn 8πR mn e -jkRmn erfc R mn E - jk 2E + e jkRmn erfc R mn E + jk 2E (18)
and

G E mn = e -jktpq•(ρ-ρ ′ ) A e -jkzpq∆ - z erfc jk zpq 2E -∆ - z E + e jkzpq∆ - z erfc jk zpq 2E + ∆ - z E (19)
for [START_REF] Chen | Glide-Symmetric Holey Leaky-Wave Antenna With Low Dispersion for 60 GHz Point-to-Point Communications[END_REF]. For the contribution of the top part of the structure ( 11)

G E mn,t = e -jkt00•ρmn 8πR mn,t e -jkRmn,t erfc R mn,t E - jk 2E + e jkRmn,t erfc R mn,t E + jk 2E (20)
and

G E mn,t = e -jktpq•(ρ-ρ ′ -ρg) A e -jkzpq∆ + z erfc jk zpq 2E -∆ + z E + e jkzpq∆ + z erfc jk zpq 2E + ∆ + z E . ( 21
)
Here,

∆ - z = |z -z ′ |, ∆ + z = |z +z ′ |, erfc is the complementary error function, k zpq = k 2 -k tpq • k tpq , A = ẑ • (s 1 × s 2 )
is the unit cell area and k tpq is the transverse phasing wave vector of the p, q indexed mode:

k tpq = k t00 + 2 p π (s 2 × ẑ) A + 2 q π ( ẑ × s 1 ) A . ( 22
)
The cost of computing the Green's function can be then further reduced by exploiting the properties of the complementary error function

erfc a -jb = erfc a + jb (23)
and

erfc -a + jb = 2 -erfc a + jb , (24)
where a and b are real numbers and the bar represents conjugation. Using ( 23) and ( 24) reduces the evaluations of the complex complementary error function, which is the most computationally expensive task in the evaluation of the Ewald sums.

III. NUMERICAL SOLUTION A. Discretization of the equation

Now we can proceed with numerically solving [START_REF]IMT traffic estimates for the years 2020 to 2030[END_REF]. First, the surface inside of the enclosed structure is discretized into triangles. During meshing, it is necessary to ensure that the edge triangle vertices coincide perfectly after a translation by a period. Then, additional triangles are added, as shown in Fig. 2. These triangles must be exact copies of their periodic counterparts as they ensure continuity of current between two adjacent unit cells. Then, we discretize (1) by expanding the surface current density J(r ′ ) into

J(r ′ ) = N n=0 I n Λ n (r ′ ), ( 25
)
where Λ n is are basis functions and I n their scaling coefficients. Here, Rao-Wilton-Glisson (RWG) basis functions are used, which are defined as linear interpolation functions over triangles pairs with a shared edge [START_REF] Rao | Electromagnetic scattering by surfaces of arbitrary shape[END_REF]. Then, a system of equations is obtained with a testing procedure. In this work, the testing functions used on (2) and ( 3) are also RWG basis functions, performing the so called Galerkin testing. Ultimately, we obtain expressions for

L mn = µ S dS Λ m (r) • S dS ′ G(r, r ′ ) • Λ n (r ′ ) ( 26
)
and

S mn = 1 ε S dS ∇•Λ m (r) S dS ′ G(r, r ′ ) ∇ ′ •Λ n (r ′ ), ( 27
)
which are combined to form the elements of the N × N impedance matrix [Z], where each element is equal to

Z mn = jωL mn + 1 jω S mn . ( 28
)
To evaluate the dispersion diagram of the structure under analysis, we assume no incident field (i.e., no impressed field is necessary for modes to exist in the structure), obtaining the following system of equations

[Z(k t00 , f )][I] = [0] , (29)
where [I] is a N × 1 vector collecting the coefficients I n of ( 25) and [0] is N × 1 vector of zeros. The system in [START_REF] Frutos | Application of Kummer's Transformation to the Efficient Computation of the 3-D Green's Function With 1-D Periodicity[END_REF] has a nontrivial solution when the value of the wave-vector k t00 at a frequency f is such that the matrix [Z] is singular. However, due to finite precision of evaluation of integrals in ( 26) and ( 27), the matrix is never truly singular. When the matrix is evaluated at frequency and wave-vector pairs close to the solutions, its determinant goes to zero. The solutions of the system are therefore the set of values of k t00 at f where the determinant is small and at a local minimum.

B. Evaluation of singular integrals

An important aspect of evaluating the impedance matrix in ( 29) is a careful treatment of integrals ( 26) and [START_REF] Vásquez | Periodic Integral Equation Formulation for the Numerical Analysis of Glide Structures[END_REF]. For basis functions that are overlapping or are adjacent, the integrals become singular and integration by standard Gaussian quadrature rules can not produce an accurate result. To efficiently and accurately compute singular integrals, a singularity cancellation procedure named radial-angular method has been proposed in [START_REF] Khayat | An Improved Transformation and Optimized Sampling Scheme for the Numerical Evaluation of Singular and Near-Singular Potentials[END_REF]. This method was developed to work with free-space Green's function, which only contains one singularity. Thus, in the HSPGF, the application of the method cancels the singularity corresponding to 1/R 00 in ( 6). Since we have a function with multiple singularities, this method must be modified such that the most important singularities are cancelled. Note that the method computes interaction between individual cells (here triangles) and not the entire basis functions. The matrix elements Z mn are thus obtained by summing all four individual interactions between the elements of the two basis functions. In the following, we describe how to apply the radial-angular method in the case of the proposed HSPGF to cancel singularities arising from periodic source replicas and quasi-singularities in ( 6) and (15) (i.e., sources from adjacent unit cells) and the quasi-singularities in the case of a small gap between the surfaces arising from the top sources, contained in [START_REF] Zeng | Fully Metallic Glide-Symmetric Leaky-Wave Antenna at Ka-Band With Lens-Augmented Scanning[END_REF]. Here, we describe below the two possibilities.

First, the periodicity of the function requires a modification which ensures that the singularity closest to the observation element is cancelled. For example, if the source element, depicted in blue in Fig. 3, is on the opposite side of the unit cell than the test element (orange), the periodic image of the source (dashed blue) is the dominant singularity. Thus, the singularity of the periodic image element 1/R 11 is not appropriately cancelled as the radial-angular method cancels only the 1/R 00 singularity. To rectify this problem, these element interactions are treated by first shifting the test elements with substitution of r → r

ρ 11 R 00 R 11 (a)
+ M 0 s 1 + N 0 s 2 . ( 30
)
where the factors of shifting M 0 and N 0 can be obtained through

M 0 = -(r -r ′ ) • (s 2 × ẑ) A + 0.5 (31)
and

N 0 = -(r -r ′ ) • ( ẑ × s 1 ) A + 0.5 . (32)
In Fig. 3, M 0 = N 0 = -1. After the translation, the dominant singularity becomes 1/R 00 and the radial-angular method can be applied to properly cancel the singularity, as seen in Fig. 3(b). Since shifting one of the element causes it to be outside of the original unit cell, the result of the integration is different to a multiplicative constant. Then, the integral is back-propagated using the Floquet theorem:

Z kl mn = Z kl mn e jkt00•(M0s1+N0s2) , (33)
where Z kl mn here represents the contribution of interaction between elements k and l to the total interaction between the two basis functions. The procedure described here could be adapted to shift the source and not the test element to ensure that the 1/R 00 is the dominant singularity. In our case, we choose to shift the test element as it corresponds to the outer loop of the impedance matrix evaluation and is thus easier to implement in our code.

Second, there are cases when a quasi-singular behavior of the top Green's function must be cancelled for accurate integration. This is particularly true for structures with a small gap between the top and bottom plate, where singularities from both contributions to Green's function of top and bottom sources are appreciable, such as the example presented in Fig. 4. Here, the blue and dashed blue triangles are the bottom R 00 R 00,t and top source elements, and the test element (orange) is sufficiently close to both such that 1/R 00 and 1/R 00,t need to be cancelled. In these cases, the radial-angular method is applied twice, once to the contribution of the bottom part of the structure, using only the Green's function in [START_REF] Chen | Glide-Symmetric Holey Leaky-Wave Antenna With Low Dispersion for 60 GHz Point-to-Point Communications[END_REF], and once to the contribution of the top part, with the Green's function in [START_REF] Zeng | Fully Metallic Glide-Symmetric Leaky-Wave Antenna at Ka-Band With Lens-Augmented Scanning[END_REF]. For the bottom part, the radial-angular method is applied directly. For the top part, the test element is first mirrored by mirroring the test element in z, as seen in Fig. 4. The integral is then evaluated with [START_REF] Chen | Glide-Symmetric Holey Leaky-Wave Antenna With Low Dispersion for 60 GHz Point-to-Point Communications[END_REF], which is the Green's function of the bottom part. This way, the roles of R 00 and R 00,t are exchanged and the radial-angular method cancels the singularity.

C. Zero-search algorithm

To find the solutions of (29), we search for all numerical singularities in the lowest eigenvalue, λ min , of the impedance matrix [Z] within a given frequency region. Here, both plus and minus branches of ( 9) and ( 12) must be considered in the search. To accelerate obtaining the results, we use interpolation of the individual matrix elements to obtain a finer sampling of the impedance matrices. We start to fill the dispersion diagram with the case of real modes (i.e., k t00 is real). The first three steps are illustrated in Fig 5 . The procedure is as follows:

Step 1

Step 2

Step 3 Fig. 5. Zero search algorithm for N init = 12 and N interpol = 600. From top to bottom, found singularities are added to the interpolation scheme. In the second step, the higher frequency solution converges, so only the lower frequency solution is added.

1) Compute N init impedance matrices in linearly equispaced frequency points (at cross markers in Fig. 5) in the frequency range of interest. 2) Use spline interpolation to obtain a finer sampling of impedance matrices with N interpol equispaced frequency points. 3) Compute smallest eigenvalue λ min of these matrices.

The magnitude of the eigenvalues is presented in full line and the phase in dashed line in Fig. 5. 4) Find solutions where the phase of λ min crosses zero and the magnitude is a local minimum. 5) Use bisection on the phase of the λ min to refine the solution up to the desired precision.

6) Compute impedance matrices at the found unconverged solutions (plus markers in Fig. 5). 7) Reobtain the solution frequencies with the procedure listed in steps from 2) to 5) and with new impedance matrices added to the set of previous initial points N init . 8) Repeat steps 6) and 7) until the scheme converges. The frequency solution is considered to have converged when:

|f next -f previous | f previous < δ search , ( 34
)
where we select the pairs of previous solution frequency f previous and next solution frequency f next as the closest two frequencies. In this work, the precision value is set to δ search = 10 -3 . The entire procedure is then repeated for the next value of k t00 .

Once real modes are tracked, we can move to the analysis of complex modes: the procedure is analogous, but a two dimensional search is needed. We start by identifying the maximum and minimum frequencies in the computed real part of the dispersion diagram. In many cases, they lie at the edge of the irreducible Brillouin zone (usually Γ, X or M). In these cases, it is possible to do only a one dimensional search as described previously by fixing the real part of components of k t00 according to the considered stopband, and finding frequencies of all singular matrices for the desired level of attenuation α. However, when the beginning of a new complex mode is not at an edge of the Brillouin zone, both real and imaginary parts are unknown. In this case, we fix a frequency value close enough to the last frequency of the mode already tracked. Then, we follow the procedure:

1) Compute N init (initial points) impedance matrices in a grid in the search region in the real and imaginary space of k t (cross markers in Figs. 6 and7). 2) Use natural interpolation [START_REF] Amidror | Scattered data interpolation methods for electronic imaging systems: a survey[END_REF] to obtain finer sampling of impedance matrices in N interpol points.

|real(k t,previous ) -real(k t,current )| |real(k t,previous )| < δ search , (35)
and imaginary parts separately:

|imag(k t,previous ) -imag(k t,current )| |imag(k t,previous )| < δ search . (36)
Thus, both ( 35) and (36) need to be satisfied for achieving convergence.

Step 1

Step 2

Step 3

Fig. 6. Magnitude in first three steps of the complex zero search algorithm with N init = 36 and N interpol = 1600.

Step 1

Step 2

Step 3

IV. NUMERICAL RESULTS

In this section, we present the results of the modelling procedure described in Sects. II and III. Two different geometries, for both mirror-and glide-symmetric unit cell, are considered. Solutions are computed for different stages of mesh refinement to ensure the convergence and stability of the scheme. The obtained results are compared to CST Eigenvalue Solver (CST ES) [START_REF]CST Studio Suite 3D EM Simulation and Analysis Software[END_REF] for propagative modes and to MMTMM [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF] for evanescent modes (including the real part of the complex modes).

Note that, differently from the CST ES and the MMTMM methods, singularities in MoM arise when modes are supported either inside or outside of the structure. Using the methodology presented in this paper, modes corresponding to solutions in the unbounded pin structure in the external region can appear. However, they can be easily identified as having zero fields inside the bounded region and non-zero fields in the external region. As they are not of interest to this work, these solutions are not shown here. Furthermore, the evanescent modes of the interior region not presented in [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF] are also not included for clarity of the plots.

A. Rectangular holey structure

First, we consider a rectangular holey structure, whose meshed bottom part is shown in Fig. 8. This geometry is chosen as it consists of flat surfaces only, hence a direct subdivision of one triangle to four results in a refined mesh with four times more elements. Such division is applied to the original mesh in Fig. 8 can be observed that the accuracy of the solution slightly deteriorates at the highest frequencies. Refining the mesh once (plus markers) and twice (square markers) improves the result. In the mirror-symmetric case, all modes are PECsymmetric, while, in the glide-symmetric case, the first mode is minus-symmetric while the second and third modes are plus-symmetric. Here, the second mode is in fact a Floquet harmonic of the first mode, as reported in [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF] and could also be obtained by extending the k t00 search space of the first mode.

B. Circular holey structure

For the second structure, we aim to verify the applicability of the code to find complex and evanescent modes. For this purpose, the geometry was chosen to compare the results with Fig. 8 from [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF]. Two meshes were generated to have ratio of wavelength to largest triangle edge of approximately 4 and 8 at 100 GHz, corresponding to 255 and 831 basis functions, as shown in Figs. 10(a) and (b). Here, due to the curved surface, the entire structure is remeshed instead of splitting each triangle into four. For the generation of the dispersion diagrams, the frequency range considered is The meshing software used is Gmsh [START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF]. 0.5 GHz to 100 GHz with N init = 20 initial points in the interpolation scheme. For evanescent modes, the frequency ranges are set to coincide with the stopbands at Γ, X and M. The number of initial points N init varies depending on the width of the stopband such that the impedance matrices are computed every 5 GHz and the number of interpolation points is N interpol = 500.

The full dispersion diagram for the mirror-symmetric unit cell can be seen in Fig. 11. The attenuation in the stopband is shown in Figs. 12(a) and (b) for the cases of ΓX and MΓ regions, respectively. Overall, the results agree very well with the reference [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF], with a reduction of accuracy as the frequency increases for the coarsest mesh. Here, the difference is not only due to the larger relative size of the unit cells when compared to the wavelength, but also due to approximation of the curvilinear surface with planar triangular elements. Only PEC-mirrored modes are found in both original (cross markers) and refined (plus markers) meshes, which is expected as the gap between top and bottom plates is small enough to suppress any PMC-mirrored modes. The dispersion diagram for the glide-symmetric unit cell can be seen in Fig. 13. The attenuation in the stopband is presented in Figs. 14(a) and (b) for the cases of ΓX and MΓ, respectively. Four modes are found, two minus-symmetric (yellow and green markers) and two plus-symmetric (orange and purple markers). As before, the original mesh is presented with cross markers and the refined mesh with plus markers. It is worth noting that, due to the extra information of the mode symmetry, utilizing the HSPGF helps with interpretation of the mode evolution. For example, the plus-mirrored evanescent mode in Fig. 14(b) is found at M for the lower frequencies and transitions to a propagative mode at about 30 GHz. Then, with increasing frequency the mode gradually moves from M to Γ. At Γ, it becomes again evanescent and does not transition into a propagative mode in the frequency range considered.

C. Computational time analysis

The computational time to evaluate the impedance matrices is presented in Table I. The structure studied is the rectangular-holey unit cell from Fig. 8 in both mirror and glide-symmetric configurations. Times for both FSPGF and the HSPGF proposed in this paper are presented. All computations were done on a HP ProBook 430 G8 Notebook with Intel i7-1165G7 and 16 GB of RAM. Note that the times presented are obtained with a single-threaded program. As the computational time has slight frequency dependence, the numbers are the average time of 10 evaluations of linearly-spaced frequencies from 0.5 GHz to 60 GHz and rounded to 1 second. For fair comparison, the mesh of the full structure (corresponding to use of FSPGF) is obtained by a mirror and a translation of ( 14) and thus has twice the number of basis functions. It can be observed from Table I that using the HSPGF results in twice faster evaluation of the matrix. While the matrix is reduced to one fourth number of entries, introducing the top part of the structure in HSPGF results in about double the computational effort in evaluating the Green's function. The average times necessary to find all the solutions (frequencies) for a given value of k t00 in Fig. 9 is presented in Table II. The difference in solution times when comparing the mirrorsymmetric and glide-symmetric structures can be attributed to a larger computational load of the zero-search algorithm since the first and the second mode are not separated in the XM region. The simulation time for obtaining the entire dispersion diagram in Figs. 9, 11 and 13 with CST ES is from 7 to 9 minutes for 30 values of k t00 .

V. CONCLUSION AND PERSPECTIVES

In this work, a novel Green's function with a modelling approach to obtain properties of modes of mirror-and glidesymmetric structures has been presented. Using the proposed Green's function with the MoM, we obtained accurate results for both propagative and evanescent modes when compared to results from CST ES and MMTMM [START_REF] Chen | Accurate Characterization and Design Guidelines of Glide-Symmetric Holey EBG[END_REF]. Moreover, the Green's function evaluation accelered via Ewald and the adapted use of the radial-angular method [START_REF] Khayat | An Improved Transformation and Optimized Sampling Scheme for the Numerical Evaluation of Singular and Near-Singular Potentials[END_REF] to correctly cancel singular integrals were presented. The solutions were found via a technique for the search of singularities of a matrix function. Using the proposed Green's function was found to reduce the computational time of one matrix evaluation to one half. Furthermore, the Green's function provides additional information of the mode parity, which was found helpful in understanding and tracking the evolution of modes. The presented method can be easily integrated to an existing MoM code and is able to obtain information not directly available from commercial solvers, such as attenuation in the stopband. In the future work, we aim to extend the code to work with dielectric materials and utilize it to study hexagonal unit cells.

where b 1 Fig. 1 .

 11 Fig. 1. Higher-symmetric periodic Green's function, represented as double array of point sources. Here, b 1 = b 2 = 0.5.

Fig. 2 .

 2 Fig. 2. Meshing process. The inside of the structure is meshed and additional triangles (orange) are added by translating their periodic counterparts (blue) in the original unit cell.

Fig. 3 .

 3 Fig. 3. Translation of the test element (orange) for proper cancellation of periodic singular elements. Before translation (a), the periodic image element indexed 11 (dashed blue) is closest to the test element. After translation (b), the source element indexed 00 (blue) becomes the closest.

Fig. 4 .

 4 Fig. 4. Mirroring of the test element (orange) for proper cancellation of singular elements from G T . The top mesh is faded as it is present through the HSPGF. Before mirroring (a), applying radial-angular method cancels the bottom singularity (full blue). After mirroring (b), the top singularity can be cancelled.

3 )

 3 Compute the minimum eigenvalues of the matrices obtained in previous step. The magnitude and phase are presented with color in Figs. 6 and 7. 4) Find the regions of the solutions by finding the coincidence of 2-D local minimums of magnitude and a change in sign of the phase. 5) Progressively refine the region of the solution around the local minimum until the desired precision in both real and imaginary part of k t is reached. The found solutions are presented as plus markers in Figs. 6 and 7. 6) Evaluate the impedance matrix in the solution point and repeat steps from 1) to 5) until the scheme converges. First three steps can be seen in Figs. 6 and 7. Each time, the finer sampling of matrices is obtained with one more initial point, located at the previous solution. 7) Move to the next frequency. Estimate the next solution with linear extrapolation and center the search region around it. The convergence criterion for termination of the search is applied to both real

Fig. 7 .

 7 Fig. 7. Phase in first three steps of the of the complex zero search algorithm with N init = 36 and N interpol = 1600.

  Fig. Meshes of the rectangular holey structure. From left to right they are (a) original, (b) refined 1 and (c) refined 2 meshes. Here, z min = -1.25 mm and zmax = -0.25 mm. The parameters periodicity are p = 4 mm, depth of hole h = 1.5 mm, the width of hole is w = 3 mm, and the gap is 0.5 mm.

Fig. 9 .

 9 Fig. 9. Dispersion diagrams for the (a) mirror-symmetric and (b) glidesymmetric rectangular holey structures with the geometry simulated and the irreducible Brillouin zone depicted in the insets.

Fig. 10 .

 10 Fig. 10. Rectangular holey structure meshes for the original mesh (a) and the refined mesh (b). Here, zmax = -0.025 mm and z min = -1.025 mm. The dimensions marked are periodicity p = 3.2 mm, hole diameter d = 2.56 mm, hole depth h = 1 mm and gap between the top and bottom plates is 0.05 mm.The meshing software used is Gmsh[START_REF] Geuzaine | Gmsh: A 3-D finite element mesh generator with built-in pre-and post-processing facilities[END_REF].

Fig. 11 .

 11 Fig. 11. Dispersion diagram for the irreducible Brillouin zone and the mirrorsymmetric unit cell, both depicted in the insets.

Fig. 12 .Fig. 13 .

 1213 Fig. 12. (a): attenuation in the ΓX region for the mirror-symmetric unit cell depicted in the inset. The solutions from 30 -60 GHz were found at X, from 76-85 GHz at Γ and the mode from 85 -89 GHz is a complex mode with the real part depicted in Fig. 11. (b): attenuation in the MΓ region for the mirrorsymmetric unit cell depicted in the inset. The solutions from 49 -60 GHz were found at M and from 75 -100 GHz at Γ.

Fig. 14 .

 14 Fig. 14. (a): attenuation in the ΓX region for the glide-symmetric unit cell depicted in the inset. The solutions from 46-59 GHz are located at Γ and from 59-87 GHz are a complex mode with the real part presented in Fig. 13. (b): attenuation in the MΓ region for the mirror-symmetric unit cell depicted in the inset. The solutions from 49-60 GHz were found at M and from 75-100 GHz at Γ.

  COMPUTATIONAL TIME (IN MINUTES) FOR OBTAINING THE SOLUTIONS (FREQUENCIES) IN FIG. 9 FOR FIXED k t00 .

  This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2023.3325207 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

TABLE I AVERAGE

 I COMPUTATIONAL TIME (IN SECONDS) FOR ONE IMPEDANCE MATRIX EVALUATION.

This article has been accepted for publication in IEEE Transactions on Antennas and Propagation. This is the author's version which has not been fully edited and content may change prior to final publication. Citation information: DOI 10.1109/TAP.2023.3325207 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/

ACKNOWLEDGMENT

The authors thank Prof. Donald R. Wilton for fruitful discussions on the topic. We also thank Dr. Qiao Chen for providing the reference data used in Figs. 11-14.

Manuscript received ; revised . This publication is based upon work from COST Action SyMat (CA18223), supported by COST (European Cooperation in Science and Technology). It is supported in part by the Horizon Europe Research and Innovation Program under the GENIUS Project, Marie Skłodowska-Curie Grant under Agreement 101072560, in part by the project PON Research and Innovation "Microwave Imaging and Detection powered by Artificial Intelligence for Medical and Industrial Applications"(DM 1062/21), funded by the Italian Ministry of University and Research (MUR), and in part by "Unite! -University Network for Innovation, Technology and Engineering".

His research interests encompass various areas, including fast integralequation methods and method-of-moment approaches in computational electromagnetics, metamaterials, nanoplasmonics, Green's function regularization, advanced quadrature integration schemes, periodic Green's functions, and high-Q cavities.

Jorge A. Tobon V. (Member, IEEE) received the Degree in electronics engineering from the Universidad de Antioquia, Colombia, in 2010, in a double degree program, the master's degree in electronic engineering and the Ph.D. degree in electronics and telecommunication engineering from the Politecnico di Torino, Italy, in 2010 and 2014, respectively, where he is an Assistant Professor. His main research activities correspond to modeling, designing, and analyzing systems for microwave imaging applications, specifically in the biomedical and food industry fields, and the numerical modeling of complex and nonhomogeneous media propagation. He received the "Premio Latmiral," an award granted by the Italian Society of Electromagnetism in 2018 and the URSI General Assembly Scientific Symposium Young Scientist Award and the Sorrentino URSI Italy National Meeting Young Scientist Paper Award in 2020. He has made scientific contributions to higher symmetries, transformation optics, lens antennas, metasurfaces, and high-impedance surfaces. He is the co-author of more than 130 papers in international journals and 200 at international conferences.

Francesca Vipiana (M'07-SM '13)


Scott Bolesta

Lisa Burry

Marc M Perreault

Céline Gélinas

Kathryn E Smith

Rebekah Eadie

Federico C Carini

Katrianna Saltarelli

Jennifer Mitchell

Jamie Harpel

Ryan Stewart

Richard R Riker

Gilles L Fraser

Brian L Erstad

International Analgesia and Sedation Weaning and Withdrawal Practices in Critically Ill Adults: The Adult Iatrogenic Withdrawal Study in the ICU*

Keywords: iatrogenic, opioid, sedation, sedative, withdrawal

OBJECTIVES:

Iatrogenic withdrawal syndrome (IWS) associated with opioid and sedative use for medical purposes has a reported high prevalence and associated morbidity. This study aimed to determine the prevalence, utilization, and characteristics of opioid and sedative weaning and IWS policies/protocols in the adult ICU population.

DESIGN:

International, multicenter, observational, point prevalence study.

SETTING:

Adult ICUs.

PATIENTS:

All patients aged 18 years and older in the ICU on the date of data collection who received parenteral opioids or sedatives in the previous 24 hours.

INTERVENTIONS:

None.

MEASUREMENTS AND MAIN RESULTS:

ICUs selected 1 day for data collection between June 1 and September 30, 2021. Patient demographic data, opioid and sedative medication use, and weaning and IWS assessment data were collected for the previous 24 hours. The primary outcome was the proportion of patients weaned from opioids and sedatives using an institutional policy/protocol on the data collection day. There were 2,402 patients in 229 ICUs from 11 countries screened for opioid and sedative use; 1,506 (63%) patients received parenteral opioids, and/or sedatives in the previous 24 hours. There were 90 (39%) ICUs with a weaning policy/protocol which was used in 176 (12%) patients, and 23 (10%) ICUs with an IWS policy/protocol which was used in 9 (0.6%) patients. The weaning policy/protocol for 47 (52%) ICUs did not define when to initiate weaning, and the policy/protocol for 24 (27%) ICUs did not specify the degree of weaning. A weaning policy/protocol was used in 34% (176/521) and IWS policy/ protocol in 9% (9/97) of patients admitted to an ICU with such a policy/protocol. Among 485 patients eligible for weaning policy/protocol utilization based on duration of opioid/sedative use initiation criterion within individual ICU policies/protocols 176 (36%) had it used, and among 54 patients on opioids and/or sedatives ≥ 72 hours, 9 (17%) had an IWS policy/protocol used by the data collection day.

CONCLUSIONS:

This international observational study found that a small proportion of ICUs use policies/protocols for opioid and sedative weaning or IWS, and even when these policies/protocols are in place, they are implemented in a small percentage of patients.

Withdrawal from opioids and sedatives administered under medical supervision is known as iatrogenic withdrawal syndrome (IWS) [START_REF] Tobias | Tolerance, withdrawal, and physical dependency after long-term sedation and analgesia of children in the pediatric intensive care unit[END_REF][START_REF] Anand | National Institute of Child Health and Human Development Collaborative Pediatric Critical Care Research Network: Tolerance and withdrawal from prolonged opioid use in critically III children[END_REF]. In the adult critically ill population, which routinely receives opioids and sedatives for supportive care, IWS prevalence has been reported to be between 17% and 100% depending on the study design and population examined [START_REF] Wang | Opioid-associated iatrogenic withdrawal in critically ill adult patients: A multicenter prospective observational study[END_REF][START_REF] Brown | Opioid and benzodiazepine withdrawal syndrome in adult burn patients[END_REF][START_REF] Arroyo-Novoa | Identifying opioid and benzodiazepine withdrawal in trauma intensive care unit (TICU) patients[END_REF][START_REF] Duceppe | Frequency, risk factors and symptomatology of iatrogenic withdrawal from opioids and benzodiazepines in critically ill neonates, children and adults: A systematic review of clinical studies[END_REF][START_REF] Cammarano | Acute withdrawal syndrome related to the administration of analgesic and sedative medications in adult intensive care unit patients[END_REF][START_REF] Hyun | Iatrogenic opioid withdrawal syndrome in critically ill patients: A retrospective cohort study[END_REF]. Risk factors associated with IWS include cumulative dose and duration of administration of 72 hours or more [START_REF] Wang | Opioid-associated iatrogenic withdrawal in critically ill adult patients: A multicenter prospective observational study[END_REF][START_REF] Duceppe | Frequency, risk factors and symptomatology of iatrogenic withdrawal from opioids and benzodiazepines in critically ill neonates, children and adults: A systematic review of clinical studies[END_REF][START_REF] Cammarano | Acute withdrawal syndrome related to the administration of analgesic and sedative medications in adult intensive care unit patients[END_REF]. Negative clinical outcomes have been associated with IWS, including prolonged mechanical ventilation and ICU length of stay [START_REF] Cammarano | Acute withdrawal syndrome related to the administration of analgesic and sedative medications in adult intensive care unit patients[END_REF].

Weaning opioids and sedatives through gradual dose reductions may be a strategy to minimize the risk of IWS in the ICU [START_REF] Sanchez-Pinto | Implementation of a risk-stratified opioid weaning protocol in a pediatric intensive care unit[END_REF][START_REF] Tiacharoen | Protocolized sedative weaning vs usual care in pediatric critically ill patients: A pilot randomized controlled trial[END_REF][START_REF] Vyas | Standardized scoring tool and weaning guideline to reduce opioids in critically ill neonates[END_REF][START_REF] Solodiuk | Effect of a sedation weaning protocol on safety and medication use among hospitalized children post critical illness[END_REF]. Protocolized weaning of opioids and sedatives in pediatric and neonatal ICUs has been shown to decrease overall drug exposure and reduce IWS and ICU length of stay [START_REF] Sanchez-Pinto | Implementation of a risk-stratified opioid weaning protocol in a pediatric intensive care unit[END_REF][START_REF] Tiacharoen | Protocolized sedative weaning vs usual care in pediatric critically ill patients: A pilot randomized controlled trial[END_REF][START_REF] Vyas | Standardized scoring tool and weaning guideline to reduce opioids in critically ill neonates[END_REF].

Despite what is known about the prevalence, risk factors, and effects of IWS in the pediatric and neonatal ICU populations, little has been published regarding opioid and sedative weaning and IWS practices currently employed in adult ICUs. In addition, methods used in the adult ICU population to monitor and assess IWS have not been adequately studied [START_REF] Duceppe | Frequency, risk factors and symptomatology of iatrogenic withdrawal from opioids and benzodiazepines in critically ill neonates, children and adults: A systematic review of clinical studies[END_REF][START_REF] Chiu | Iatrogenic opioid withdrawal in critically ill patients: A review of assessment tools and management[END_REF]. The aim of this study was to determine the prevalence, utilization, and characteristics of opioid and sedative weaning and IWS policies/protocols in the adult ICU population. We hypothesized that opioid and sedative weaning and IWS practices are not standardized and are infrequently employed in adult ICUs.

MATERIALS AND METHODS

Study Design

We conducted an international, multicenter, observational, 1-day point prevalence study. The study is registered with ClinicalTrials.gov (NCT04422808). Sites were recruited through an iterative process via professional networks, Listservs, a study website, and social media. The study was determined to be exempt by the Wilkes University Institutional Review Board (no. 116); each individual site was required to obtain local research ethics approval before study initiation. Study procedures were in accordance with the ethical standards of the responsible committee on human experimentation (institutional or regional) and with the Helsinki Declaration of 1975, as most recently amended. We report study findings in accordance with the STrengthening the Reporting of Observational studies in Epidemiology checklist [START_REF] Von Elm | STROBE Initiative: The Strengthening the Reporting of Observational Studies in Epidemiology (STROBE) statement: Guidelines for reporting observational studies[END_REF].

Study Population

Each participating ICU selected 1 day for data collection between June 1 and September 30, 2021. Criteria for enrollment included all patients aged 18 years or older in the ICU on the data collection date who received any parenteral opioids or sedatives in the previous 24 hours. Parenteral opioid or sedative use was defined as drug administration by either the IV, subcutaneous, intramuscular, or epidural routes. Enrolled patients were assessed for regularly dosed opioids or sedatives, which we defined as either continuous IV infusion, scheduled intermittent injection, or as needed doses with at least half of the possible doses from the previous 24 hours being administered.

Study Data and Clinical Information

All data were collected from the patient medical record and managed using Research Electronic Data Capture (REDCap) electronic data capture tools hosted at Wilkes University [START_REF] Harris | Research electronic data capture (REDCap)-a metadata-driven methodology and workflow process for providing translational research informatics support[END_REF][START_REF] Harris | The REDCap consortium: Building an international community of software platform partners[END_REF]. REDCap is a secure, web-based software platform designed to support data capture for research studies, providing 1) an intuitive interface for validated data capture; 2) audit trails for tracking data manipulation and export procedures; 3) automated export procedures for seamless data downloads to common statistical packages; and 4) procedures for data integration and interoperability with external sources.

Demographic data, opioid and sedative medication, and weaning and IWS assessment data were collected from the medical records of enrolled patients. Order sets, policies, and protocols were considered equivalent regarding their use to guide opioid and sedative utilization, weaning, and assessment and treatment of IWS. Weaning and IWS were not defined since each study site may have used their own definition. Thus, utilization of opioidsedative weaning and IWS policies/protocols were self-reported by study sites. Although the policies/protocols themselves were not evaluated by study investigators, the following data on these policies/protocols were reported by individual sites: provider(s) allowed to perform assessment, timing when weaning policy/protocol initiated, degree of opioid-sedative dose reduction, and IWS assessment tool(s) used. If an ICU reported having a weaning policy/protocol in place, patients admitted to that ICU were considered to be eligible for utilization of the policy/protocol if they met the timing for initiation criteria reported for the policy/protocol. Patients were coded as receiving mechanical ventilation if it was used on the date of data collection. Duration of mechanical ventilation and duration of ICU admission was the time up to the date of data collection.

Outcomes

The primary outcome was the proportion of patients who had parenteral opioids and/or sedatives weaned by an institutional policy/protocol on the data collection day. Secondary outcomes included the proportion of patients who were assessed for IWS with an institutional policy/protocol on the data collection day, proportion of patients in the ICU on the data collection day receiving any parenteral opioids and sedatives, and specific weaning practices.

Statistical Analysis

Data were summarized using frequencies and percentages for categorical variables and means with sd or medians with interquartile ranges (IQRs) for continuous variables depending on data distribution. Characteristics for ICUs and patients were grouped based on utilization of a policy/protocol for weaning or IWS and summarized since this was directly related to the primary outcome. Comparisons were made for ICU and patient demographic data based on the use of either a weaning or IWS policy/protocol. For comparisons between two parametric categorical variables either the chi-square of Fisher exact test was used, and for nonparametric continuous variables the Mann-Whitney U test was applied. Denominators are reported for categorical variable fields with missing data, and the number (n) is noted for continuous variable fields with missing data. Data entry for some fields was limited due to linked responses from other fields.

All descriptive and statistical analyses were performed using IBM SPSS Statistics version 28.0.0 (IBM, Armonk, NY). For comparisons between groups, the a priori alpha level was set at 5%.

RESULTS

A summary of hospital and patient enrollment is provided in Figure 1. Among 104 hospitals, 88 screened 2,437 patients admitted to 236 ICUs for use of any parenteral opioid and/or sedative in the previous 24 hours. Data for one hospital were not included in the final analysis due to a protocol violation where only one or two patients were enrolled from each ICU; thus, it appeared not all patients in this hospital's ICUs on the day of data collection were screened for opioid and sedative use. In addition, 22 patients were excluded due to missing data, and 4 due to age less than 18 years. This left 2,402 patients from 229 ICUs at 87 hospitals that were screened for enrollment, of which 1,506 received any parenteral opioid or sedative in the previous 24 hours and were enrolled.

Study Sites Characteristics

Most hospitals (56.3%) were from the United States and the United Kingdom (Fig. 2), 60.9% were academic medical centers, 21.7% were community not-for-profit hospitals, 4.6% were community for-profit, and 13.8% were government or other nonpublic institutions. The median number of adult ICUs among all hospitals was 3 (IQR, 1-5), and the median number of total adult ICU beds per hospital was 46 (IQR, . The predominant ICU patient population was mixed medico-surgical (40.2%) with interprofessional bedside rounds occurring daily in the majority (53.3%) of ICUs (Supplement Table 1, http://links.lww.com/CCM/ H358). Most participating ICUs (67.7%) employed a closed practice model [START_REF] Valentin | ESICM Working Group on Quality Improvement: Recommendations on basic requirements for intensive care units: Structural and organizational aspects[END_REF][START_REF] Brilli | American College of Critical Care Medicine Task Force on Models of Critical Care Delivery. The American College of Critical Care Medicine Guidelines for the Definition of an Intensivist and the Practice of Critical Care Medicine: Critical care delivery in the intensive care unit: Defining clinical roles and the best practice model[END_REF] and had a mean of 19 (sd ± 9) beds.

There were no significant differences in demographic data between ICUs that reported having an opioid-sedative weaning or IWS policy/protocol and those that did not, except a significantly higher percentage of ICUs that use a policy/protocol also used an order set to guide opioid-sedative therapy and daily interruption (Supplement Table 1, http://links.lww.com/CCM/H358).

Patient Characteristics and Opioid-Sedative Use

The majority of patients (70.6%) screened for enrollment were from North America (Fig. 3). Demographic results for the 1,506 (62.7%) patients enrolled who received any parenteral opioid or sedative in the previous 24 hours are provided in Supplement Table 2 (http://links.lww.com/CCM/H358). The mean age of patients was 56.4 (sd ± 16.0) years, 38.0% were female, 40.8% were admitted primarily for respiratory disease, 12.9% were using opioids before admission, and 29.8% had a positive PCR test for COVID-19. On the day of data collection, the median duration of ICU stay was 6 (IQR; 2-14) days, 30.6% had active acute respiratory distress syndrome, and 74.4% were mechanically ventilated. There was a significant imbalance in various ethnic groups among patients admitted to ICUs that reported having or not having an opioid-sedative weaning or IWS policy/protocol. Also, patients in ICUs that use a policy/protocol were admitted for trauma more frequently and spent more time on mechanical ventilation, while patients in those that do not use a policy/protocol were more often admitted for circulatory system diseases.

Among enrolled patients 1,305 (86.7%) received regularly dosed opioids, 1,151 (76.4%) received regularly dosed sedatives, and 950 (62.1%) received regularly dosed opioids and sedatives (Table 1). Thus, among the 1,506 patients enrolled 201 did not receive regularly dosed opioids, and 355 did not receive regularly dosed sedatives.

Opioid and Sedative Weaning Policies/Protocols

On the day of data collection, among all participating ICUs, a weaning policy/protocol was used in 176 of the 1,506 (11.7%) patients enrolled. Among enrolled patients, 521 (34.6%) were admitted to 90 (39.3%) ICUs with an opioid-sedative weaning policy/protocol. In ICUs with a weaning policy/protocol in place, it was used in 33.8% (176/521) of patients on the day of data collection (Table 2). Although clinical criteria to describe appropriateness for initiation of a weaning policy/protocol were not assessed, 485 (93.1%) patients admitted to ICUs with a weaning policy/protocol were eligible for policy/protocol initiation based on their duration of opioid/sedative use. Therefore, the policy/protocol was being used in 36.3% (176/485) of patients who met the criterion related to duration of opioid/sedative use on the day of data collection.

IWS Policies/Protocols

Among all participating ICUs, a policy/protocol for IWS was being used in 9 of the 1,506 (0.6%) patients enrolled. A total of 97 (6.4%) enrolled patients were admitted to 23 (10.0%) ICUs with an IWS policy/protocol. In ICUs with an IWS policy/protocol in place, it was used in 9.3% (9/97) of patients on the day of data collection. There were 55.7% (54/97) of patients in these ICUs who had been receiving parenteral opioids or sedatives for 72 hours of more by the day of data collection, which has been associated with IWS risk [START_REF] Wang | Opioid-associated iatrogenic withdrawal in critically ill adult patients: A multicenter prospective observational study[END_REF][START_REF] Duceppe | Frequency, risk factors and symptomatology of iatrogenic withdrawal from opioids and benzodiazepines in critically ill neonates, children and adults: A systematic review of clinical studies[END_REF]. Although clinical criteria to describe appropriateness for initiation of an IWS policy/protocol were not assessed, 16.7% (9/54) of patients who received opioids or sedatives for 72 hours or more were utilizing the policy/protocol on the day of data collection.

Weaning and IWS Practices

ICUs with weaning policies/protocols in place generally delegate the responsibility of patient assessment to nurses and physicians, followed by pharmacists and other providers (Table 3). Also, the policies/protocols in the majority of ICUs do not define when to initiate weaning, and most that do define a time indicate less than 72 hours after beginning opioids or sedatives.

The policies/protocols in the majority of ICUs use an assessment scale or tool to determine the degree of weaning, while about one-quarter use a percent or absolute dose reduction, or do not specify.

Among ICUs using an IWS policy/protocol, 13 (56.5%) incorporated an IWS assessment tool: six used the Opioid and Benzodiazepine Withdrawal Score, one the Sophia Observation Symptoms Scale, one the Sophia Benzodiazepine and Opioid Withdrawal Checklist, and five an institutional tool.

Comparison of Opioid and Sedative Utilization Among ICUs

Significantly more patients in ICUs without an opioid-sedative weaning or IWS policy/protocol received regularly dosed sedatives (Table 1). There was no significant difference in use of an opioid-sedative daily interruption order set, opioid-sedative dose reduction in the 24 hours before data collection, percent dose reduction, or initiation of enteral opioids or sedatives between patients in an ICU with an opioid-sedative weaning or IWS policy/ protocol and those in an ICU without a policy/protocol.

DISCUSSION

This international point prevalence study describes current opioid-sedative weaning and IWS management practices in adult ICUs. The results show that opioid-sedative weaning and IWS policies/protocols are not prevalent and are infrequently employed in adult ICUs. This is despite the fact that more than three-quarters of the patients who received parenteral opioids or sedatives had them provided on a regularly dosed basis, with more than two-thirds receiving them as continuous IV infusions. Also, the results demonstrate that in ICUs with a policy/protocol in place about one in three patients who met the ICU's duration of opioid/sedative use criterion for weaning policy/protocol initiation, and one in six who were at risk of IWS due to receiving opioids or sedatives for 72 hours or more [START_REF] Wang | Opioid-associated iatrogenic withdrawal in critically ill adult patients: A multicenter prospective observational study[END_REF][START_REF] Duceppe | Frequency, risk factors and symptomatology of iatrogenic withdrawal from opioids and benzodiazepines in critically ill neonates, children and adults: A systematic review of clinical studies[END_REF], had the policy/protocol used at the time of data collection. Finally, one-quarter of ICUs utilizing a weaning policy/protocol do not specify the degree of opioid-sedative weaning, and the majority of ICUs with a weaning policy/protocol do not specify when to initiate the policy/protocol. This study also found that almost two-thirds of patients admitted to adult ICUs received parenteral opioids or sedatives in the 24 hours before data collection. This provides an estimate of the current prevalence of parenteral opioid and sedative use in critically ill adults. Utilization of parenteral opioids and sedatives in the current study was similar to that previously reported [START_REF] Iapichino | Use of sedative and analgesic drugs in the first week of ICU stay in high-level-ofcare[END_REF][START_REF] Gill | Perceived versus actual sedation practices in adult intensive care unit patients receiving mechanical ventilation[END_REF].

Standardization of analgesia and sedation practices in the adult ICU utilizing protocols, including daily sedation interruption, is recommended by clinical practice guidelines and has demonstrated improved patient outcomes [START_REF] Devlin | Clinical practice guidelines for the prevention and management of pain, agitation/ sedation, delirium, immobility, and sleep disruption in adult patients in the ICU[END_REF]. Current guidelines do not offer insight about best practices for weaning opioids and sedatives or for assessing IWS in adult ICU patients. Current data about weaning protocols and IWS management are largely derived from the pediatric population. Utilization of such protocols in the pediatric population has been shown to decrease exposure to opioids and benzodiazepines, facilitate more calm and awake days, reduce the occurrence of IWS, and shorten hospital stays [START_REF] Curley | Protocolized sedation vs usual care in pediatric patients mechanically ventilated for acute respiratory failure: A randomized clinical trial[END_REF][START_REF] Amirnovin | Implementation of a risk-stratified opioid and benzodiazepine weaning protocol in a pediatric cardiac ICU[END_REF][START_REF] Sneyers | Strategies for the prevention and treatment of iatrogenic withdrawal from opioids and benzodiazepines in critically ill neonates, children and adults: A systematic review of clinical studies[END_REF]. The results of the current study show that ICUs utilizing an opioid-sedative order set, including daily interruption, often incorporate a policy/protocol for opioid-sedative weaning or IWS assessment; however, consistent with other reports, the availability of such policies/protocols does not guarantee their use in individual patients [START_REF] Burry | Analgesic, sedative, antipsychotic, and neuromuscular blocker use in Canadian intensive care units: A prospective, multicentre, observational study[END_REF][START_REF] Borkowska | Nurses' Sedation practices during weaning of adults from mechanical ventilation in an intensive care unit[END_REF]. Policies/protocols for opioid-sedative weaning and IWS should be used along with clinical assessment of individual patients during interprofessional rounds, and not implemented strictly based on generalized predetermined parameters. Additional research in critically ill adults is needed that focuses on enhancing adherence to available policies/protocols and measuring their clinical benefits when incorporated at the bedside.

Interestingly, our results show that more than half of the ICUs that use an IWS policy/protocol incorporate an IWS assessment tool. Although several tools for assessment of IWS have been validated in critically ill children, none have been validated in the adult ICU population. The Withdrawal Assessment Tool-1 was tested in critically ill adults and demonstrated low sensitivity and specificity, leading the authors to conclude that it is not useful to assess IWS in the adult ICU [START_REF] Capilnean | Validation of the withdrawal assessment tool-1 in adult intensive care patients[END_REF]. Until a valid tool is available, it is difficult to provide a rigorous standardized approach to IWS assessment, and this may have contributed to the results reported from the current study. This represents another important gap that warrants further investigation.

Although this study included data from a large number of hospitals in multiple countries, there are several limitations. Data from many countries were provided by a small number of institutions, which limits generalizability of the data for those countries. In addition, data provided by institutions may not necessarily be representative of the practice within a particular country or geographic region. There also may have been bias in study site participation based on the recruitment methods used, with limited participation of sites whose investigators may not utilize Listservs, internet resources, or social media to identify opportunities for participation in multicenter research. Another important study design limitation is that reasons for the abrupt discontinuation of opioids and sedatives, and patient clinical criteria for readiness to wean, were not collected. Lack of data on why opioids and sedatives were abruptly discontinued, and on the clinical readiness of individual patients for weaning, regardless of policy/protocol triggers, could have biased the results when there may have been clinically justifiable reasons for acutely discontinuing these agents and not initiating weaning per policy/protocol triggers. While the terms "weaning" and "IWS" were clearly contextually tied to the use of analgesics and sedatives in the study protocol and case report form, the terms were not specifically defined for site investigators. This limitation may have led to bias in sites indicating if they used a weaning or IWS policy/protocol. An additional limitation is that only 24 hours of opioid and sedative data were collected at one point in time. Although this approach in study design made its completion more feasible, longitudinal data collection would have captured data that better represented the entire ICU stay of patients. Another limitation is that nearly 30% of included patients had a positive PCR test for COVID-19, and published reports suggest that these patients require higher cumulative doses of opioids and sedatives for longer durations, often with the need for more than one agent, compared to non-COVID-19 admissions [START_REF] Tapaskar | Sedation usage in COVID-19 acute respiratory distress syndrome: A multicenter study[END_REF][START_REF] Dreucean | Approach to sedation and analgesia in COVID-19 patients on venovenous extracorporeal membrane oxygenation[END_REF][START_REF] Balakrishna | An examination of sedation requirements and practices for mechanically ventilated critically ill patients with COVID-19[END_REF][START_REF] Wongtangman | SICU Optimal Mobilization Team (SOMT) Group: Association of sedation, coma, and in-hospital mortality in mechanically ventilated patients with coronavirus disease 2019-related acute respiratory distress syndrome: A retrospective cohort study*[END_REF][START_REF] Flinspach | High sedation needs of critically ill COVID-19 ARDS patients-a monocentric observational study[END_REF]. Inclusion of these patients may have biased the results, especially duration of opioid and sedative use, receipt of regularly dosed agents, dose reduction, and initiation of enteral agents. Finally, data were only captured for patients who received parenteral opioids or sedatives. Although we do report on those patients who began enteral opioids and/or sedatives in the 24 hours before data collection, those who were administered only nonparenteral agents, or who received both parenteral and nonparenteral agents, may have been utilizing weaning or IWS policies/protocols, but were not represented in the results of this study.

CONCLUSIONS

Very few ICUs utilize an opioid-sedative weaning or IWS policy/protocol, and for those that do these policies/ protocols are implemented in a small proportion of patients. This is despite that more than three-quarters of adult ICU patients received regularly dosed opioids or sedatives in the 24 hours before data collection. These results demonstrate that there is underutilization of opioid-sedative weaning and IWS assessment practices in a significant proportion of critically ill adults. This is likely related to the absence of standardized methods for opioid-sedative weaning and a validated IWS assessment tool for critically ill adults. Such methods and a tool are required components for any effort trying to establish the clinical value of opioid-sedative weaning and IWS strategies. Current guidelines should advocate for research regarding valid methods and tools to guide opioid and sedative weaning and IWS practices, including the prevalence, assessment, prevention, and treatment of IWS, and clinicians should be aware of the potential risk of IWS with regularly dosed opioids and sedatives.

ACKNOWLEDGMENTS

List of ALERT-ICU collaborators: site principal investigators-

Sara Alcántara, Hospital Universitario Puerta de Hierro Majadahonda, Madrid, Spain. Namareq Aldardeer, King Faisal Specialist Hospital and Research Centre, Riyadh, Saudi Arabia. Abdullah Alhammad, King Saud University Medical City-King Khalid University Hospital, Riyadh, Saudi Arabia. Sulaiman Almohaish, Virginia Commonwealth University, Richmond, VA. Mahmoud Ammar, Yale New Haven Hospital, New Haven, CT. Carmen Mabel Arroyo-Novoa, Puerto Rico Trauma Hosptial, San Juan, Puerto Rico. Don-Kelena Awassi, Hôpital Maisonneuve-Rosemont, Montreal, Canada. Mariana Badano, Hospital Italiano La Plata, La Plata, Argentina. Belinda Badman, Princess Alexandra Hospital, Harlow, United Kingdom. Sanchia Barnes, Manchester University NHS Foundation Trust, Manchester, United Kingdom. Matias Bertozzi, Sanatorio Anchorena San Martin, Buenos Aires, Argentina. Shady Botros, Ninewells Hospital, NHS Tayside, Dundee, Scotland. Gretchen Brophy, Virginia Commonwealth University, Richmond, VA. Caitlin Brown, Mayo Clinic, Rochester, MN. Glen Brown, St. Paul's Hospital, Vancouver, Canada. Mitchell Buckley, Banner University Medical Center Phoenix, Phoenix, AZ. Sarah Burgess, Queen Elizabeth II Health Sciences Centre, Halifax, Canada. Joshua Campbell, Guthrie Robert Packer Hospital, Sayre, PA. Chancey Carothers, Orlando Regional Center, Orlando, FL. Patricio Celario, Sanatorio 9 de Julio, Buenos Aires, Argentina. Lauren Chambers, Vidant Medical Center, Greenville, SC. Gérald Chanques, CHU of Montpellier, Montpellier, France. Sofia Contreras, Vall d'Hebron University Hospital, Barcelona, Spain. Jennifer Cortes, Memorial Hermann-Texas Medical Center, Houston, TX. Gabrielle Cozzi, Beth Israel Deaconess Medical Center, Boston, MA. Michaelia Cucci, Cleveland Clinic Akron General, Akron, OH. Eleonora Cunto, Hospital de Infecciosas F J Muñiz, Buenos Aires, Argentina. Denis Dao, Jewish General Hospital, Montreal, Canada. Aubrey Defayette, Roswell Park Comprehensive Cancer Center, Buffalo, NY. Alexandra Dey, Prince of Wales Hospital, Randwick, New South Wales. Anatilde Diaz, Hospital Misericordia Nuevo Siglo, Buenos Aires, Argentina. Deepali Dixit, Rutgers University Hospital, Newark, NJ. Jean-Noel Dube, Centre hospitalier affilié universitaire régional de Trois-Rivières, Trois-Rivières, Canada. Marc-Alexandre Duceppe, Royal Victoria Hospital-McGill University Health Centre, Montreal, Canada. Amy Dzierba, New York-Presbyterian Hospital, New York, NY. Odran Farrell, University Hospital Plymouth NHSFT, Plymouth, United Kingdom. Pietro Ferranti, Lister Hospital, Stevenage, East and North Herts NHS Trust, Stevenage, United Kingdom. Alexander Flannery, University of Kentucky HealthCare, Lexington, KY. Gabriel Fontaine, Intermountain Medical Center, Murray, KY. Christy Forehand, Augusta University Medical Center, Augusta, GA. Nydia Funes Nelson, Sanatorio del Oeste, Buenos Aires, Argentina. Manon Gillet (CHU Liege, Be), Centre Hospitalier Universitaire de Liege, Liege, Belgium. Gustavo Nicolas Gonzalez, Hospital Italiano De Buenos Aires, Buenos Aires, Argentina. Emma Graham-Clarke, Sandwell and West Birmingham NHS Trust, West Bromwich, United Kingdom. Sarah Green, Royal North Shore Hospital, St. Leonards, New South Wales. Christine Groth, University of Rochester Medical Center, Rochester, NY. Hytham Hamid, AlMoalem Medical City, Jartum, Sudan. Leslie Hamilton, University of TN Medical Center, Knoxville, TN. Joanne Hanley, Antrim Area Hospital, Antrim, Northern Ireland. Jessica

Drs. Bolesta, Burry, Perreault, Gélinas, Riker, Fraser, and Erstad contributed to the study conception and design. Communication and coordination of study sites were performed by Drs. Bolesta, Eadie, Carini, Saltarelli, Mitchell, and Ms. Harpel and Mr. Stewart. Data collection was performed by Drs. Bolesta, Burry, Perreault, Gélinas, Smith, Eadie, Carini, and The remaining authors have disclosed that they do not have any potential conflicts of interest.

For information regarding this article, E-mail: scott.bolesta@ wilkes.edu