Joffret Ml
email: laurent.dacheux@pasteur.fr

Genomic characterization of Sebokele virus 1 (SEBV1) reveals a new candidate species among the genus Parechovirus

Keywords: Sebokele virus 1, picornavirus, parechovirus, Ljungan virus, high throughput sequencing, Illumina sequencing, complete genome, phylogenetic analysis

.

Summary

We determined the genomic features and the taxonomic classification of Sebokele virus 1 (SEBV1), a previously unclassified arbovirus isolated in 1972 from rodents collected in Botambi, Central African Republic. The complete genome sequence was obtained using a deep sequencing approach (Illumina technology) and dedicated bioinformatics workflows for data analysis. Molecular analysis identified SEBV1 as a picornavirus, most closely related to Ljungan viruses of the genus Parechovirus. The genome has a typical Ljungan virus-like organization, including the presence of two unrelated 2A protein motifs. Phylogenetic analysis confirmed that SEBV1 belongs to the parechovirus phylogroup and was most closely related to the Ljungan virus species. However, it appeared clearly distinct from all members of this phylogroup, suggesting that it represents a novel species of the genus Parechovirus.

Main text

The family Picornaviridae encompasses small nonenveloped and positive single-stranded RNA viruses, and includes many human and animal pathogens. The typical genome organization of picornaviruses includes a single open reading frame (ORF) encoding a large protein precursor (polyprotein), preceded by a 5' untranslated region (5'UTR) and followed by a 3'UTR and a polyadenylated tract of variable length [START_REF] Knowles | Picornaviridae. In Virus taxonomy: classification and Nomenclature of Viruses: Ninth Report of the International Committee on Taxonomy of Viruses[END_REF][START_REF] Racaniello | Picornaviridae: the viruses and their replication[END_REF]. The polyprotein contains, in the following order, a non structural leader (L) protein in some cases, then a first domain P1 encoding three structural proteins (VP0, which is autocatalytically cleaved further into VP4 and VP2 in most picornaviruses, VP3 and VP1), and two distinct non structural coding regions (P2 encoding proteins 2A, 2B and 2C ATPase , a multifunctional ATPase, and P3 encoding proteins 3A, 3B VPg , 3C pro , a cysteine protease, and 3D pol , a RNA-dependent RNA polymerase). The family Picornaviridae is currently divided into 12 genera: Aphtovirus, Avihepatovirus, Cardiovirus, Enterovirus, Erbovirus, Hepatovirus, Kobuvirus, Parechovirus, Sapelovirus, Senecavirus, Teschovirus and Tremovirus. However, the diversity of this viral family is far from being fully documented, and numerous novel picornavirus species and genera have recently been described [START_REF] Knowles | Picornaviridae. In Virus taxonomy: classification and Nomenclature of Viruses: Ninth Report of the International Committee on Taxonomy of Viruses[END_REF] (http://www.picornaviridae.com/).

The genus Parechovirus includes two species, human parechovirus (HPeV) and Ljungan virus (LV), each including various types and/or genotypes [START_REF] Johansson | Molecular analysis of three Ljungan virus isolates reveals a new, close-to-root lineage of the Picornaviridae with a cluster of two unrelated 2A proteins[END_REF][START_REF] Joki-Korpela | Parechoviruses, a novel group of human picornaviruses[END_REF][START_REF] Stanway | Parechoviruses[END_REF]. HPeV is a pathogen infecting mainly the gastrointestinal and the respiratory tracts of children [START_REF] Baumgarte | Prevalence, types, and RNA concentrations of human parechoviruses, including a sixth parechovirus type, in stool samples from patients with acute enteritis[END_REF][START_REF] Harvala | Human parechoviruses: biology, epidemiology and clinical significance[END_REF][START_REF] Harvala | Epidemiology and clinical associations of human parechovirus respiratory infections[END_REF][START_REF] Joki-Korpela | Parechoviruses, a novel group of human picornaviruses[END_REF]. Sixteen different types of HPeV and several different genotypes have been identified [START_REF] Knowles | Picornaviridae. In Virus taxonomy: classification and Nomenclature of Viruses: Ninth Report of the International Committee on Taxonomy of Viruses[END_REF] (http://www.picornaviridae.com/). The genome organization of HPeV is typical of the family, except that the VP0 protein appears not to be cleaved and its N-terminus not myristoylated [START_REF] Knowles | Picornaviridae. In Virus taxonomy: classification and Nomenclature of Viruses: Ninth Report of the International Committee on Taxonomy of Viruses[END_REF][START_REF] Stanway | Parechoviruses[END_REF].

Ljungan virus (LV) was first isolated from bank voles (Myodes glareolus) trapped in Sweden [START_REF] Niklasson | Could myocarditis, insulindependent diabetes mellitus, and Guillain-Barre syndrome be caused by one or more infectious agents carried by rodents?[END_REF][START_REF] Niklasson | A new picornavirus isolated from bank voles (Clethrionomys glareolus)[END_REF]. Few isolates of the LV species have been described and only five whole genome sequences are available. They include three Swedish LV: strains 87-012 (the prototype strain) and 174F, both representing genotype 1, and 145SL representing genotype 2 [START_REF] Johansson | Molecular analysis of three Ljungan virus isolates reveals a new, close-to-root lineage of the Picornaviridae with a cluster of two unrelated 2A proteins[END_REF]. Two American LV isolates (M1146 and 64-7855) have also been sequenced [START_REF] Johansson | Molecular characterization of M1146, an American isolate of Ljungan virus (LV) reveals the presence of a new LV genotype[END_REF][START_REF] Tolf | Molecular characterization of a novel Ljungan virus (Parechovirus; Picornaviridae) reveals a fourth genotype and indicates ancestral recombination[END_REF][START_REF] Whitney | Two viruses isolated from rodents (Clethrionomys gapperi and Microtus pennsvlvanicus) trapped in St. Lawrence County, New York[END_REF], and represent distinct genotypes.

Molecular characterization of LV isolates revealed genomic features unusual in picornaviruses: they have only three different structural proteins like all parechoviruses [START_REF] Ekstrom | Replication of Ljungan virus in cell culture: the genomic 5'-end, infectious cDNA clones and host cell response to viral infections[END_REF][START_REF] Johansson | Cell culture propagation and biochemical analysis of the Ljungan virus prototype strain[END_REF][START_REF] Tolf | Characterization of polyclonal antibodies against the capsid proteins of Ljungan virus[END_REF] and two different 2A protein motifs [START_REF] Johansson | Molecular characterization of M1146, an American isolate of Ljungan virus (LV) reveals the presence of a new LV genotype[END_REF][START_REF] Johansson | Molecular analysis of three Ljungan virus isolates reveals a new, close-to-root lineage of the Picornaviridae with a cluster of two unrelated 2A proteins[END_REF]. LV has been suggested to be the aetiological agent of myocarditis and possibly other human diseases [START_REF] Niklasson | Could myocarditis, insulindependent diabetes mellitus, and Guillain-Barre syndrome be caused by one or more infectious agents carried by rodents?[END_REF][START_REF] Niklasson | Association of zoonotic Ljungan virus with intrauterine fetal deaths. Birth defects research Part A[END_REF].

Sebokele virus 1 (SEBV1) was originally isolated in 1972 from the crushed organs of Hylomyscus sp.

(African wood mice) collected in Botambi, Central African Republic [START_REF] Digoutte | Sebokele (An B 1227d), nouvel Arbovirus non groupé[END_REF][START_REF] El Mekki | Characterization of some ungrouped viruses[END_REF][START_REF] Zeller | Electron microscopy and antigenic studies of uncharacterized viruses. I. Evidence suggesting the placement of viruses in families Arenaviridae, Paramyxoviridae, or Poxviridae[END_REF]. The virus was isolated by intracerebral inoculation of suckling mice [START_REF] Digoutte | Sebokele (SEB)[END_REF]. The pathogenic agent appeared to be resistant to chloroform and could be filtered through 220 nm-pore size filters. Attempts to identify the virus using sera against arboviruses and herpesvirus, and by electron microscopy after inoculation of cell cultures, were unsuccessful. The isolate was therefore described to be a non-identified arbovirus. Three similar viral isolates were reported from Muridae in the Central African Republic: two from Hylomyscus sp. and one from Praomys sp. [START_REF] Digoutte | Sebokele (An B 1227d), nouvel Arbovirus non groupé[END_REF].

We characterized the viral genome of SEBV1 by deep and classical Sanger sequencing approaches.

Total RNA (TRI Reagent) was extracted from a vial, dated 1974, of a lyophilized 10% suspension of the original inoculated suckling mouse brain. The cDNA was synthesised and subjected to sequenceindependent amplification as described previously [START_REF] Dacheux | Application of broadspectrum resequencing microarray for genotyping rhabdoviruses[END_REF]. Amplified DNA was then used for Illumina single-read sequencing (69 nt single reads, GAIIx). The Illumina Analysis Pipeline version 1.6 was used for image analysis, base calling, error estimation and demultiplexing.

A total of 8,509,420 reads were obtained, trimmed and mapped to the Mus musculus genome by using CLC Assembly Cell (v. 3.11) with default parameters. Non-mapped reads (147,826 reads) were extracted and assembled into contigs by three different assemblers (CLC novo assembler, VelvetOptimiser and SOAP de novo). The contigs were used independently for BLASTn and BLASTx identity searches against EMBL and UniProt databases, respectively. The blast2taxoclass program (http://mobyle.pasteur.fr/cgi-bin/portal.py#forms::blast2taxoclass) was used for taxonomic classification, and demonstrated that most of the viral contigs were related to the family Picornaviridae, strongly suggesting that SEBV1 belonged to this family. The blast2genoclass program (http://mobyle.pasteur.fr/cgi-bin/portal.py#forms::blast2genoclass) was used for genomic annotation of contigs and high-scoring segment pairs (hsps) matching with picornaviruses were extracted.

Extracted hsps were assembled using Sequencher 5.0 software (Gene Codes Corporation) with default parameters, resulting in the generation of three main contigs covering nearly 99% of the predicted genomic sequences (Fig 1a). Following BLASTx analysis, contigs presented a maximum amino acid (aa) identity of 61 to 65% with parechoviruses.

Attempts to culture the virus in several cell lines were all unsuccessful, so a second vial of SEBV1 was used to finalise the full-length genome sequence determination. Total RNA was extracted, amplified by RT-PCR with primers corresponding to contig sequences, and subjected to Sanger sequencing (Supplementary Table 1). The final full-length sequence served as a reference for mapping reads generated by deep sequencing (CLC Assembly Cell). The average coverage for each nucleotide position was 104 times (lower at the extremities), without any discrepancies being detected (Fig. 1b). SEBV1 genomic sequences were compared to reference nucleotide and protein sequences from the different viral genera downloaded from GenBank. BioEdit software [START_REF] Hall | BioEdit: a user-friendly biological sequence alignment editor and analysis program for Windows 95/98/NT[END_REF] was used to translate DNA and the ClustalX program, version 2.0 (Larkin et al., 2007;[START_REF] Thompson | The CLUSTAL_X windows interface: flexible strategies for multiple sequence alignment aided by quality analysis tools[END_REF] for sequence alignments. Pairwise alignments were performed using CLC Main Workbench 5.7.2 software (CLC bio). Phylogenetic analyses of aa sequences were conducted by the neighbor-joining algorithm using MEGA 4.0.2 [START_REF] Tamura | MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0[END_REF] with a Poisson correction model and 1,000 bootstrap replicates.

Trees were visualised using the FigTree program, version 1.3.1 (available from http://tree.bio.ed.ac.uk/software/figtree). The GenBank accession numbers of the viral sequences used in the phylogenetic analyses are indicated in the trees.

The genome of SEBV1 was 7,537 nucleotides (nt) long (excluding the poly(A) tail), and is organised as follows: a 5'UTR sequence of at least 754 nt, followed by an ORF of 6,702 nt predicted to encode a polyprotein of 2,233 aa, terminated by a 3'UTR of 81 nt and a polyadenylated tract. The exact length of the 5'UTR sequence remains however to be confirmed by functional analysis, as the use of lyophilised stored vials as the starting material did not allow us to assume that the 5' end of the genome sequence was complete. More than 30% of the 5'UTR sequence (nucleotide position 400 to 735) presented nt identities with 5'UTR sequences of parechoviruses, according to BLASTn analysis, with maximum coverage (43%) and nt identity (72%) obtained with LV strain 64-7855. The 3'UTR sequence of SEBV1 seems to be unique, and BLASTn searches found no related sequence. The G+C content of the SEBV1 genome is 45.7%, higher than the G+C contents of LV (41.6-42.5%) and human parechoviruses (39.0-39.5%). BLASTp analysis of the complete deducted polyprotein identified the best matches as LV species (the best being LV strain 87-012), with 62-63% aa identity and 77-78% aa similarity, and human parechoviruses HPeV-1 and -6 exhibited nearly 46% aa identity and 63% aa similarity. The polyprotein sequence of SEBV1 was then aligned with its closest relatives: the Ljungan viruses (strains 87-012, 174F, 145SLG, M1146 and 64-7855) and the human parechoviruses HPeV-1 (strain Harris) and HPeV-6 (strain NII561-2000). With NetPicoRNA prediction [START_REF] Blom | Cleavage site analysis in picornaviral polyproteins: discovering cellular targets by neural networks[END_REF], this allowed the identification of potential cleavage sites and processing of this polyprotein (see Table 1 and Fig. 1c). The organisation of the SEBV1 genome is typical of picornaviruses, with a P1 region encoding the structural proteins and two non structural regions, P2 and P3.

The P1 region is the genomic region of SEBV1 most divergent from the most closely related picornaviruses, with nt identities of 56-57% with LV isolates and 36-37% with the most closely related HPeV (Table 1). This region encodes the classical capsid proteins with, in order, VP0, VP3 and VP1. The P1 region contains the picornavirus capsid protein domain-like (pfam sequence cluster cd00205 -Rhv-like) sequence found in other species of the genera Parechovirus, Avihepatovirus and Aquamavirus, and other unclassified picornaviruses including turdivirus. VP0 of SEBV1 is probably not cleaved into VP4 and VP2. Similarly to the other parechoviruses, the N-terminal end of this capsid protein does not contain the myristylation motif (GxxxS/T) found in most picornaviruses [START_REF] Chow | Myristylation of picornavirus capsid protein VP4 and its structural significance[END_REF]. Like LV and pasivirus 1 (SPaV1), VP0 of SEBV1 has a short N-terminal extremity [START_REF] Sauvage | A member of a new Picornaviridae genus is shed in pig feces[END_REF]. The VP3 of SEBV1 has a long N-terminal extremity enriched with basic residues, like those of parechoviruses [START_REF] Johansson | Molecular analysis of three Ljungan virus isolates reveals a new, close-to-root lineage of the Picornaviridae with a cluster of two unrelated 2A proteins[END_REF]. However, the highly conserved KxKxxRxK motif in the parechovirus P1 region was not found, and is replaced by the motif R 270 FKWTRNN [START_REF] Johansson | Molecular analysis of three Ljungan virus isolates reveals a new, close-to-root lineage of the Picornaviridae with a cluster of two unrelated 2A proteins[END_REF][START_REF] Williams | Evolution and conservation in human parechovirus genomes[END_REF]. VP1 is predicted to be cleaved at the canonical cleavage site DxExNPG 812 P (Fig. 1c) [START_REF] Ryan | Virus-encoded proteinases of the picornavirus supergroup[END_REF]. The parechovirus (PS)ALxAxETG motif was not present in this capsid protein. As in LV, SPaV1 and some HPeVs, SEBV1 lacked the RGD tri-peptide involved in integrin receptor binding, raising questions about the mechanism of cellular entry and the host tropism of this virus [START_REF] Knowles | Picornaviridae. In Virus taxonomy: classification and Nomenclature of Viruses: Ninth Report of the International Committee on Taxonomy of Viruses[END_REF][START_REF] Merilahti | Endocytosis of integrin-binding human picornaviruses[END_REF][START_REF] Sauvage | A member of a new Picornaviridae genus is shed in pig feces[END_REF][START_REF] Williams | Evolution and conservation in human parechovirus genomes[END_REF]. The VP1 protein of SEBV1, like those of LV and SPaV1, contained one N-terminal insertion (8 aa long) and a unique C terminal extension (48 aa long). Similarly to LV, this protein contains the 3C pro cleavage site RRQ 794 GCR predicted by NetPicoRNA at the end of its N-terminus, suggesting the existence of a short 2A-like protein (18 aa) (Fig. 1c).

The non structural P2 region of SEBV1 shows 70-72% aa sequence identity with LV and nearly 48% with HPeV. This region was deduced to encompass the non-structural proteins 2A, 2B and 2C (Table 1 andFig 1c). The 2A protein possesses the highly conserved picornaviral H-box/NC regions, involved in the control of cell proliferation [START_REF] Hughes | The 2A proteins of three diverse picornaviruses are related to each other and to the H-rev107 family of proteins involved in the control of cell proliferation[END_REF]. The absence of the conserved GxCG motif suggests that the 2A protein of SEBV1 does not have trypsin-like proteolytic activity [START_REF] Lamphear | Mapping the cleavage site in protein synthesis initiation factor eIF-4 gamma of the 2A proteases from human Coxsackievirus and rhinovirus[END_REF]. The 2C protein, like those of other parechoviruses, carries the SxxGxGKx NTase motif and the D 1271 DLxQ helicase motif (Gorbalenya et al., 1989a).

The P3 region encompasses the non-structural proteins 3A, 3B, 3C and 3D, displaying 60-62% amino acid identity with LV and nearly 42% with HPeV (Table 1 andFig 1c). The protein 3B contains the conserved tyrosine (Y) residue at the third position of the predicted N-terminus, present in all known picornaviruses and necessary for the priming function of the VPg [START_REF] Ambros | Protein is linked to the 5' end of poliovirus RNA by a phosphodiester linkage to tyrosine[END_REF]. The conserved catalytic triad H 1610 -D 1648 -C 1723 identifies the 3C protein as a cysteine-protease [START_REF] Bazan | Viral cysteine proteases are homologous to the trypsin-like family of serine proteases: structural and functional implications[END_REF] with the active site motif G 1721 xCG and G 1739 xH (Gorbalenya et al., 1989b). The 3D protein contained the five highly conserved RNA polymerase motifs: K 1930 DELR, G 2057 GxPSG, S 2061 GX 3 TX 3 N, Y 2096 GDD and F 2144 LKR [START_REF] Koonin | The phylogeny of RNA-dependent RNA polymerases of positivestrand RNA viruses[END_REF].

Phylogenetic analysis based on the most informative picornaviral protein sequence, i.e. VP1, of SEBV1 and representative members of the family Picornaviridae confirmed that this virus is related to the genus Parechovirus (Fig. 2). SEBV1 appears to be distinct from the two species constituting this genus, albeit more closely related to LV. It appears at the basal position of the LV phylogroup in the phylogenetic tree, suggesting an ancestral position. Phylogenetic analyses of 3D pol sequences (data not shown) and of the complete structural region P1 in representative picornaviral sequences gave similar results (Supplementary Fig. 1).

This study illustrates the value of deep sequencing approaches (Illumina technology) for the identification and the taxonomic classification of unclassified or misclassified viruses. This approach is useful when no genomic information is available and even when using the sample originally collected as the starting material. We developed and validated a workflow for data analysis, from raw read sequences to the generation of large contig sequences covering almost all of the genome. The aa identity between SEBV1 and LV isolates was over 55% in the polyprotein, and over 50% in P1, P2 or P3. SEBV1 exhibited the same genome organization as LV, with some common genetic features. Both SEBV1 and LV were originally isolated from rodents, suggesting that they share similar host tropism.

These various findings indicate that SEBV1 is a member of the genus Parechovirus, closely related to Ljungan virus species Despite these common features, the aa identity between SEBV1 and other known parechoviruses was below 70% for the polyprotein and in particular for the P1 protein. The G+C content of the SEBV genome differs from those of parechoviruses by more than 1%. Moreover, SEBV1 exhibits features not shared with those of LV and other parechoviruses. We therefore suggest that SEBV1 belongs to a new distinct species of parechovirus, in line with the criteria for species demarcation defined by the International Committee for Taxonomy of Viruses.

Identification of this putative third species in the genus Parechovirus highlights the viral and genetic diversity of this genus. It is interesting also to note that SEBV1 was the first LV-like virus to be isolated in Africa. However, further information is needed about this newly characterised picornavirus, especially regarding its host range, geographical distribution and genetic diversity. The tree is midpoint rooted for purposes of clarity only.

Legends to the figures

Fig. 1 .

 1 Fig. 1. Schematic representation of predicted genome organisation of SEBV1. (a) The three large sequence contigs obtained after de novo assembly and BLAST analysis are indicated with black lines. The positions of primers designed to generate overlapping PCR products and for genome resequencing are shown above the sequence contigs. (b) Representation of the genome coverage of SEBV1 after Illumina sequencing, with sequence reads mapped on the full-length genome sequence. The number of times that each nucleotide position was sequenced is indicated on the ordinate. (c) Predicted genome organization of SEBV1. Positions of predicted cleavage sites are indicated along the polyprotein, in nucleotides (upper numbers) and amino acids (lower numbers). Conserved picornaviral amino acid motifs are presented below the bar.

Fig. 2 .

 2 Fig. 2. Phylogenetic analysis based on the SEBV1 VP1 protein (310 aa) with representative members of the genera of the family Picornaviridae, and also some unassigned picornaviruses. The scale bar indicates branch length, and bootstrap values ≥ 60% are shown next to the relevant nodes. The tree is midpoint rooted for purposes of clarity only.

Fig. S1 .

 S1 Fig. S1. Phylogenetic analysis based on the P1 region (812 aa) of SEBV1 and representative members of selected genera in the family Picornaviridae, with some additional unassigned picornaviruses. The scale bar indicates branch length, and bootstrap values ≥ 70% are shown next to the relevant nodes.

Table 1 .

 1 Pairwise amino acid identities between the predicted proteins of SEBV1 and related members of the genus

	374 375 376	Parechovirus.
		Genomic features				Pairwise amino acid identity (%) (length of protein in aa)
			SBKV				Related parechoviruses
					LV1	LV1	LV2	LV3	LV4	HPeV-1	HPeV-6
	Protein		Position	Length	87-012	174F	145SLG	M1146	64-7855	Harris	NII561-2000
					(EF202833)	(AF327921)	(FJ384560)	(AF538689)	(EU854568)	(L02971)	(AB252582)
	VP0		1 Met -259 Gln		65.3 (259)	65.6 (259)	63.3 (259)	65.3 (259)	64.9 (259)	37.7 (289)	37.7 (289)
	VP3		260 Ala -502 Gln		61.1 (244)	61.5 (244)	62.3 (244)	60.7 (244)	62.3 (244)	48.8 (253)	48.0 (252)
	VP1		503 Gly -812 Gly		46.3 (317)	46.3 (317)	46.6 (317)	44.6 (319)	45.8 (319)	25.3 (234)	26.6 (234)
	P1		1 Met -812 Gly		56.6 (820)	56.9 (820)	56.5 (820)	55.8 (822)	56.7 (822)	36.7 (776)	36.5 (775)
	2A		813 Gly -946 Gln		72.6 (135)	72.6 (135)	71.9 (135)	69.6 (135)	70.4 (135)	42.6 (147)	42.7 (150)
	2B		947 Ser -1083 Gln		72.9 (138)	73.6 (138)	75.7 (140)	74.3 (140)	74.3 (140)	47.9 (122)	48.6 (122)
	2C		1084 Gly -1416 Glu		69.5 (333)	68.6 (333)	69.8 (334)	68.9 (334)	70.0 (334)	49.9 (329)	49.6 (329)
	P2		813 Pro -1416 Glu		70.9 (606)	70.6 (606)	71.6 (609)	70.4 (609)	69.8 (609)	48.2 (598)	48.1 (601)
	3A		1417 Ser -1540 Glu		43.9 (130)	43.1 (130)	46.2 (130)	46.2 (128)	44.6 (128)	23.2 (117)	24.0 (117)
	3B		1541 Arg -1569 Glu	29	37.9 (29)	37.9 (29)	37.9 (29)	41.4 (29)	37.9 (29)	20.7 (20)	20.7 (20)
	3C		1570 Ala -1763 Glu		61.1 (198)	61.6 (198)	61.1 (198)	62.8 (196)	62.2 (196)	38.9 (200)	38.9 (200)
	3D		1764 Gly -2233 Glu		67.5 (470)	67.7 (470)	67.7 (470)	66.0 (470)	65.3 (470)	47.9 (469)	48.3 (469)
	P3		1417 Ser -2233 Glu		61.2 (827)	61.3 (827)	61.7 (827)	61.0 (823)	60.4 (823)	42.4 (806)	42.0 (806)
	Polyprotein	1 Met -2233 Glu	2233	56.6 (2253)	56.9 (2253)	56.5 (2256)	55.7 (2254)	56.5 (2254)	39.8 (2180)	40.2 (2182)
	377 378

Acknowledgements

We are grateful to Laurence Ma (Plate-forme Génomique) for her excellent technical expertise and Nick J. Knowles for information about new parechovirus species.

We are grateful for the financial support of the Agence Nationale pour la Recherche (ANR 09 MIEN 019) and the Fondation pour la Recherche Médicale (FRM DMI20091117313).


M B Doost

Resonant state expansion for Transverse Electric modes of two-dimensional open optical systems

Keywords: numbers: 03.50.De, 42.25.-p, 03.65.Nk

The resonant state expansion (RSE), a rigorous perturbative method in electrodynamics, is formulated for Transverse Electrodynamic modes of an effectively 2-dimensional system. The RSE is a perturbation theory based on the Lippmann Schwinger Green's function equation and requires knowledge of the Green's function of the unperturbed system constructed from Resonant-states. I use the analytic Green's function for the magnetic field to normalize the modes appearing in the corresponding spectral Green's function. This use of the residue and cut of the analytic Green's function is a solution to the failure of the flux volume integral method of normalization for continuum states, a problem which is discussed in detail in this manuscript. In brief, the flux volume integral fails for continuum states because they are not true resonance. An analytic relation between normalized magnetic and electric modes is developed as part of the solution to these difficulties. The complex eigenfrequencies of modes are calculated using the RSE for the case of a homogeneous perturbation.

I. INTRODUCTION

Dielectric microcavities have attracted significant interest since 1989 when they were found to support whispering gallery modes (WGMs) [1], particularly with regards to the development of hybrid optoelectronic circuits [2]. Perturbed microcavities have proved to be the most promising direction of research since their emissions are not isotropic [3][4][5][6][7][8][9]. Other uses of 2-dimensional optical microcavities include miniature lasing devices [START_REF] Vahala | Optical microcavities[END_REF].

Recently there has been developed the Resonant-stateexpansion (RSE) a theory for calculating the resonances of open optical systems [15][16][17][18][19][20]. That the RSE is suitable for use with the RSE Born approximation, to calculate the far field scattering, has been demonstrated in Ref. [21][22][23]. The RSE Born approximation is a rigorous perturbation method for electrodynamic scattering, which becomes the exact solution of Maxwell's wave equations in the far field, in the limit of an infinite number of resonance being taken into account in the Born approximation expansion.

In this paper I report extending the RSE to the Transverse Electric (TE) modes of effectively two-dimensional (2D) systems (i.e. 3D systems translational invariant in one dimension) which are not reducible to effective 1D systems. I only treat systems with zero wavevector along the translational invariant direction at present. I use a dielectric cylinder with uniform dielectric constant in vacuum as unperturbed system and calculate perturbed Resonant-states (RSs) for homogeneous perturbations in order to demonstrate the method.

The treatment in this manuscript is restricted to 2dimensional systems, however for the modes of a cavity with mode wavelength only a small fraction of the thickness of the cavity it is possible to make the mathematical approximation to an effectively 2-dimensional system [START_REF] Smotrova | [END_REF]12]. A consequence of this approximation for the mathematical treatment is that the refractive index n r has to be replace with an appropriately calculated n ef f which takes into account the thickness of the resonator, for the detailed mathematics, see Ref. [START_REF] Smotrova | [END_REF]12].

In this manuscript dispersion in the refractive index is not treated directly, however it is straight-forward to add dispersion using the method from Ref. [20], which I developed in collaboration with E. A. Muljarov, based on my successful initial independent work incorporating dispersion into the RSE.

Interestingly the RSE is a near identical translation to electrodynamics of a much earlier theory from Quantum Mechanics by More, Gerjuoy, Bang, Gareev, Gizzatkulov and Goncharov Ref. [13,14]. The only difference between the two approaches is the choice of RS normalisation method. The RSE [13,14] is a perturbation theory based on the Lippmann Schwinger Green's function equation and requires knowledge of the Green's function (GF) of the unperturbed system constructed from RSs. I was able to show in Ref. [21] that the general normalisation of RSs which I derived in Ref. [17] is the most numerically stable available normalisation method. The general normalisation which I derived in Ref. [17] is based on a prototype normalisation which appeared in Ref. [19].

One difficulty in formulating an RSE for 2D systems, is the presence of a one-dimensional continuum in the manifold of RSs. This continuum is specific to 2D systems and is required for the completeness of the basis and thus for the accuracy of the RSE applied in 2D, as discussed when similar observations were previously made for Transverse Magnetic (TM) modes of a microcylinder in my Ref. [16].

In Ref. [16,17] it was shown that the continuum states in the RSE cannot be normalised by the flux volume nor-malisation,

1 + δ kn,0 = A E n (ṙ) • ∂ω 2 ε(ṙ, ω) ∂ω 2 ω=ωn E n (ṙ)dṙ (1) + lim k→kn LA E n • ∇E -E • ∇E n k 2 -k 2
n dL I will discuss the reason for this in more detail in Sec. III. Instead of using Eq. ( 1), it is necessary to make a comparison of the residue and contour integral between the spectral GF Ĝk (r,

r ′ ) = n E n (r) ⊗ E n (r ′ ) 2k(k -k n ) ≈ n φ nE n(r) ⊗ E n(r ′ ) 2k(k -k n)
(2) and the analytically derived GF. The electrodynamic wave function E n is the solution to

∇ × ∇ × E n(r) = k 2 n ε(r)E n(r) (3)
with outgoing boundary conditions (BCs) and Ĝk (r, r ′ ) is the corresponding GF for Eq. ( 3). Unfortunately it is not possible to derive the analytic GF for the electric E-field since for 2-dimensional TE modes of an open system Maxwell's wave equation is not reducible to one dimension whilst in the electric field description. However the magnetic H-field GF is available and in Appendix A I use physical arguments to construct an analytic formula linking normalised magnetic modes H n to normalised electric modes E n. The analytic Green's function for the homogeneous 2-dimensional microcylinder is given in Appendix B.

This manuscript is organised as follows Sec. II outlines the derivation of the RSE perturbation theory, Sec. III explains the special provisions made for the normalisation of the E n modes, Sec. IV gives the basis modes of the 2-dimensional RSE, Sec. V contains the numerical validation of the rigorous analytics provided in this manuscript, Appendix A and B derive the normalisation of the TE E n modes, Appendix C gives perturbation matrix elements required to reproduce the numerical results of Sec. V. Appendix D gives the evaluation of an important expression in Appendix A. Appendix E gives a new and valuable numerical recipe used in the production of numerical results for this manuscript.

II. THE RSE PERTURBATION METHOD

In this section I give a brief account of the derivation of the RSE [13,14]. The method of derivation is built upon the Lippmann Schwinger Green's function equation. For a more detailed account of these derivations I direct the interested reader to Ref. [13][14][15][16][17][18][19][20][21][22][23].

The spectral representation of the 2D system GF is modified to include this cut contribution as integral [16] Ĝk

(r, r ′ ) = n E n (r) ⊗ E n (r ′ ) 2k(k -k n ) ≈ n φ nE n(r) ⊗ E n(r ′ ) 2k(k -k n) (4)
In practice, the continuum of non-resonant states are discretised and included as cut poles in the perturbation basis. This problem is treated in depth in Appendix A, B and C. The combined index n is used to denote both real poles k n and cut poles k α simultaneously. The weighting factors φ n are defined as follows φ n = φ n = 1 for real poles, φ α for cut poles (5) After modifying the expansion of the perturbed wave functions to include the cut poles to,

E ν (r) = n φ nb nν E n(r) (6)
which is the solution to the equation,

∇ × ∇ × E ν (r) = κ 2 ν [ε(r) + ∆ε(r)] E ν (r) (7)
with outgoing BCs, I can repeat the derivation of Ref. [16] and make use of the Lippmann Schwinger Green's function equation

E ν (r) = -κ 2 ν Ĝk (r, r ′ )∆ε(r ′ )E ν (r ′ )dr ′ (8) in conjuncture with the spectral GF Ĝk (r, r ′ ) = n E n (r) ⊗ E n (r ′ ) 2k(k -k n ) ≈ n φ nE n(r) ⊗ E n(r ′ ) 2k(k -k n) (9
) to arrive at the modified the RSE Matrix equations,

κ ν n′ δ nn ′ + φ nφ ′ nV nn ′ /2 b n′ ν = k n b nν (10)
and

n′   δ nn ′ k n′ + V nn ′ 2 φ nφ ′ n k nk n′   c n′ ν = 1 κ ν c nν (11)
As always V nn ′ is given by

V nn ′ = E n(r)∆ε(r)E n′ (r) dr (12)
The perturbation method which I have briefly outlined was recently translated from Quantum Mechanics [13,14] and has now become known as the Resonant-stateexpansion.

III. MAXWELL'S WAVE EQUATION FOR MAGNETIC FIELD

In this section I discuss the numerical and analytic difficulties associated with constructing an electrodynamic RSE perturbation theory for magnetic resonances.

If we examine Maxwell's wave equation for magnetic field H(r), -∇×∇×H(r)+ε(r) -1 ∇ε(r)×(∇×H(r))+k 2 ε(r)H(r) = ∇ × J(r) -ε(r) -1 ∇ε(r) × J(r) (13) it would at first sight seem of no use to the RSE. In Eq. ( 13) we can see ∇ε(r) which has prevented me from formulating an equation analogous to Eq. ( 8) relating the perturbed eigenmodes of magnetic field to their unperturbed GF and perturbation in the dielectric profile. Therefore as I am considering perturbations in the dielectric profile I am forced to formulate the RSE in terms of electric modes E n(r).

However even in light of the previous paragraph I do need to consider Eq. ( 13) in this paper because the cut poles cannot be normalised using [17,22]

1 + δ kn,0 = A E n (ṙ) • ∂ω 2 ε(ṙ, ω) ∂ω 2 ω=ωn E n (ṙ)dṙ(14) + lim k→kn LA E n • ∇E -E • ∇E n k 2 -k 2 n dL
The method of normalisation will thus make use of Eq. ( 13). The reason I can't use Eq. ( 14) for the cut poles is that they are not true resonances of the system. This problem can be intuitively understood if we assume we can normalise a single cut pole with Eq. ( 14) and observe the resulting logical contradiction. If the discretisation of the cut is made finer by increasing the number of cut poles by several orders of magnitude, logically the normalisation of the chosen cut pole should be drastically reduced as it's weight is further shared between many of these extra poles. This final remark gives the contradiction because the normalisation calculated from Eq. ( 14) cannot change, it is fixed.

In light of the problems normalising cut poles I am forced to normalise the basis modes in this chapter by comparing spectral GFs of the form [16] Ĝk (r,

r ′ ) = n E n (r) ⊗ E n (r ′ ) 2k(k -k n ) ≈ n φ nE n(r) ⊗ E n(r ′ ) 2k(k -k n) (15
) which will still hold true for cut poles, with the analytic GF of magnetic field calculated in Appendix B. Unlike the TM case where the analytic GF for the electric field ĜE k (r, r ′ ) is available, I can only derive the analytic GF for the magnetic field ĜH k (r, r ′ ) in the case of transverse electric (TE) polarisation. The reason for this is that a scalar GF equation is available for ĜH k (r, r ′ ) but not for ĜE k (r, r ′ ) in the TE case. Fortunately in Appendix A I have been able to derive a method of expressing the electric field basis modes E n in terms of the magnetic field basis modes H n.

The Green's function for the magnetic component of an electrodynamic system is a tensor ĜH k which satisfies the outgoing wave boundary conditions and Maxwell's wave equation with a delta function source term

-∇ × ∇ × ĜH k (r, r ′ ) + ε(r) -1 ∇ε(r) × (∇ × ĜH k (r, r ′ )) + k 2 ε(r) ĜH k (r, r ′ ) = 1δ(r -r ′ ) (16)
following the derivation for the spectral representation of ĜE k in Appendix B we find ĜH k can be written as

ĜH k (r, r ′ ) = n φ nH n(r) ⊗ H n(r ′ ) 2k(k -k n) ( 17
)
where the field H n(r) satisfies,

-∇ × ∇ × H n(r) + ε(r) -1 ∇ε(r) × (∇ × H n(r)) + k 2 ε(r)H n(r) = 0 (18)
with outgoing BCs. In Appendix A I will make use of H n defined in Eq. ( 16), Eq. ( 17) and Eq. ( 18) to normalise E n. I will briefly explain the idea behind the analytics in Appendix A and direct the interested reader to that relevant appendix where I give my detailed mathematics. The E and H fields are both generated by an electric current J. By using an identical current to generate E and H fields and letting the electric current amplitude tend to zero as its harmonic wave number of oscillation tends to k n we can compare residues of ĜH k and ĜE k in order to relate normalised E n and H n. The H n modes, including the continuum modes, discussed in Appendix B, can be normalised by comparing the residue of the spectral GF and the residue of the analytic GF, derived in Appendix B.

IV. UNPERTURBED BASIS FOR 2D SYSTEMS

In this section I will use the approach outlined in Sec. III to formulate the basis states for both TE, TM polarizations. For completeness of the basis I include longitudinal electric modes which are curl free static modes satisfying Maxwell's wave equation for k n = 0. Longitudinal electric modes have zero magnetic field.

Splitting off the time dependence ∝ e -iωt of the electric fields E and D and magnetic field H, the first pair of Maxwell's equations can be written in the form

∇ × E = ikH , ∇ × H = -ikD -J (19)
where k = ω/c and D(r) = ε(r)E(r). Combining them leads to

-∇ × ∇ × E(r) + k 2 ε(r)E(r) = ikJ(r) (20)
and

∇ × ∇ × E n (r) = k 2 n ε(r)E n (r) (21)
which satisfy the outgoing wave boundary condition, and to Ĝk

(r, r ′ ) = n E n (r) ⊗ E n (r ′ ) 2k(k -k n ) ≈ n φ nE n(r) ⊗ E n(r ′ ) 2k(k -k n) (22)
for the corresponding GF.

For k z = 0, these three groups of modes of a homogeneous dielectric cylinder can be written as

TM : E = f e z , TE : H = f e z , LE : E = -∇f ,
where f (r) is a scalar function satisfying the Helmholtz equation

∇ 2 f + k 2 εf = 0 ( 23
)
with the dielectric susceptibility of the sphere in vacuum

ε(ρ, ϕ) = n 2 r for ρ R , 1 for ρ > R (24)
and

f (r) = R m (ρ, k n )χ m (ϕ) (25)
The angular parts are defined by

χ m (ϕ) =    π -1/2 sin(mϕ) if m < 0 , ( 2π
) -1/2 if m = 0 , π -1/2 cos(mϕ) if m > 0 (26)
where there is no m = 0 mode for the TE case. The χ m (ϕ) are orthonormal according to

2π 0 χ m (ϕ)χ m ′ (ϕ)dϕ = δ mm ′ (27)
The radial components satisfy,

∂ 2 ∂ρ 2 + 1 ρ ∂ ∂ρ - m 2 ρ 2 + ε(ρ)k 2 R m (ρ, k) = 0 (28)
and have the form

R m (ρ, k) = J m (n r kρ)/J m (n r kR) for ρ R , H m (kρ)/H m (kR) for ρ > R (29) in which J m (z) and H m (z) ≡ H (1)
m (z) are, respectively, the cylindrical Bessel and Hankel functions of the first kind. H m (z) is a multiple valued function, with multiple values corresponding to a range of possible boundary conditions at infinity. In order to obtain the correct boundary conditions for the RSE I provide a numerical recipe in Appendix E.

I treat in this work all polarisations, however only TM and LE polarisations mix due to the perturbations treated being strictly scalar in dielectric permittivity. If the perturbation left the dielectric profile as a tensor then all polarisations would mix.

The unperturbed RS wave functions factorise as TM :

E n (r) = A TM R m (ρ, k n )χ m (ϕ)e z TE : H n (r) = A TE m (k n )R m (ρ, k n )χ m (ϕ)e z (30)
I normalized the wave functions from the analytic TM (TE) Green's function in Appendix B with normalisation constants A(k)

A TM = 1 R 2 n 2 r -1 (31)
A TE m = 2k [J m (n r kR)] 2 (n 2 r -1) m 2 k [J m (n r kR)] 2 + R 2 k [J ′ m (n r kR)] 2
(32) The two boundary conditions at the surface of the cylinder, the continuity of the electric field and its radial derivative, produce a secular equation for the RS wave number eigenvalues k n , which has the form

D TM m (k n R) = 0 (33) D TE m (k n R) = 0 (34)
where

D TM m (z) = n r J ′ m (n r z)H m (z) -J m (n r z)H ′ m (z) (35) D TE m (z) = J ′ m (n r z)H m (z) -n r J m (n r z)H ′ m (z) (36)
and J ′ m (z) and H ′ m (z) are the derivatives of J m (z) and H m (z), respectively. Here z represents a complex argument, as opposed to the spatial coordinate used earlier.

In cylindrical coordinates, the electric vector field E(r) can be written as

E(ρ, ϕ, z) = E ρ e r + E ϕ e ϕ + E z e z =   E ρ E ϕ E z   (37)
Therefore the TM modes can be written as

E TM n (r) = A TM    0 0 R m (ρ, k n)χ m (ϕ)    (38)
Since I am considering a homogeneous micro-cylinder and TE modes, I can now calculate the normalised E n from the H n using the following simple relations derived in Appendix A

∇ × H n(r) = M kn k 2 nE n(r) (39)
which gives, inside the homogeneous cylinder only,

E TE n (r) = A TE M kn k 2 n       R m (ρ, k n) ρ dχ m (ϕ) dϕ - ∂R m (ρ, k n) ∂ρ χ m (ϕ) 0       (40)
The unknown M kn is calculated analytically in Appendix A and Appendix D.

The k = 0 LE modes are given by, inside the homogeneous cylinder only,

E LE n (r) = |m|(n 2 r -1) lim kn→0 E TE n (r) (41)
which follows from Eq. ( 14), in the limit of k = 0, noting the fast attenuation of the RS far field in this k = 0 case. For both the TM and TE case, treated here, the full GF of the homogeneous dielectric cylinder, which is defined via Maxwell's equation with a line current source term, is given by Ĝk = G k e z ⊗ e z , in which

G k (ρ, ρ ′ ) = m G m (ρ, ρ ′ ; k)χ m (ϕ)χ m (ϕ ′ ) (42)
and the radial components have the following spectral representation

G m (ρ, ρ ′ ) = n [A m (k n )] 2 R m (ρ, k n )R m (ρ ′ , k n ) 2k(k -k n ) + 0 -i∞ [A m (k)] 2 R m (ρ, k ′ )R m (ρ ′ , k ′ ) 2k(k -k ′ ) σ m (k ′ )dk ′ ≡ n [A m (k )] 2 R m (ρ, k n )R m (ρ ′ , k n ) 2k(k -k n ) (43)
derived in Appendix B. For TE modes I take A m (k) = A TE m (k) and for the TM modes I take

A m (k) = A T M m , then σ m (k) = 8J m (n r kR) 2 π 2 kR 2 D + m (kR)D - m (kR) 1 A m (k) 2 (44)
Here D ± m (z) are the two limiting values of D m (z ′ ) for z ′ approaching point z on the cut from its different sides Re z ′ ≷ 0. D ± m (z) are taken to be D TE m or D TM m in the case of TE or TM cut pole respectively.

I use the same method to discretised this continuum into cut poles k α and cut pole strengths φ (m) α as in my Ref. [18]

k (m) α = q (m) α+1 q (m) α kσ m (k)dk φ (m) α (45)
where the cut pole strength φ (m) α is defined as An example of cut poles assigned for m = 0, 10, and 20 is given in Fig. 1(a). The cut poles contribute to the RSE in the same way as the normal poles.

φ (m) α = q (m) α+1 q (m) α σ m (k)dk (46)
In the numerical calculation in Sec. V, I use N is the number of cut poles, in order to demonstrate the convergence towards the exact solution.

The next section discusses results of the RSE for different effective 2D systems. We consider a homogeneous dielectric cylinder of radius R and refractive index n r = 2 (ε = 4) with a homogeneous perturbation of the whole cylinder in Sec. V. Explicit forms of the matrix elements for these perturbations and details of their calculation are given in Appendix C.

V. HOMOGENEOUS CYLINDER PERTURBATION

In this section I present my numerical results which support the validity of the analytic results relating to the 2-dimensional TE RSE which I have provided in this manuscript.

The perturbation we consider in this section is a ho-mogeneous change of ε over the whole cylinder, given by ∆ε(ρ, ϕ) = ∆εθ(R -ρ) = ∆ε for ρ R , 0 for ρ > R (47) with the strength ∆ε = 4 used in the numerical calculation. For ϕ-independent perturbations, modes with different azimuthal number m are decoupled, and so are even and odd (cosine and sine). I show only the sine modes here, and use for illustration m = 20. The matrix elements of the perturbation are calculated analytically and given by Eq. ( 12). The homogeneous perturbation does not change the symmetry of the system, so that the perturbed modes obey the same secular equation Eq. ( 34) with the refractive index n r of the cylinder changed to n 2 r + ∆ε, and thus the perturbed wave numbers κ ν calculated using the RSE can be compared with the exact values κ (exact) ν . We choose the basis of RSs for the RSE in such a way that for the given azimuthal number m and the given number of normal RSs N we find all normal poles |k n | < k max (N ) with a suitably chosen maximum wave vector k max (N ) and then add the cut poles. We find that as we increase N , the relative error κ ν /κ (exact) ν -1 decreases as N -3 . Following the procedure described in Ref. [18] we can extrapolate the perturbed wave numbers. The resulting perturbed wave numbers are shown in Fig. 2. The perturbation is strong, creating 3 additional WGMs with m = 20 having up to 4 orders of magnitude narrower linewidths. For N = 800, the RSE reproduces about 100 modes to a relative error in the 10 -7 range, which is decreasing by one or two orders of magnitude after extrapolation. The contribution of the cut and zero frequency mode is significant: Ignoring these modes leads to a relative error of the poles in the 10 -3 range. The fact that the relative error improves by 4-5 orders of magnitude after taking into account the cut in the form of the cut poles shows the validity of the reported analytical treatment of cuts in the RSE, and the high accuracy of the discretisation method into cut poles.

VI. SUMMARY

I have applied the resonant state expansion (RSE) to effective two-dimensional (2D) open optical systems, such as dielectric micro-cylinders and micro-disks with perturbations, by finding a formula linking the normalisation of resonant states appearing in the magnetic and electric Green's functions which describe identical systems. These analytic formulas allow us to take advantage of an analytic Green's function of the magnetic field of the homogeneous micro-cylinder to normalise the modes appearing in the spectral representation of the Green's function for the electric field, electric modes which when normalised form the basis of the RSE.

Using the analytically known basis of resonant states (RSs) of an ideal homogeneous dielectric cylinder -a complete set of eigenmodes satisfying outgoing wave boundary conditions -I have treated an effectively 1D homogeneous perturbations. I investigated the convergence for this perturbations and compared the RSE with analytic solutions. I have found agreement between the RSE and analytic solutions. generated in collaboration with Dr E. A. Muljarov and presented in my PhD Thesis. I do not make use of any part of these reformulations of my original work in this manuscript, hence I am the sole author of this work.

Appendix A: Normalised electric modes expressed in terms of normalised magnetic modes

In this appendix I provide the analytic formula relating the normalised E n and normalised H n modes. The basic principals of my derivation are as follows, let us consider the case where the E and H fields are both generated by a oscillating electric current J. By using an identical current to generate E and H fields and letting the current tend to zero as its harmonic wavenumber of oscillation tends to the resonant wavenumber k n we can compare residues of ĜH k and ĜE k in order to relate normalised E n and H n. The H n modes, including the continuum modes, discussed in Appendix B, can be normalised by comparing the residue of the spectral GF and the residue of the analytic GF, derived in Appendix B. The detailed mathematics will now be thoroughly explained.

Our starting point for finding a relation between the correct normalisation of the modes appearing in the spectral Green's function for the magnetic H(r)-field and electric E(r)-field of a homogeneous dielectric microcylinder is Maxwell's equations of electrodynamics,

∇ × E = ikH , ∇ × H = -ikD -J (A1)
where we have assumed harmonic oscillation of the driving current with respect to time, so that the time derivative in Eq. (A1) can be replaced with ik.

For the E(r)-field with a current source J(r) oscillating at frequency k we have according to Maxwell,

-∇ × ∇ × E(r) + k 2 εE(r) = ikJ(r) (A2)
and also the corresponding wave equation with a source for the H(r)-field,

-∇ × ∇ × H(r) + ∇ε(r) ε(r) × (∇ × H(r)) + k 2 ε(r)H(r) = ∇ × J(r) - ∇ε(r) ε(r) × J(r) (A3)
The GF for Eq. (A2) is,

-∇× ∇× ĜE k (r, r ′ )+ k 2 ε(r) ĜE k (r, r ′ ) = 1δ(r-r ′ ) (A4)
and for Eq. (A3)

-∇ × ∇ × ĜH k (r, r ′ ) + ε(r) -1 ∇ε(r) × (∇ × ĜH k (r, r ′ )) + k 2 ε(r) ĜH k (r, r ′ ) = 1δ(r -r ′ ) (A5)
The formulas for the E(r)-field and H(r)-field in terms of their respective Green's functions are given as the following convolutions with the source,

E(r) = ĜE k (r, r ′ )ikJ(r)dr ′ (A6) H(r) = ĜH k (r, r ′ ) ∇ × J(r ′ ) - ∇ε(r ′ ) ε(r ′ ) × J(r ′ ) dr ′
(A7) Since at this stage J(r) is arbitrary we can set J(r) = 0 when ρ > R -where R -is some radius R -< R. R is the radius of the micro-cylinder. Also since we choose our basis system to be a homogeneous micro-cylinder, importantly for Eq. (A11) ∇ε(ρ, ϕ) = 0 when ρ < R. At this stage k is also arbitrary so let it be k = k n + ∆k, ∆k → 0. To simplify our discussion let W(r) be a function related to J(r) by,

ikJ(r) = 2k n∆kW(r) (A8)
Eq. (A8) allows J(r) to meet the requirements that I have just set out.

In light of Eq. (A8) Eq. (A6) becomes

lim ∆k→0 R- m φ mE m(r) ⊗ E m(r ′ ) 2k m(k n -k m + ∆k) 2k n∆kW(r ′ )dr ′ = φ nE n(r) E n(r ′ ) • W(r ′ )dr ′ = E(r) (A9)
where we have replaced ĜE k (r, r ′ ) with its spectral representation Ĝk (r,

r ′ ) = n E n (r) ⊗ E n (r ′ ) 2k(k -k n ) ≈ n φ nE n(r) ⊗ E n(r ′ ) 2k(k -k n) (A10) Similarly Eq. (A7) becomes H(r) = lim ∆k→0 R- m φ mH m(r) ⊗ H m(r ′ ) 2k m(k n -k m + ∆k) • 2k n∆k∇ × W(r ′ ) ik dr ′ = φ nH n(r) R- H n(r ′ ) • (∇ × W(r ′ )) ik n dr ′ (A11)
where we have replaced ĜH

k (r, r ′ ) with its spectral rep- resentation ĜH k (r, r ′ ) = n φ nH n(r) ⊗ H n(r ′ ) 2k(k -k n) (A12)
Eq. (A11) and Eq. (A9) can be substituted in Eq. (A1) to give,

[∇ × E n(r)] R- E n(r ′ ) • W(r ′ )dr ′ = H n(r) R- H n(r ′ ) • (∇ × W(r ′ ))dr ′ (A13)
We now make the substitution W(r) = Γ(ρ)E n(r), where,

Γ(ρ) = 1 for ρ R -, 0 for ρ > R - (A14)
We can calculate the curl of W(r) analytically. First note for a given mode we can evaluate Eq. (A1) to be equivalently (in this case)

∇ × E n(r) = n 2 r H n(r) M kn (A15)
where M kn is an unknown to be found. Using Eq. (A15) we can write the following identity for W(r):

∇ × W(r) = ∇Γ(ρ) × E n(r) + Γ(ρ)∇ × E n(r) = -δ(ρ -R -)e ρ × E n(r) + Γ(ρ) n 2 r H n(r) M kn (A16)
Substituting Eq. (A16) into Eq. (A13) and cancelling common factors of H n(r) gives

n 2 r M kn R- E n(r ′ ) 2 dr ′ = - R- H n(r ′ ) • δ(ρ ′ -R -)e ρ × E n(r ′ )dr ′ + + R- Γ(ρ ′ ) n 2 r [H n(r ′ )] 2 M kn dr ′ (A17)
Substituting Eq. (40) into Eq. (A17)

R-

E n(r ′ ) 2 dr ′ = π -π R -dϕχ m (ϕ)χ m (ϕ)• A 2 nJ m (nk nR -) k 2 nn 2 r ∂J m (nk nR -) ∂ρ + R- H n(r ′ ) 2 dr ′ (A18)
Eq. ( 40) was generated by Eq. (A21), where M kn is the unknown we are trying to find. In Eq. (A18) we took

H n(r) = A nJ m (nk nρ)χ m (ϕ)e z (A19)
and re-calculated E n according to Eq. (A21).

Calculating the first integral in Eq. (A18) by making substitutions from Eq. ( 40)

R- E n(r ′ ) 2 dr ′ = A(k n) 2 M 2 kn k 4 n R- m 2 J 2 m (nk nρ ′ ) ρ ′ + ρ ′ ∂J m (nk nρ ′ ) ∂ρ 2 dρ ′ (A20)
having made use of

∇ × H n(r) = M kn k 2 nE n(r) (A21)
which comes from taking the curl of Eq. (A15) and using the definition of RSs, specifically, inside the homogenous cylinder,

∇ × ∇ × E n(r) = k 2 nn 2 r E n(r) (A22)
Eq. (A20) is not instantly recognisable as a standard integral, however if we make use of the two following general properties of Bessel function,

J m-1 (z) -J m+1 (z) = 2J ′ m (z) (A23) J m-1 (z) + J m+1 (z) = 2m z J m (z) (A24)
we see that

R- E n(r ′ ) 2 dr ′ = A(k n) 2 n 2 r 2M 2 kn k 2 n R- J 2 m-1 (n r k nρ ′ ) + J 2 m+1 (n r k nρ ′ ) ρ ′ dρ ′ (A25)
which is a standard integral given in Appendix C.

In the last integral in Eq. (A18) we can make substitutions from Eq. ( 23) and Eq. ( 30) to obtain the following standard integral given in Appendix C, R-

H n(r ′ ) 2 dr ′ = A(k n) 2 R- J 2 m (nk nρ ′ )ρ ′ dρ ′ (A26)
Finally we use Eq. (A25) and Eq. (A26) to solve Eq. (A18) for M kn

2M 2 kn n 4 r = = J 2 m-1 (n r k nρ ′ ) + J 2 m+1 (n r k nρ ′ ) ρ ′ dρ ′ R -J m (n r k nR -) ∂J m (n r k nR -) ∂r + k 2 nε J 2 m (nk nρ ′ )ρ ′ dρ ′ = 2 n 2 r k 2 n (A27) or M 2 kn = n 2 r k 2 n (A28)
I evaluate analytically the integral in Eq. (A27) in Appendix D. Hence our analytics in this section have reached their conclusion. We can now mathematically relate normalised E n modes to normalised H n modes with the analytic formula,

∇ × H n(r) = M kn k 2 nE n(r) (A29)
I make use of this valuable derived Eq. (A29) in Appendix B, Sec. III, Sec. IV and Sec. V, to translate the correctly normalised H n modes therein derived, into correctly normalised E n modes.

Appendix B: Green's function of a homogeneous cylinder

In this appendix I nearly exactly repeat the derivation I made with E. A. Muljarov in Ref. [16] and derive an analytic GF for 2-dimensional cylindrical homogeneous electromagnetic resonators.

The only difference between the analytics here and the analytics of Appendix B Ref. [16] is that here I am restricted to the magnetic H-field GF of the 2-dimensional TE modes of the cylinder, while in Appendix B of Ref. [16] we were restricted to the E-field TM modes of a mathematically identical system.

The reason for the difference in approach between this appendix and Ref. [16] is that to formulate the E-field TE modes into effectively 1-dimensional equations is currently an unsolved problem. To generate the analytic GF of the electrodynamic problem, using the methods in this appendix, requires the recasting of the equation into a set of decoupled 1-dimensional problems.

The TE component of the GF of a homogeneous cylinder in vacuum satisfies the following equation for magnetic field

-∇ × ∇ × ĜH k (r, r ′ ) + ∇ε(r) ε(r) × (∇ × ĜH k (r, r ′ )) + k 2 ε(r) ĜH k (r, r ′ ) = 1δ(r -r ′ ) (B1) with ε(ρ) = n 2 r for ρ R , 1 for ρ > R (B2)
Using the angular basis Eq. ( 26) the GF can be written as

G H k (ρ, ρ ′ ) = 1 √ ρρ ′ m GH m (ρ, ρ ′ ; k)χ m (ϕ)χ m (ϕ ′ ) (B3)
similar to Eq. ( 42). Note that we redefined here the radial part as GH m (ρ,

ρ ′ ; k) = √ ρρ ′ G H m (ρ, ρ ′ ; k) which satisfies d 2 dρ 2 - m 2 -1/4 ρ 2 + 1 ε(ρ) dε(ρ) dρ 1 2ρ GH m (ρ, ρ ′ ; k)+ + k 2 ε(ρ) - 1 ε(ρ) dε(ρ) dρ d dρ GH m (ρ, ρ ′ ; k) = δ(ρ -ρ ′ ) (B4)
Using two linearly independent solutions f m (ρ) and g m (ρ) of the corresponding homogeneous equation which satisfy the asymptotic boundary conditions

f m (ρ) ∝ ρ m+1/2 for ρ → 0 , g m (ρ) ∝ e ikρ for ρ → ∞
the GF can be expressed as

GH m (ρ, ρ ′ ; k) = f m (ρ < )g m (ρ > ) W (f m , g m ) (B5)
in which ρ < = min{ρ, ρ ′ }, ρ > = max{ρ, ρ ′ }, and the Wronskian W (f, g) = f g ′ -f ′ g . For TE polarization, a suitable pair of solutions is given by

f m (ρ) = √ ρ • J m (n r ρk) , ρ R , a m J m (ρk) + b m H m (ρk) , ρ > R g m (ρ) = √ ρ • c m J m (n r ρk) + d m H m (n r ρk) , ρ R , H m (ρk) , ρ > R
where

a m (k) = H m (x)J ′ m (n r x) -n r H ′ m (x)J m (n r x) πix/2n r , b m (k) = n r J ′ m (x)J m (n r x) -J m (x)J ′ m (n r x) πix/2n r , c m (k) = n r H ′ m (x)H m (n r x) -H m (x)H ′ m (n r x) πin r x/2 , d m (k) = H m (x)J ′ m (n r x) -n r H ′ m (x)J m (n r x) πin r x/2
with x = kR. The Wronskian is calculated to be

W (f m , g m ) = 2in 2 r a m (k) = -nxD m (x)/π, ( B6
)
with D m (x) defined in Eq. ( 36). Inside the cylinder, the GF then takes the form

GH m (ρ, ρ ′ ; k) = π 2i ρρ ′ J m (n r kρ < )H m (n r kρ > ) - c m (k) a m (k) J m (n r kρ < )J m (n r kρ > ) (B7)
The GF has simple poles k n in the complex k-plane which are the wave vectors of RSs, given by a m (k n ) = 0, an equation equivalent to Eq. (34). The residues Res n of the GF at these poles are calculated using

r m (k n ) = - c m (k) d dk a m (k) k=kn = 2i π(n 2 r -1) m 2 k n [J m (nk n R)] 2 + R 2 k n [J ′ m (nk n R)] 2 (B8)
which can be seen by noting that inside the resonantor Eq. (B5) and Eq. (A12) are equal as k → k n as they are equal for all k, at every point in space inside the resonator.

In addition to the poles, the GF has a cut in the complex k-plane along the negative imaginary half-axis. The cut is due to the Hankel function H m (z).

Owing to the cut of the Hankel function H m (z) the GF also has a cut along the negative imaginary halfaxis in the complex k-plane, so that on both sides of the cut Gm takes different values: G+ m one the right-hand side and Gm one the left-hand side of the cut. The step ∆ Gm = G+ m -Gm over the cut can be calculated using the corresponding difference in the Hankel function:

∆H m (z) = H + m (z) -H - m (z) = 4J m (z)
The result is

∆ GH m (ρ, ρ ′ ; k) = π 2i ρρ ′ J m (n r kρ < )J m (n r kρ > )∆Q m (k) (B9) where ∆Q m (k) = 4 - c + m a + m + c - m a - m = - 4 πkR 2 1 D + m (k)D - m (k) (B10) with D m (k) given by Eq. (36). GH m (ρ, ρ ′ ; k ′ ) k -k ′ dk ′ = -i∞ 0 GH + m dk ′ k -k ′ + 0 -i∞ GH - m dk ′ k -k ′ = 2πi n Res n k -k n -2πi GH m (ρ, ρ ′ ; k) (B11)
Note that in second part of the above equation we have made use of the residue theorem, expressing the closedloop integral in the left-hand side in terms of a sum over residues at all poles inside the contour. Using Eq. (B11) the GF can be expressed as

GH m (ρ, ρ ′ ; k) = n Res n k -k n + 1 2πi 0 -i∞ ∆ GH m (ρ, ρ ′ ; k ′ )dk ′ k -k ′ (
B12) which is a generalization of the Mittag-Leffler theorem.

In conclusion then, the residues Res n of the GF contributing to Eq. (B12) are calculated as

Res n = π 2i ρρ ′ J m (n r kρ < )J m (n r kρ > )r m (k n ) (B13)
with r m (k n ) found in Eq. (B8). Given that the spatial dependence of the GF, as described by Eqs. (B12), (B13), and (B9), is represented by products of the RS wave functions R m (ρ, k n ) and their analytic continuations R m (ρ, k) with k-values taken on the cut, we arrive at the GF in the form of Eqs.(43) and (44) which are then used in the RSE, once the modes have been converted into E n form using the analytics in Appendix A.

In order to evaluate Eq. (D2) let us consider the numerator in the above integral, taking into account standard integrals,

2 R 0 J 2 m-1 (n r k nρ ′ ) + J 2 m+1 (n r k nρ ′ ) ρ ′ dρ ′ = R 2 - m -1 k n n r 2 J 2 m-1 (n r k nR)+R 2 J ′ m-1 (n r k nR) 2 + R 2 - m + 1 k n n r 2 J 2 m-1 (n r k nR)+R 2 J ′ m+1 (n r k nR) 2 (D2)
with rearrangement and use of the following identities

J m-1 (z) -J m+1 (z) = 2J ′ m (z) (D3) J m-1 (z) + J m+1 (z) = 2m z J m (z) (D4) we arrive at, R 0 J 2 m-1 (n r k nρ ′ ) + J 2 m+1 (n r k nρ ′ ) ρ ′ dρ ′ = R 2 2 J 2 m-1 (n r k nR) + J 2 m+1 (n r k nR) - R 2 J m (n r k nR) n r k nR [(m -1) J m-1 (n r k nR) + (m + 1) J m+1 (n r k nR)] + J 2 m (n r k nR)R 2 (D5)
Further use of the identities Eq. (D3), Eq. (D4) and rearrangement gives

R 0 J 2 m-1 (n r k nρ ′ ) + J 2 m+1 (n r k nρ ′ ) ρ ′ dρ ′ = R 2 2m J 2 m-1 (n r k nR) -J 2 m+1 (n r k nR) -R 2 J m+1 (n r k nR)J m-1 (n r k nR) + (n r k nR) 2 4m 2 [J m-1 (n r k nR) + J m+1 (n r k nR)] 2 (D6)
Next we move on to the denominator

= RJ m (n r k nR) ∂J m (n r k nR) ∂ρ + k 2 nε R 0 J 2 m (nk nρ ′ )ρ ′ dρ ′ = (n r k nR) 2 4m J 2 m-1 (n r k nR) -J 2 m+1 (n r k nR) + (n r k nR) 4 8m 2 [J m-1 (n r k nR) + J m+1 (n r k nR)] 2 - (n r k nR) 2 2 J m-1 (n r k nR)J m+1 (n r k nR) (D7)
where again I have made use of standard integrals, the identities Eq. (D3), Eq. (D4) and rearrangement.

Hence making use of these evaluations we arrive at the important result, Where M 2 kn is the factor that links E n to H n in Eq. (A29). where Ñm (z) = z m F m (z 2 ) is a single-valued polynomial [START_REF] Gradshtein | Table of Integrals, Series, and Products[END_REF] while ln z is a multiple-valued function defined on an infinite number of Riemann sheets. However considering Eq. (E2), since we know,

ln z 2 = ln ze i2πM 2 = ln z 2 + 2πM i (E3)
where M is an integer, we can also say that,

H m (z) = [J m (z) + iN m (z)] -4M J m (z) (E4)
Hence by Eq. (E4) we can add any multiple of 4J m (z) to H m (z) and still obtain a valid solution to the following defining equation we can see from Eq. (E4) that changing the value of the integer M corresponds to altering the boundary condition for the solution of Eq. (E5) such that as M → ∞ we have an equal balance between incoming and outgoing waves, infinitely far from the cylinder. Hence in order to find the Hankel function with the correct boundary conditions we must consider the general solution to Eq. (E5) in the limit of ρ → ∞, where Eq. (E5) becomes is the sum of incoming and outgoing plane waves. We require in Eq. (E8) A = 0. We achieve this by first calculating the ratio A/B and then adding or subtract-ing multiples of 4J m (ρk) to H m (ρk) in order to obtain a solution with the correct boundary conditions. The calculation of A/B is as follows, Hence I have formulated a numerical recipe to rotate the cut and find the correct sheet of the Hankel function. The correct sheet of the Hankel function is the sheet which provides us with outgoing boundary condition for use in the calculation of the RS eigenmodes (see pairs Eq. (33) and Eq. (35), Eq. (34) and Eq. ( 36)).

FIG. 1 :

 1 FIG. 1: (a): Cut poles kα (stars) representing the cut of the GF of a homogeneous dielectric cylinder with nr = 2, in the complex wave-number plane for m = 0, 11, and 20. Normal poles kn (open squares) are also shown. (b): Cut pole density σm(k) (solid curves) and the cut pole strength φ (m) α (stars), for the same values of m.

  ) , where N (m) is the number of normal poles in the basis for the given m and N (m) c

FIG. 2 :

 2 FIG. 2: (a): Perturbed RS wave numbers for the homogeneous perturbation Eq. (47) calculated via the RSE with N = 800 (only sine modes are shown). The perturbed poles with (+) and without (×) the cut contribution are compared with the exact solution (open squares). Unperturbed wave numbers are also show (open circles with dots). Inset: Dielectric constant profile for the unperturbed and perturbed systems. (b): Relative error in the calculated perturbed wave numbers with (heptagons) and without (triangles) contribution of the cut. Relative error for a simulation including the cut and improved by extrapolation is also shown (crossed circles).

1

 1 (n r k nρ ′ ) + J 2 m+1 (n r k nρ ′ ) ρ ′ dρ ′ R -J m (n r k nR -) ∂J m (n r k nR -) ∂r + k 2 nε J 2 m (nk nρ ′ )ρ ′ dρ ′

  Appendix E: Rotating the cutIn this section I explain how to find the correct sheet of the Hankel function, this corresponds to rotating the cut in the Hankel function on the complex plane. The correct sheet of the Hankel function gives purely outgoing plane cylindrical waves, infinitely far from the resonator.First of all let us consider how the Hankel function is defined. The Hankel function can be expressed as[START_REF] Gradshtein | Table of Integrals, Series, and Products[END_REF] Hm (z) = J m (z) + iN m (z)(E1) using a multiple-valued Neumann function N m (z) = Ñm (z

ρ 2 +

 2 k 2 H m (ρk) = 0 (E5) Let us investigate the meaning of this result further. Since lim ρ→∞ J m (ρk) = 2 πkρ e iρk + e -iρk (E6)

  lim ρ→∞ H m (-ρk n ) H m (+ρk n ) = Ae +iρkn + Be -iρkn Ae -iρkn + Be +iρkn = A B (E9) noting ℑ(k n ) ≤ 0.If A/B = 0 we have the correct sheet of the Hankel function. The Eq. (E9) should be evaluated numerically on the computer for large values of ρ, convergence will be quick due to the exponential nature of the functions involved.

Acknowledgments

The work in this paper is completely my own original mathematical derivations and numerical data. This work was completed in 2013, and shortly after its completion I presented a draft of this manuscript to my then collaborator Dr E. A. Muljarov. Upon Dr E. A. Muljarov's insistence (as my PhD supervisor and manager), a reformulated version of the mathematics in Appendix A was

Appendix C: Homogeneous cylinder perturbation

In this appendix is given the perturbation matrix elements for the homogeneous perturbation of the TE modes of a homogeneous cylinder. I give the analytic expression for the perturbation matrix elements required to produce the numerical results in Sec. V, numerical results which validate the analytics which I developed for this manuscript.

The homogeneous perturbation Eq. (47) does not mix different angular momentum m-values. The matrix elements between RS with the same azimuthal number m are given by the radial overlap integrals

where,

and

yielding for identical basis states (n = n′ )

and for different basis states (n = n′ )

When k n = 0 we can use the asymptotic form of the Bessel function,

from which we see k n → 0 cancels out everywhere in

There is an important transformation to note,

Hence with the analytic integrals in this appendix we are able to completely reproduce the analytic results in Sec. V.

Appendix D: Evaluating Equation A26

In this appendix I analytically evaluate the important integral equation Eq. (A27)