Sanjiarizaha Randriamaherijaona

Jacky Raharinjatovo

Sébastien Boyer

Durability monitoring of long-lasting insecticidal (mosquito) nets (LLINs) in Madagascar: physical integrity and insecticidal activity

Keywords: Malaria, Long-lasting insecticide-treated nets, Fabric integrity, Bio-efficacy, Madagascar

Background: Long-lasting insecticidal mosquito nets (LLINs) are highly effective for malaria prevention. However, it is also clear that durability monitoring is essential to predict when, post-distribution, a net population, no longer meets minimum WHO standards and needs to be replaced. Following a national distribution campaign in 2013, we tracked two durability indicators, physical integrity and bio-efficacy at six and 12 months post-distribution. While the loss of net integrity during this period was in line with expectations for a one-year net life, bio-efficacy results suggested that nets were losing insecticidal effect faster than expected. The rate of bio-efficacy loss varied significantly between different net brands. Methods: We tested 600 randomly selected LLINs, 200 from each of three net brands. Each brand came from different eco-epidemiological zones reflecting the original distribution scheme. Fabric integrity (size and number of holes) was quantified using the proportional hole index (pHI). A subsample of the nets, 134 new nets, 150 at six months and 124 at 12 months, were then tested for bio-efficacy using the World Health Organization (WHO) recommended method.

Results: Three net types, Netprotect®, Royalsentry® and Yorkool®, were followed. After six months, 54%, 39% and 45%, respectively, showed visible loss of integrity. The median pHI by type was estimated to be one, zero and one respectively. The percentage of damaged nets increased after 12 months such that 83.5%, 74% and 68.5%, had holes. The median pHI for each brand of nets was 47.5, 47 and 23. No significant difference in the estimated pHI at either six or 12 months was observed. There was a statistically significant difference in the proportion of hole size category between the three brands (χ 2 = 15.761, df = 4, P = 0.003). In cone bio-assays, mortality of new Yorkool® nets was surprisingly low (48.6%), mortality was 90.2% and 91.3% for Netprotect® and Royalsentry® (F (2, 131) = 81.59, P < 0.0001), respectively. At 12 month use, all tested nets were below the WHO threshold for replacement.

Conclusion:

These findings suggest that there is a need for better net quality control before distribution. More frequent replacement of LLINs is probably not an option programmatically. Regardless of prior approval, LLIN durability monitoring for quality assessment as well as net loss following distribution is necessary to improve malaria control efforts.

Background

The use of long-lasting insecticide treated nets (LLINs) is a key malaria prevention measure. Unlike conventional bed nets, LLINs kill susceptible Anopheles spp. vectors. LLINs use in sub-Saharan Africa has increased by as much as 30% since 2010 [START_REF] Who | World malaria report[END_REF]. The World Health Organization (WHO) estimates that in this part of the world, 53% of the population at risk of malaria slept under an LLINs in 2015 (95% confidence interval [CI]: 50-57).

LLINs prevent human-mosquito contact [START_REF] Zaim | Safety of pyrethroid-treated mosquito nets[END_REF] by creating both a physical and an insecticidal barrier which, theoretically, will remain in effect despite repeated washing (at least 20 times), and during extended use (up to three years of use under field conditions). However, there is mounting evidence that these assumptions about the effective life of LLINs in the field are overly optimistic in some settings. Therefore, WHO has invited the National Malaria Control Programmes (NMCP) such as the Ministry of Public Health of Madagascar to assess the durability of LLINs in operational conditions. The recommended design includes three indicators: bioefficacy, a measure of insecticidal effect, fabric integrity, a measure of physical damage sustained by the net and coverage, the loss of nets from houses where they had previously been hung. Threshold values for these indicators, below which LLINs are judged to need replacement, provide a reference that can be used to assess net loss from a population of nets, viz. nets of the same type, distributed at the same time.

In Madagascar among the 23 million inhabitants, about 2 million confirmed malaria cases and 6000 deaths were reported in 2015 [START_REF] Who | World malaria report[END_REF], which occur across different ecoepidemiological zones (Additional file 1). The two coastal regions exhibit hyperendemic patterns with transmission lasting all year in the East and more than six months per year in the West. In the South, the period of transmission is short and episodic. Fringe, Central Highlands and South are prone to outbreaks. In the fringe areas, i.e. at intermediate altitudes, the transmission pattern is seasonal, lasting from November to May (rainy season). In the Central Highlands, the transmission is unstable, and episodic or epidemic [START_REF] Madagascar | République de Madagascar Enquête sur les Indicateurs du Paludisme[END_REF]. There are four principal malaria vectors: Anopheles funestus, An. gambiae (s.s.), An. arabiensis and An. mascarensis [START_REF] Chauvet | Résultats de 5 années d'étude: Contribution à l'étude biologique des vecteurs du paludisme à Madagascar[END_REF][START_REF] Fontenille | Anopheles mascarensis a new malaria vector in Madagascar?[END_REF][START_REF] Ralisoa | Biogéographie du complexe Anopheles gambiae de Madagascar, vecteur du paludisme[END_REF]. Recently, An. coustani, infected with Plasmodium spp., was identified in the Central Highlands. The infection could be considered Plasmodium vivax or P. falciparum [START_REF] Nepomichene | Malaria case in Madagascar, probable implication of a new vector. Anopheles coustani[END_REF].

The current strategy of the Malagasy NMCP is based on effective case management and vector control, using LLINs and indoor residual spraying (IRS). In 2015, more than 11 million nets were distributed in Madagascar [START_REF] Who | World malaria report[END_REF]. Thus, LLIN coverage is thought to be relatively high, approaching the universal coverage target of one net for every two people of the population at risk as recommended by the WHO. Nets fulfilling the criteria of ≥ 80% mortality in cone test were still effective as described in the World Health Organization Pesticide Evaluation Scheme guidelines (WHOPES) [START_REF]For laboratory and field-testing of long-lasting insecticidal nets[END_REF]. However, concerns about premature loss of ITN bio-efficacy exist [START_REF] Yewhalaw | Bioefficacy of selected long-lasting insecticidal nets against pyrethroid resistant Anopheles arabiensis from south-western Ethiopia[END_REF] have raised the question of whether the bioefficacy of many of these nets is adequate (given that nets with reduced fabric integrity i.e. more holes) rely on insecticidal effect as a compensatory mechanism, a problem seen elsewhere. In Uganda, ITN bio-efficacy testing indicated that Permanet® 2.0 LLINs produced 74% functional mortality after two years of household use in rural conditions [START_REF] Kilian | Long-term field performance of a polyester-based long-lasting insecticidal mosquito net in rural Uganda[END_REF]. Similarly, in Ethiopia, Permanet® 2.0 LLINs, used for two years, showed 67-72% mortality against Anopheles arabiensis [START_REF] Fettene | Utilization, retention and bio-efficacy studies of PermaNet® in selected villages in Buie and Fentalie districts of Ethiopia[END_REF]. In Cambodia, only 73% of tested LLINs fulfil the WHO criteria against An. dirus (s.s.) susceptible strain [START_REF] Van Roey | A phase III trial to evaluate the efficacy, fabric integrity and community acceptance of Netprotect® using a recommended long-lasting insecticidal net as positive control[END_REF].

To guide planning around the timing of LLIN campaigns to achieve sustainable impact, net programs support LLIN durability monitoring. In this study, we report on physical integrity and the bio-efficacy of three brands of nets distributed in 2013 in six districts of Madagascar.

Methods

Study areas

The study took place in six districts (Fig. 1). The study sites represent the environments and cultural settings in which LLINs were distributed during the mass campaign. Selected sites (Table 1) were chosen based on criteria such as epidemiology, LLIN brand distributed to local households and accessibility.

LLINs

Three LLINs brands followed during the study are described in Table 2. Netprotect® ITNs manufactured by Bestnet Europe LTD, is an LLIN made of polyethylene with deltamethrin incorporated into the 118-denier monofilament with a target dose of 1.8 g active ingredient (a.i.)/kg, corresponding to 68 mg a.i./m 2 . Royal-sentry® ITNs manufactured by Disease control Technologies LLC, is an LLIN made of polyethylene with alpha-cypermethrin incorporated into the 150denier monofilament with a target dose of 5.8 g (a.i.)/ kg corresponding to 261 mg/m 2 . Yorkool® ITNs manufactured by Yorkool International Company LTD, is an ITN made of polyester fibre coated with deltamethrin 55 mg (a.i.)/m 2 . All LLINs tested were rectangular.

Study design and sample size

A phase III field trial, designed as a prospective longitudinal study, was set up to study the performance of three brands of nets. There were two rounds of data collection during which the same LLINs were evaluated after six months and 12 months in the field (initial distribution occurred in late 2013). The 6-month survey was done from April to June 2014 for five sites (Ambanja, Morondava, Antsiranana, Mandoto and Sakaraha) and in July for Toamasina.

Toamasina was different from the other sites because continuous LLIN distribution started there in January 2014. The 12 month (in 2014) surveys were conducted from September to November for the five sites and in February for Toamasina. This study adopted WHO guidelines for the minimum sample size required per product advised for better precision. During this study, 100 LLINs per site, selected randomly according to site accessibility, were used for physical integrity assessment using the WHO recommended hole assessment method [START_REF]For laboratory and field-testing of long-lasting insecticidal nets[END_REF]. A subsample of 30 LLINs per site/time point was randomly selected for further bio-efficacy test. Nets collected for the bioassays analysis were replaced with new nets at each point in time. Collected nets were labelled and stored in separate plastic bags which were transported to a central laboratory for physical integrity and bioassay testing.

A total of 128 LLINs newly removed from plastic storage bag, including Netprotect® (n = 40), Royalsentry® (n = 46) and Yorkool® (n = 48) were used to assess and develop a baseline, for LLIN bio-efficacy profile.

LLIN physical integrity

LLINs sampled during the bio-efficacy assessment were scored for fabric-integrity (physical integrity). Integrity was quantified by as described in WHOPES [START_REF]For laboratory and field-testing of long-lasting insecticidal nets[END_REF]. Hole sizes were categorised into four groups; holes smaller than a thumb (0.5-2.0 cm), holes between a thumb and a closed fist (2-10 cm), holes between a closed fist and a head (10-25 cm) and holes bigger than a head (> 25 cm). A proportionate hole index (pHI), which characterises the midpoint diameter to an estimated hole size [START_REF] Allan | An observational study of material durability of three World Health Organization recommended long-lasting insecticidal nets in eastern Chad[END_REF][START_REF] Wills | Physical durability of PermaNet 2.0 long-lasting insecticidal nets over three to 32 months of use in Ethiopia[END_REF][START_REF] Bhatt | Effectiveness and durability of interceptor long-lasting insecticidal nets in a malaria endemic area of central India[END_REF], was calculated by making the sum of the holes weighted by size for each net. For these groups, the weights used to calculate the pHI were 1, 23, 196 and 576. To better translate the hole index to an integrity status (net condition) for each sampled net, the pHI is categorised into "good" (0-64), "damaged" (65-642), and "so torn" that protection from mosquitoes was judged to be compromised (≥ 643) [START_REF]For laboratory and field-testing of long-lasting insecticidal nets[END_REF].

ITN bio-efficacy assessment: WHO cone test method

Standard WHO cone bioassays were performed with a susceptible laboratory strain of Anopheles arabiensis Abbreviation: a.i. active ingredient [START_REF] Randriamaherijaona | Susceptibility status of Anopheles arabiensis (Diptera: Culicidae) commonly used s biological materials for evaluations of malaria vector control tools in Madagascar[END_REF], following the recommendations of the WHOPES [START_REF]For laboratory and field-testing of long-lasting insecticidal nets[END_REF]. For each LLIN, five 30″ × 30″ sub-samples were cut from the LLIN selected for testing. The subsamples were cut from the top and each of the four sides of the net. Each sub-sample was placed in an aluminium foil envelope, labelled, and kept individually in a 4 °C refrigerator before conducting the bio-assay.

For each sub-sample, four cone tests were conducted at a time following standard WHO procedure [START_REF]For laboratory and field-testing of long-lasting insecticidal nets[END_REF] (Fig. 2). Five non-blood-fed, two-to-five-day-old female An. arabiensis were introduced into each cone and exposed to LLIN samples for 3 min, before being transferred to paper cups, covered with netting, and held for 24 h at 28 °C and 80% humidity with access to 10% sugar solution. Mortality was recorded 24 h post-exposure. By following this methodology, a total of 100 mosquitoes were tested for each net. Each day of testing, four cones, each with 10 An. arabiensis were fixed on a non-impregnated net as a negative control. If the mortality in the control was < 10% for a given day, the data were adjusted with Abbott's formula [START_REF] Abbott | Method of computing the effectiveness of an insecticide[END_REF]. If the mortality in the control was > 10%, all the tests for that day were repeated. The standard protocol recommends using a mixed outcome, i.e. mortality ≥ 80% or KD ≥ 95% to consider a net as effective. Results from a study on nets bioassay suggest that mortality outcome was better than the KD outcome at predicting the validity of LLINs.

Data analysis

All data were recorded on standard forms, before being entered into an Excel database and then imported into R, version R-3.1.3, for statistical analysis [START_REF] Core | R: A language and environment for statistical computing[END_REF]. Categorical variables were compared using Chi-square test, and continuous, discrete variables were tested using Analysis of Variance (ANOVA) to test the significance among the three study areas or the three brands of nets. A P-value < 0.05 was considered significant. The Tukey's HSD test was applied to assess the significance of the differences.

Results

Physical condition of LLINs

A total of 600 nets, 200 for each of the three brands, were sampled and scored for integrity six months after distribution and hanging. There were 47%, 62% and 46% of nets examined were completely intact (no holes) for Netprotect®, Royalsentry® and Yorkool®, respectively. The percentages of intact nets dropped to 16.5% (n = 200), 26% (n = 200) and 31.5% (n = 200) for Netprotect®, Royalsentry® and Yorkool®, respectively, at 12 months post-distribution (Table 3). The median and interquartile range (IQR 0.25-0.75) was used to present the pHI in Fig. 3 and Table 4. After six months of use, the median pHI was 1 (IQR 0-66.25) for Netprotect®, 0 (IQR 0-27) for Royalsentry® and 1 (IQR 0-25) for Yorkool®. In Toamasina and Morondava, where Netprotect® was distributed and in Antsiranana and Ambanja where Royalsentry® was distributed, significant difference of pHI median values was observed (F (1,198) = 5.08, P = 0.02) (Table 5). After 12 months, the pHI increased to 47.5 (IQR 2-271.2), 47 (IQR 0-162.5) and 23 (IQR 0-123) for Netprotect®, Royalsentry® and Yorkool®, respectively. The median and interquartile range for each district are shown in Table 4.

At six months, the mean pHI for polyethylene and polyester nets was 71.6 and 101, respectively. At 12 months, this increased to 271 and 251 for polyethylene and polyester nets, respectively. There was no significant statistical difference between net threads neither at six or 12 months.

The proportion of LLINs judged to be in "good", "damaged", or "too torn" categories at different ages of follow-up are summarised in Table 6. After six months use, more than 70% of distributed nets from the three bands were in "good" condition. At 12 months post-distribution, 55.6%, 56.8% and 69.2% of Netprotect®, Royalsentry® and Yorkool®, respectively, were in "good" condition. There was a statistically significant difference in the proportion of hole size category between the three brands (χ 2 = 15.761, df = 4, P = 0.003) with Yorkool® showing greater loss of integrity followed by NetProtect® and RoyalSentry®.

Bio-efficacy

Results from WHO cone bioassays are presented in Table 7. Mortality in the negative control never exceeded 2%. At baseline, Yorkool® LLINs were already "not fully effective" according to the threshold established by WHO (mortality < 80%). There was a significant difference in mortality between the three net brands (F (2, 131) = 81.59, P < 0.0001). Significant difference between Yorkool® and Netprotect®, both of which were treated with deltamethrin, was observed (Tukey's HSD test, P < 0.0001). However, there was no statistical difference between mortality induced by Netprotect® and Royal Sentry®. New Netprotect® LLINs caused an average mortality of 91% with 90% (36/40) being above the minimum threshold. For Royalsentry®, the average induced mortality was 90%, and 91.3% (42/46) were above the threshold for acceptable insecticidal effect. Yorkool®, in contrast, presented the lowest induced mortality. The average mortality was 48.6%, with only 20.8% (10/38) being above the minimum threshold (Fig. 4). After six months, mortality decreased significantly for all three net types enrolled in the assessment (F (2, 147) = 6.33, P = 0.002). The average mortality was 37.4% for Netprotect®, 32 and 23.1%, respectively, for Royalsentry® and Yorkool ®. Only one Royalsentry® and one Netprotect® net scored above the threshold level; none of the Yorkool® nets tested was found to be above the minimum threshold. The difference between Royalsentry® and Netprotect®, both made of polyethylene, was not significant (Tukey's HSD test, P = 0.23). Even though both are treated with the same pyrethroid insecticide, deltamethrin, Yorkool® and Net-protect® LLINs caused significantly different mortality rates (Tukey's HSD test, P = 0.001). No significant difference was noted between mortality caused by Yorkool® and Royalsentry® (Tukey HSD, P = 0.07).

At 12 months, the mortality associated with all LLINs tested decreased considerably. For LLINs treated with deltamethrin (Netprotect® and Yorkool®), the average mortality rates were 11% and 14%, respectively, with no significant difference between those brand nets. Alphacypermethrin-treated LLINs (Royalsentry®) showed an average mortality rate of 23.1%. A significant difference was observed between Royalsentry® -Netprotect® and Royalsentry® -Yorkool®. However, none of the 12-month nets exceeded the minimum threshold characteristic of nets described as "in need of replacement". Differences in mortality associated with the location of the LLIN were also observed. At six months, there was a significant difference between percent mortality of Netprotect® distributed in Toamasina and Morondava (F (1, 58) = 4.18, P = 0.04), 41.9 vs 33%, respectively. In Antsiranana and Ambanja, where Royalsentry® was distributed, mortalities were 43.8 and 19.7%, respectively (F (1, 58) = 32.49, P < 0.0001). In Mandoto and Sakaraha, where Yorkool® was distributed, mortality was 25.4 and 21.1%, respectively (F (1, 28 ) = 0.83, P = 0.37) (Table 8, Fig. 5). At 12 months, the average mortality rate for all three brands ranged between 6.9 and 25.9% (Fig. 6). No significant difference was found between the values (Table 8).

Discussion

This study is one of only a few studies that examine bio-efficacy as well as fabric integrity under operational conditions in Madagascar. In the present assessment, more than half of nets were still in a good physical condition after 12 months use. This result parallels observations made in an earlier (2012) investigation involving epidemiological and entomological methods in the South-East region of Madagascar, to identify factors that could have caused a malaria outbreak. Among 39 LLINs collected during the study, 15%, 42.5% and 42.5% were observed to be "good", "damaged" and "too torn" condition, respectively, after two years of use [START_REF] Kesteman | Multiple causes of an unexpected malaria outbreak in a high-transmission area in Madagascar[END_REF].

For LLINs made of polyester, the percentage of nets having any holes after 12 months was higher than in Western Uganda where 33.7% of nets were assessed as having holes after one year [START_REF] Kilian | Evidence for a useful life of more than three years for a polyester-based long-lasting insecticidal mosquito net in western Uganda[END_REF]. In Zambia, 9.6% of polyester and polyethylene nets were classified as "too torn" after 12 months in the field [START_REF] Tan | A longitudinal study of the durability of long-lasting insecticidal nets in Zambia[END_REF], which is a relatively low proportion compared to observation in this study where 11.9% of polyester nets and 12.1% of polyethylene nets were torn. There is evidence that at the household level, LLINs can inhibit blood feeding even when they are in damaged condition (65 < pHI < 642). This is partially due to the repellent effect of pyrethroids incorporated or coated on the net [START_REF] Kilian | Evidence for a useful life of more than three years for a polyester-based long-lasting insecticidal mosquito net in western Uganda[END_REF]. Nonetheless, a high value of pHI (> 643), would easily allow for a mosquito to enter an LLIN to bite a sleeping human and then take rest outside the LLIN.

The findings on net bioassays were surprising, especially given that most programs assume that LLINs retain their insecticidal activity for at least three years [START_REF]Insecticide treated nets: a position statement[END_REF]. At the baseline of the current study, 10% of new Netprotect® nets and more than 75% of new Yorkool® nets did not meet the WHO cut-off value. The proportion of new nets which did not meet the WHO criteria was probably due to a problem in the manufacturing process. In Cambodia, 100% of Netprotect® used as baseline met the WHOPES criteria but, 43% of them had a deltamethrin content below the target dose [START_REF] Van Roey | A phase III trial to evaluate the efficacy, fabric integrity and community acceptance of Netprotectis® using a recommended long-lasting insecticidal net as positive control[END_REF]. These results suggest that cone bioassays alone might be not adequate to assess the comparative efficacy of these nets. Therefore, more elaborate tests such as a ring-net bioassay, which measures the median knockdown time (MKDT) of mosquitoes on three different LLINs should be performed to understand the bioavailability of the insecticide on the LLINs. The MKDT is expected to be directly correlated to the insecticide concentration on the surface for fast acting pyrethroids. It could be helpful in assessing the amount of insecticide left on the net thread [34]. Added to that, tunnel test was not done here due to lack of materials. Such a test is recommended by the WHOPES, for each net that fails to meet the criteria of the WHO cone test [START_REF] Madagascar | République de Madagascar Enquête sur les Indicateurs du Paludisme[END_REF]. Bioassay results from this current study showed one very low mortality rate at six and 12 months uses. At 12 months, average percent mortality was 11%, 14% and 23% for Netprotect®, Yorkool® and Royalsentry®, respectively. However, other studies in another country showed that at 12 months, the average percent of mortality could be relatively high depending on net type [START_REF] Tan | A longitudinal study of the durability of long-lasting insecticidal nets in Zambia[END_REF][START_REF] Dutta | Long-term insecticidal activity and physical integrity of Olyset nets in Tafea Province, Vanuatu[END_REF][START_REF] Anshebo | Estimation of insecticide persistence, biological activity and mosquito resistance to PermaNet® 2 long-lasting insecticidal nets over three to 32 months of use in Ethiopia[END_REF]. A significant difference of percent mortality between two localities after six months may be explained by user behavior. In Madagascar, some reasons could lead the household to frequently use or not-use of bed net, for example the feeling of suffocation during the night, the skin irritation due to insecticide on the net, the room becomes darker [START_REF] Mattern | Tazomoka is not a problem". Local perspectives on malaria, fever case management and bed net use in Madagascar[END_REF].

The highest proportion of nets that needed replacement was found in an area where Netprotect® was distributed; 52.6% (10/19) of them were too torn and whose protective bio-efficacy for the user was in serious doubt. In this study, the finding was that more than 90% and 84% of nets collected after six and 12 months were in good or damaged condition but did not meet the WHO bioassay cut-off value. These results are important to better understand the link between the standard measure of holes, by weighting the number of each hole size and insecticide bio-availability.

The results from this study could have several implications for the LLIN strategy of the NMCP. The finding was that percent mortality induced by Yorkool® was lower than those induced by Netprotect® and Royalsen-try® even if newly removed from plastic storage bag. This information may be relevant to the Malagasy NMCP regarding the effectiveness of LLINs brand, based on bio-efficacy, for future LLIN mass distribution through campaigns. Our findings are limited to the three LLIN brands, Netprotect®, Royalsentry® and Yorkool®, distributed during the mass campaign in 2013; other LLIN brands could perform differently under the same or different conditions. Hence, results from this study may not be extrapolated to other LLIN brands, even in similar settings. Moreover, more research still needs to be conducted to determine how the physical integrity and the residual bio-efficacy, of any brand, affect its skill to prevent and reduce malaria transmission. In the one hand, monitoring of new nets is needed before an LLIN mass campaign, on their arrivals in port. Before starting such mass distribution, it would be essential that nets are checked for their compliance with WHO specifications. Results from this study also show the importance of the quality control along the supply chain, right through hanging of the LLIN. On the other hand, the manufacturers need to clarify the quality-assessment/ quality-control (QA/QC) of their product and the guarantee to prospective LLIN customers of the "performance" of their products for three years. This guarantee is currently based on laboratory results from insecticidal activity and burst strength tests [START_REF]For laboratory and field-testing of long-lasting insecticidal nets[END_REF]. As more field monitoring of LLINs durability is conducted by public or private institutions, more evidence of LLIN insecticidal activity will inform these changes needed.

Conclusion

This study is the first to report on the performance of nets under operational conditions in Madagascar by checking the physical integrity and the insecticidal efficacy of new nets and nets post distribution. The findings from this study highlight the low insecticide efficacy of Yorkool® even if unused. This recommends that there is a need for better net quality control before LLIN mass deployment.

Fig. 1

 1 Fig. 1 Study sites with different brand of nets

Fig. 3

 3 Fig. 3 Median and interquartile range (IQR 0.25-0.75) of proportionate hole index (pHI)

Fig. 4

 4 Fig. 4 Bio-efficacy results on baseline nets

Fig. 6

 6 Fig. 6 Bio-efficacy results on 12 months of use, by locality. Localities: A, Toamasina (n = 17); B, Morondava (n = 21); C, Antsiranana (n = 17); D, Ambanja (n = 30); E, Mandoto (n = 18); F, Sakaraha (n = 21)

Table 1

 1 Characteristics of net brands used during the study

	Brands	Manufacturer	Thread	Insecticide	Concentration	Fabrication
	Netprotect®	Best Net Europe LTD	Polyethylene	Deltamethrin	68 mg a.i./m 2	Incorporated
	Royalsentry®	Disease control Technologies LLC	Polyethylene	Alpha-cypermethrin	261 mg a.i./m 2	Incorporated
	Yorkool®	Yorkool International Company LTD	Polyester	Deltamethrin	55 mg a.i./m 2	Coated

Table 2

 2 Study sites with epidemiological characteristics

	Characteristics	District	Netprotect®	Royalsentry®	Yorkool®
	Perennial transmission	Toamasina	100 nets	-	-
		Ambanja	-	100 nets
	Long transmission sites	Morondava	100 nets	-	-
		Antsiranana	-	100 nets	-
	Seasonal transmission	Mandoto	-	-	100 nets
		Sakaraha	-	-	100 nets

Fig. 2 Sampling locations used for a rectangular type bed net and bioassay

Table 3

 3 Percentage of unholed nets

	Follow-up time	District	6 months		12 months
			n	%	n	%
	Netprotect®	Morondava	100	33	100	19
		Toamasina	100	60	100	14
	Royalsentry®	Ambanja	100	80	100	30
		Antsiranana	100	43	100	22
	Yorkool®	Mandoto	100	47	100	28
		Sakaraha	100	44	100	35

Table 4

 4 Median and interquartile range (IQR 0.25-0.75) of proportionate hole index (pHI)

	Net brand	Months	n	Average pHI	Median pHI	Interquartile range
	Netprotect®	6	200	121.5	1	0-66.25
		12	200	281.5	47.5	2-271.2
	Royalsentry®	6	200	79.65	0	0-27
		12	200	220.9	47	0-162.5
	Yorkool®	6	200	71.64	1	0-25
		12	200	271	23	0-121

Table 5

 5 Median and interquartile range (IQR 0.25-0.75) of proportionate hole index (pHI) by locality

	Net brand	District	Months	n	Average pHI	Median pHI	Interquartile range
	Netprotect®	Toamasina	6	100	81.47	0 a	0-34
			12	100	279.6	46	2-271.8
		Morondava	6	100	163.1	23	0-77
			12	100	283.2	50	2-256.2
	Royalsentry®	Antsiranana	6	100	114.2	23 b	0-77
			12	100	242.5	51.5	4-193.5
		Ambanja	6	100	44.75	0	0
			12	100	199	28	0-127
	Yorkool®	Mandoto	6	100	54.3	1	0-25
			12	100	355.2	23	0-164.2
		Sakaraha	6	100	88.99	1	0-25
			12	100	187.7	23	0-97

a Significant difference between districts with the same brand nets (F (1, 198) = 3.45, P = 0.05) b Significant difference between districts with the same brand nets (F (1, 198) = 5.08, P = 0.02)

Table 6

 6 Physical condition of nets by locality

	LLIN physical condition			Good		Damaged		Too torn
		District	n	6 months	12 months	6 months	12 months	6 months	12 months
	Netprotect®	Toamasina	100	78	58	20	28	2	14
		Morondava	100	70	53	21	31	9	16
	Royalsentry®	Antsiranana	100	72	52	22	39	6	9
		Ambanja	100	92	61	5	29	3	9
	Yorkool®	Mandoto	100	87	68	12	18	1	14
		Sakaraha	100	79	70	20	20	1	10

Table 7

 7 Bio-efficacy results comparing three LLINs products.

	There was no significant difference between values which share
	the same letters
		Netprotect®	Royalsentry®	Yorkool®
		n	Mortality	n	Mortality	n	Mortality
	Beseline	40	91.1 a	46	90.2 a	48	48.6 b
	6 months	60	37.4 c	60	32.0 d	30	23.1 e
	12 months	38	11.0 f	47	23.1 g	39	14.0 f

Table 8

 8 Bio-efficacy of LLINs products by locality

	Follow-up time		6 months		12 months
			No. of tested nets	Average mortality (%)	No. of tested nets	Average mortality (%)
	Netprotect®	Toamasina	30	41.9 a	17	6.9
		Morondava	30	33	21	14.3
	Royalsentry®	Antsiranana	30	43.8 b	17	18.1
		Ambanja	30	19.7	30	25.9
	Yorkool®	Mandoto	14	25.4	18	16.2
		Sakaraha	16	21.1	21	12.1

a Significant difference between districts with the same brand nets (F (1, 58) = 4.18, P = 0.04) b Significant difference between districts with the same brand nets (F (1, 58) = 32.49, P < 0.0001)

Acknowledgements

We especially thank data collectors from Population Service International Madagascar (PSI Madagascar) and Institut Pasteur de Madagascar (IPM). Special thanks to the Malagasy National Malaria Control Programme (NMCP), the President's Malaria Initiative (PMI), Ray Beach from the Centers for Disease Control and Prevention (CDC) and Judith Hedge from PMI. The opinions expressed herein are those of the authors and do not necessarily reflect the views of the USAID.

Funding

This research was supported by the United States Agency for International Development (USAID) and funded by the President's Malaria Initiative (PMI).

Availability of data and materials

All data generated or analysed during this study are included in this published article and its additional file.

Additional file

Additional file 1: Malaria transmission patterns in the districts of Madagascar (TIFF 310 kb) Abbreviations ANOVA: analysis of variance; CI: confidence interval; IQR: interquartile range; ITN: insecticide treated net; KD: knock-down; LLINs: long-lasting insecticidal nets; MKDT: median knock-down time; NMCP: National Malaria Control Programme; pHI: proportional hole index; QA/QC: quality-assessment/qualitycontrol; WHO: World Health Organization; WHOPES: World Health Organization pesticide evaluation scheme Authors' contributions SR, JR and SB designed the study. SR and JR were involved in data analysis. SR wrote the first draft of the manuscript. JR and SB critically revised and substantially contributed to drafting the manuscript. All authors read and approved the final manuscript.

Ethics approval and consent to participate

This study was reviewed by both CDC Atlanta and Madagascar local ethics board and was determined quality improvement that did not need ethical clearance. However, written informed consent was sought for participants in this study.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


Songsong Ji
email: songsong.0211@163.com

Yibo Yang
email: ybyang@uchicago.edu

Gang Pang

Xavier Antoine
email: xavier.antoine@univ-lorraine.fr

Accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations on rectangular domains

Keywords: Schrödinger equation, heat equation, semi-discretization, rectangular boundary, artificial boundary conditions, Green's function 2010 MSC: 81Q05, 35Q41, 65M06

The aim of this paper is to design some accurate artificial boundary conditions for the semi-discretized linear Schrödinger and heat equations in rectangular domains. The Laplace transform in time and discrete Fourier transform in space are applied to get the Green's functions of the semi-discretized equations in unbounded domains with single-source. An algorithm is given to compute these Green's functions accurately through some recurrence relations. Furthermore, the finite-difference method is used to discretize the reduced problem with accurate boundary conditions. Numerical simulations are presented to illustrate the accuracy of our method in the case of the linear Schrödinger and heat equations. It is shown that the reflection at the corners is correctly eliminated.

Introduction

The Schrödinger and heat equations in infinite domains are standard models with many interesting applications in computational physics and engineering. The Schrödinger equation is one of the fundamental equations arising in quantum mechanics [START_REF] Schrödinger | An undulatory theory of the mechanics of atoms and molecules[END_REF], but has also many other important applications including electromagnetic/acoustic wave propagation [START_REF] Craig | Molecular Quantum Electrodynamics[END_REF] through the so-called parabolic equation, plasma physics [START_REF] Hazeltine | The Framework of Plasma Physics[END_REF], seismic migration [START_REF] Yilmaz | Seismic Data Analysis: Processing, Inversion, and Interpretation of Seismic Data[END_REF] or optical fibers with the Fresnel equation [START_REF] Senior | Optical Fiber Communications: Principles and Practice[END_REF]. The heat equation defined on an unbounded spatial domain arises as a model for fluid dynamics [START_REF] Ferziger | Computational Methods for Fluid Dynamics[END_REF] and also for financial applications [START_REF] Black | The valuation of options and corporate liability[END_REF]. The numerical solution and analysis of these linear PDEs on unbounded domains has so far received a great attention [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF]. To get a reliable computational method, a standard approach is to truncate the infinite domain around the region of interest and to implement an artificial boundary condition on the fictitious boundary (see e.g. [START_REF] Antoine | A Review of artificial boundary conditions for the Schrödinger equation[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] for Schrödinger-like equations). Alternative strategies include the Perfectly Matched Layer (PML) (and also Complex Absorbing Scaling (CAS)) [START_REF] Antoine | A Review of artificial boundary conditions for the Schrödinger equation[END_REF][START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF].

There are rather extensive studies on artificial boundary conditions for the one-dimensional case. For instance and without being exhaustive, Fevens and Jiang [START_REF] Fevens | Absorbing boundary conditions for the Schrödinger equation[END_REF] constructed some local artificial boundary conditions. In [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF][START_REF] Antoine | Unconditionally stable discretization schemes of non-reflecting boundary conditions for the one-dimensional Schrödinger equation[END_REF], the authors derive a stable quadrature rule for approximating the transparent operator at the boundary and localizations through Padé approximations of the exact operator. Wu and Zhang [START_REF] Wu | High-order local absorbing boundary conditions for heat equation in unbounded domain[END_REF] constructed some high-order local absorbing boundary conditions, and avoided high-order derivatives by using auxiliary variables on the fictitious boundary. For the theoretical analysis, Arnold et al. [START_REF] Arnold | Discrete transparent boundary conditions for the Schrödinger equation: Fast calculation, approximation, and stability[END_REF] proposed an efficient way to treat the Schrödinger equation, and proved the stability of the resulting initial boundary-value scheme. For the heat equation, Greengard and Lin [START_REF] Greengard | On the numerical solution of the heat equation in unbounded domains (Part I)[END_REF] proposed an algorithm based on the evolution of the continuous spectrum for the heat equation on unbounded domains. Wu and Sun [START_REF] Wu | Convergence of difference scheme for heat equation in unbounded domains using artificial boundary conditions[END_REF] proved the unconditional stability for the heat equation when a well-designed finite-difference scheme is used.

For the two-dimensional case, Baskakov and Popov [START_REF] Baskakov | Implementation of transparent boundaries for numerical solution of the Schrödinger equation[END_REF] obtained the exact transparent boundary conditions for the Schrödinger equation in a half-plane. Han and Huang [START_REF] Han | Exact and approximating boundary conditions for the parabolic problems on unbounded domains[END_REF][START_REF] Han | Exact artificial boundary conditions for Schrödinger equation in R 2[END_REF] constructed an exact artificial boundary condition to reduce the original heat equation to an initial boundary-value problem on a finite computational domain. They also derived the exact artificial boundary condition for the two-dimensional Schrödinger equation by expressing the solution with Hankel's functions, implemented with a cut-off. In [START_REF] Arnold | Discrete transparent boundary conditions for the Schrödinger equation on circular domains[END_REF], Arnold et al. proposed some discretized transparent boundary conditions for the time-dependent Schrödinger equation on a circular computational domain, and illustrated the accuracy, stability, and efficiency of the proposed method. Zhang et al. [START_REF] Li | Local artificial boundary conditions for Schrödinger and heat equations by using high-order azimuth derivatives on circular artificial boundary[END_REF] designed high-order artificial boundary conditions for the two-dimensional Schrödinger equation on a circular boundary by rationally approximating the kernel functions. For more general cases, Antoine et al. [START_REF] Antoine | Numerical schemes for the simulation of the two-dimensional Schrödinger equation using nonreflecting boundary conditions[END_REF] constructed and studied a Crank-Nicolson-type discretization for some increasing order approximations of the non-reflecting boundary condition. They also proposed some absorbing boundary conditions for general nonlinear and two-dimensional Schrödinger equations with an exterior potential (see [START_REF] Antoine | Absorbing boundary conditions for general nonlinear Schrödinger equations[END_REF][START_REF] Antoine | Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential part I: construction and a priori estimates[END_REF][START_REF] Antoine | Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential. Part II: discretization and numerical results[END_REF]). Their approach relies on pseudo-differential operators to construct truncated boundary conditions, which can be proved to be stable. Most of these boundary conditions are proposed on a circular or smooth boundary and need a truncation rule in the final implementation. We refer to [START_REF] Antoine | A Review of artificial boundary conditions for the Schrödinger equation[END_REF] for an extensive survey for the Schrödinger equation, and to [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] for a recent overview for absorbing boundary conditions and PMLs for the Schrödinger, Klein-Gordon and Dirac equations.

Considering a rectangular/square computational domain is natural in practice, but deriving an associated accurate artificial boundary condition is non trivial most particularly because of the presence of the singular corner points. Until now, there are only a few results for building accurate boundary conditions for the Schrödinger, wave and heat equations on a rectangular computational domain. Some contributions can be found e.g. in the references [START_REF] Engquist | Absorbing boundary conditions for the numerical simulation of waves[END_REF][START_REF] Higdon | Absorbing boundary conditions for difference approximations to the multi-dimensional wave equation[END_REF]. To deal with the Schrödinger and heat equations on unbounded domains, we propose first to introduce a second-order spatial discretization of the associated operator on a uniform rectangular grid. Inspired by our previous work in lattice dynamics [START_REF] Pang | Time history kernel functions for square lattice[END_REF], an efficient algorithm is designed to accurately compute the Green's function of the semi-discretized Schrödinger equation with a single-source. Based on this computation, we are able to design some accurate boundary conditions that eliminate the corner reflection when a rectangular boundary is considered.

The paper is organized as follows. In Section 2, the Green's functions of the semi-discretized equations on unbounded domains with a single-source are derived and an accurate algorithm is designed for computing the Green's functions. They are next used to build the accurate artificial boundary conditions for both the time-dependent Schrödinger and heat equations. In Section 3, the fully discretized schemes are presented for the Schrödinger and heat equations with the accurate artificial boundary conditions set on a rectangular domain. Numerical examples are reported in Section 4 to demonstrate the accuracy of our method. We end this paper (Section 5) with some concluding remarks and perspectives.

Building accurate artificial boundary conditions for the semi-discretized Schrödinger and heat equations

We first start in subsection 2.1 by the basic ideas to design accurate artificial boundary conditions for the linear Schrödinger equation. The heat equation is treated shortly in subsection 2.2.

The case of the linear Schrödinger equation

To design some exact artificial boundary conditions, let us rewrite the Schrödinger equation as:

find u := u(x, t) such that i ∂ ∂t u = - ∂ 2 ∂x 2 + ∂ 2 ∂y 2 u, (x, t) in R 2 ×]0, T ], (1)
u(x, t = 0) = 0, x ∈ Ω ext = R 2 \Ω int , (2)
u(x, •) → 0, for ||x|| → +∞. (3)
The initial data u 0 (x) is supposed to be compactly supported into the interior square domain Ω int =]0, a[ 2 (the extension to a rectangular domain is straightforward up to some additional notations), Ω ext being the associated exterior domain.

The maximal time of computation is T . The complex number i is defined by i := √ -1. Let us remark that adding a purely time-dependent potential can be directly considered by eliminating the potential effect through a gauge change [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF].

Let us introduce ∆x as the uniform mesh size in the x-and y-directions to discretize Ω int with (L + 1) segments in each direction, i.e. (L + 1)∆x = a. After scaling the time by (∆x) 2 , the spatial semi-discretization of (1) with a five-points stencil finite difference scheme gives

iu m,n = -(u m+1,n + u m-1,n + u m,n+1 + u m,n-1 -4u m,n ), in Z 2 ×]0, T ], (4)
u m,n | t=0 = 0, in (m, n) ∈ Z 2 \{(p, q)/1 ≤ p ≤ L, 1 ≤ q ≤ L}, (5)
u m,n → 0, when |m|, |n| → +∞, (6)
where u m,n (t) approximates u(m∆x, n∆x, t), and um,n is the time derivative of u m,n . The computation of ( 4)-( 6) should be confined within a finite computational domain {(m, n)/1 ≤ m ≤ L, 1 ≤ n ≤ L} where the value of u m,n is needed for m = 0, L + 1, and n = 0, L + 1, to complete system (4) through a boundary relation.

Here we derive the Green's function f m,n (t) of the semi-discretized linear Schrödinger equation with single-source. Let us consider the following system

i ḟm,n = -( f m+1,n + f m-1,n + f m,n+1 + f m,n-1 -4 f m,n ), in Z 2 ×]0, T ], (7)
f m,n | t=0 = 0, in Z 2 \{(0, 0)}, (8)
f m,n → 0, when |m|, |n| → +∞, (9)
where f 0,0 (t) = δ(t) is the Dirac distribution, and let us compute the related discrete Green's function. By defining the Laplace transform as fm,n (s)

= ∞ 0 f m,n (t)e -st dt,
we obtain, for (m, n) (0, 0), A fm,n = fm+1,n + fm-1,n + fm,n+1 + fm,n-1 ,

setting A = 4si. By using the discrete Fourier transform

F(x, y) = ∞ m=-∞ ∞ n=-∞
e imx e iny fm,n together with [START_REF] Antoine | A friendly review of absorbing boundary conditions and perfectly matched layers for classical and relativistic quantum waves equations[END_REF] gives

F(x, y) = A -4 f0,1 A -2 cos(x) -2 cos(y) . (11)
By inverse discrete Fourier transform, one gets

fm,n (s) = gm,n (s) g0,0 (s) = 1 4π 2 2π 0 2π 0 e -imx-iny A -2 cos(x) -2 cos(y) dxdy 1 4π 2 2π 0 2π 0 1 A -2 cos(x) -2 cos(y) dxdy , ( 12
)
where one has [START_REF] Katsura | Lattice Green's functions for the rectangular and the square lattices at arbitrary points[END_REF] gm,n (s

) = 1 4π 2 2π 0 2π 0 e -imx-iny A -2 cos(x) -2 cos(y) dxdy = 1 (4 -si) m+n+1 (m + n)! m!n! 4 F 3              1 2 (m + n + 1), 1 2 (m + n + 1), 1 2 (m + n) + 1, 1 2 (m + n) + 1; 1 (1 - si 4 ) 2 m + 1, n + 1, m + n + 1              . ( 13
)
By inverse Laplace transform, we obtain

g m,n (t) = (m + n)! m!n! 4i √ π(it) m+n e -4it 2 m+n+2 Γ m + n + 1 2 Γ m + n 2 + 1 2 F 3         m + n + 1 2 , m + n 2 + 1; -4t 2 m + 1, n + 1, m + n + 1         . ( 14
)
The function Γ is the Gamma special function and the generalized hypergeometric function 2 F 3 is

2 F 3 a 1 , a 2 ; z b 1 , b 2 , b 3 = ∞ k=0 (a 1 ) k (a 2 ) k k!(b 1 ) k (b 2 ) k (b 3 ) k z k , (15)
with (a) k = a(a + 1)...(a + k -1) for k ≥ 1 and (a) 0 = 1. We display in Fig. 1 the

(L + 2) × (L + 2) square lattice (0 ≤ m, n ≤ L + 1) arising in (4) for L = 4. 1 ) 1 , 1 ( I 2 ) 2 , 1 ( I 3 ) 3 , 1 ( I 4 ) 4 , 1 ( I 5 ) 4 , 2 ( I 6 ) 4 , 3 ( I 7 ) 4 , 4 ( I 8 ) 3 , 4 ( I 9 ) 2 , 4 ( I 10 ) 1 , 4 ( I 11 ) 1 , 3 ( I 12 ) 1 , 2 ( I 1 ) 1 , 0 ( J 2 ) 2 , 0 ( J 3 ) 3 , 0 ( J 4 ) 4 , 0 ( J 5 ) 5 , 1 ( J 6 ) 5 , 2 ( J 7 ) 5 , 3 ( J 8 ) 5 , 4 ( J 9 ) 4 , 5 ( J 10 ) 3 , 5 ( J 11 ) 2 , 5 ( J 12 ) 1 , 5 ( J 13 ) 0 , 4 ( J 14 ) 0 , 3 ( J 15 ) 0 , 2 ( J 16 ) 0 , 1 ( J Figure 1
: Discrete computational lattice (for L = 4): the red bullets represent the outer layer points, the black bullets correspond to the interior layer points, the yellow bullets are the interior points and the green bullets represent the corner points.

The numerical boundary points are indexed anticlockwise, with the interior layer as

I for = 1, • • • , 4L -4. More explicitly, we set: I = (1, ) for 1 ≤ ≤ L, I +L-1 = ( , L) for 1 ≤ ≤ L, l +2L-2 = (L, L + 1 -) for 1 ≤ ≤ L and I +3L-3 = (L + 1 -, 1) for 1 ≤ ≤ L.
The outer layer points are denoted by J k for k = 1, • • • , 4L, with the same orientation. To this end, let u I (t), for = 1, • • • 4L -4, be given. We first decompose them into (4L -4) independent single-sources κ I k (t), namely,

u I k (t) = 4L-4 =1 f I k -I (t) * κ I (t), 1 ≤ k ≤ 4L -4. ( 16
) Now for (m, n) such that (m, n) ∈ Z 2 \{(p, q)/1 ≤ p ≤ L, 1 ≤ q ≤ L}, we define the Laplace transform of u m,n (t) as ũm,n (s) = 4L-4 =1 g(m,n)-I (s) g0,0 (s) κI (s). (17)
A direct computation shows that, for |m| > 0 or |n| > 0,

isg m,n = -(g m-1,n + gm+1,n + gm,n-1 + gm,n+1 -4g m,n ). ( 18
)
Plugging this expression into (17) yields

(4 -si)ũ m,n (s) -ũm+1,n (s) -ũm-1,n (s) -ũm,n+1 (s) -ũm,n-1 (s) = 4L-4 =1 (4 -is)g (m,n)-I -(g (m-1,n)-I + g(m+1,n)-I + g(m,n-1)-I + g(m,n+1)-I ) g0,0 (s) κI (s) = 0. ( 19
)
It is also easy to check that u m,n (0) = 0, which leads to

iu m,n = -(u m+1,n + u m-1,n + u m,n+1 + u m,n-1 -4u m,n ). (20)
Therefore, (u m,n (t)) m,n is exactly the solution of the semi-discretized system (4) outside the domain {(p, q)/1 ≤ p ≤ L, 1 ≤ q ≤ L}. Once the sources κ I k (t) are obtained, an exact boundary condition reads as follows

u J k (t) = 4L-4 =1 f J k -I (t) * κ I (t), 1 ≤ k ≤ 4L. ( 21
)
Defining λI k (s) = κI k (s) g0,0 (s) , for 1 ≤ k ≤ 4L -4, the boundary conditions can be transformed into

u I k (t) = 4L-4 =1 g I k -I (t) * λ I (t), 1 ≤ k ≤ 4L -4, (22)
and

u J k (t) = 4L-4 =1 g J k -I (t) * λ I (t), 1 ≤ k ≤ 4L. ( 23
)
Combining the boundary conditions ( 22) and ( 23) and the interior equation ( 4), we obtain the reduced semi-discretized problem. The values g m,n (t) are needed in the boundary conditions ( 22) and [START_REF] Antoine | Numerical schemes for the simulation of the two-dimensional Schrödinger equation using nonreflecting boundary conditions[END_REF]. For large m and n, it is not easy to compute g m,n (t) precisely through hypergeometric functions for large t. Consequently, we need to develop an efficient algorithm for computing g m,n (t) accurately. To this end, one easily shows that

gm-1,n -gm+1,n = 1 2π 2 2π 0 2π 0 sin(mx) sin(x) cos(ny) A -2 cos(x) -2 cos(y) dxdy, (24)
leading to

i d ds (g m-1,n -gm+1,n ) = d dA (g m-1,n -gm+1,n ) = - 1 2π 2 2π 0 2π 0 sin(mx) sin(x) cos(ny) (A -2 cos(x) -2 cos(y)) 2 dxdy = 1 4π 2 2π 0 2π 0 sin(mx) cos(ny)d 1 A -2 cos(x) -2 cos(y) dy = -m 1 4π 2 2π 0 2π 0 cos(mx) cos(ny) A -2 cos(x) -2 cos(y) dxdy = -mg m,n . (25)
We obtain

g m+1,n = g m-1,n + img m,n t . ( 26
)
We also have

i(sg 0,0 -i) = -(g -1,0 + g1,0 + g0,-1 + g0,1 -4g 0,0 ), (27)
and

g m,n = g n,m , g -m,n = g m,n , g m,-n = g m,n . (28)
Now collecting ( 18), ( 26), ( 27) and ( 28), we can construct a system of ordinary differential equations to compute g m,n (t) in the domain

{(m, n)|0 ≤ m ≤ L, 0 ≤ n ≤ L}.
The equations in the interior domain can be written as

iġ m,n = -(g m+1,n + g m-1,n + g m,n+1 + g m,n-1 -4g m,n ), (29)
g m,n | t=0 = 0, for {(m, n)|1 ≤ m ≤ L -1, 1 ≤ n ≤ L -1}. ( 30
)
The lower boundary condition is

iġ 0,n = -(g 0,n+1 + g 0,n-1 + 2g 1,n -4g 0,n ), (31)
g 0,n | t=0 = 0, 1 ≤ n ≤ L -1. ( 32
)
Similarly, the left boundary condition reads

iġ m,0 = -(g m+1,0 + g m-1,0 + 2g m,1 -4g m,0 ), ( 33
)
g m,0 | t=0 = 0, 1 ≤ m ≤ L -1. ( 34
)
By ( 26) and ( 28), one has the upper boundary condition

iġ L,n = -( iLg L,n t + g L,n+1 + g L,n-1 + 2g L-1,n -4g L,n ), (35)
g L,n | t=0 = 0, 1 ≤ n ≤ L -1. ( 36
)
In the same way, the right boundary condition writes

iġ m,L = -( iLg m,L t + g m+1,L + g m-1,L + 2g m,L-1 -4g m,L ), ( 37
)
g m,L | t=0 = 0, 1 ≤ m ≤ L -1. ( 38
)
By [START_REF] Baskakov | Implementation of transparent boundaries for numerical solution of the Schrödinger equation[END_REF], we get the lower-left corner condition

iġ 0,0 = -(4g 0,1 -4g 0,0 ), g 0,0 | t=0 = i. ( 39
)
Following a similar derivation, the other corner conditions can be obtained as

iġ L,0 = -( iLg L,0 t + 2g L,1 + 2g L-1,0 -4g L,0 ), g L,0 | t=0 = 0, ( 40
)
iġ 0,L = -( iLg 0,L t + 2g 1,L + 2g 0,L-1 -4g 0,L ), g 0,L | t=0 = 0, ( 41
)
iġ L,L = -( 2iLg L,L t + 2g L,L-1 + 2g L-1,L -4g L,L ), g L,L | t=0 = 0. ( 42
)
In practice, to compute an accurate numerical approximation of g m,n (t), we use the Runge-Kutta 4 (RK4) scheme. If u 0 is known, we then compute u 1 by RK4 on a large enough domain, and use the zero boundary condition. If u 0 is compactly supported and the domain is large enough, the boundary points of u 1 are still 0 and we can use RK4 and closed relations to compute u 1 and so on. Let us also remark that other schemes than RK4 could be used, but of order larger than 2. Nevertheless, stability problems may arise and therefore the ODE solver must be carefully chosen. To illustrate the results, we compare the numerical and exact solutions for g 0,0 (t) in Fig. 2. As the time step ∆t decreases, the approximation quality of g m,n (t) increases. The computed convergence rates of the numerical approximation of g 60,60 [START_REF] Antoine | Absorbing boundary conditions for the one-dimensional Schrödinger equation with an exterior repulsive potential[END_REF] and g 60,60 (36) are reported in Fig. 3. As expected, it is around 4. The algorithm is then able to efficiently and accurately compute g m,n (t) even for large values of the indices m and n.

We also remark that the proposed approach extends to the N-dimensional Schrödinger equation, for which the single-source kernel function for the Schrödinger equation reads

gm 1 ,m 2 ,...m N (s) = 1 (2π) N 2π 0 2π 0 • • • 2π 0 exp        -i N =1 m x        4 -is -2 N =1 cos(x ) dx 1 • • • dx N . ( 43
)
The recurrence relation is then

g m 1 ,m 2 ,...m N (t) = g m 1 ,m 2 ,...m -1,...m N (t) + im g m 1 ,m 2 ,...m N (t) t . (44)
Therefore, one can also build an associated system of ordinary differential equations to calculate numerically g m 1 ,m 2 ,...m N (t). t ℜ(g r e f 0 , 0 (t)) ℜ(g ∆ t 0 , 0 (t)) ℑ(g r e f 0 , 0 (t)) ℑ(g ∆ t 0 , 0 (t))

Figure 2: Numerical and reference solutions for g 0,0 (t) and the semi-discretized Schrödinger equation.

Extension to the heat equation

We extend now the above method to the heat equation

∂ ∂t u = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 u, (x, t) in R 2 ×]0, T ], (45) u
(x, t = 0) = 0, x ∈ Ω ext = R 2 \Ω int , (46)
u(x, •) → 0, for ||x|| → +∞. ( 47
)
We remark that the heat equation can be a transformed Black-Scholes equation. Hence, some applications in finance can also be considered from the present developments. For (L + 1)∆x = a, the semi-discretization of (45) gives

um,n = (u m+1,n + u m-1,n + u m,n+1 + u m,n-1 -4u m,n ), in Z 2 ×]0, T ], (48)
u m,n | t=0 = 0, in (m, n) ∈ Z 2 \{(p, q)|1 ≤ p ≤ L, 1 ≤ q ≤ L}, (49)
u m,n → 0 for |m|, |n| → +∞. ( 50
)
If the computation of (48), ( 49) and ( 50) is confined within the finite domain

{(m, n)|0 ≤ m ≤ L + 1, 0 ≤ n ≤ L + 1},
the boundary conditions read

u I k (t) = 4L-4 =1 g I k -I (t) * λ I (t), 1 ≤ k ≤ 4L -4, u J k (t) = 4L-4 =1 g J k -I (t) * λ I (t), 1 ≤ k ≤ 4L, (51)
with

g m,n (t) = (m + n)! m!n! 4 √ π(t) m+n e -4t 2 m+n+2 Γ m + n + 1 2 Γ m + n 2 + 1 2 F 3         m + n + 1 2 , m + n 2 + 1; 4t 2 m + 1, n + 1, m + n + 1         . (52)
According to [START_REF] Antoine | Absorbing boundary conditions for the two-dimensional Schrödinger equation with an exterior potential part I: construction and a priori estimates[END_REF], one has

g m+1,n = g m-1,n - mg m,n t . ( 53
)
By the same procedure as in Section 2.1, a system of ODEs can be constructed to compute g m,n (t) in the domain

{(m, n)|0 ≤ m ≤ L, 0 ≤ n ≤ L}. The interior equations are ġm,n = g m+1,n + g m-1,n + g m,n+1 + g m,n-1 -4g m,n , in {(m, n)|1 ≤ m ≤ L -1, 1 ≤ n ≤ L -1}, (54)
g m,n | t=0 = 0, in {(m, n)|1 ≤ m ≤ L -1, 1 ≤ n ≤ L -1}. ( 55
)
Finally, the corner and boundary conditions can be obtained by extension. The values of g m,n (t) are computed accurately by RK4. For t ∈ [0, 10], we report in Fig. 4 (left) the approximate solution of g 0,0 (t) and the corresponding exact solutions which are close. In Fig. 4 (right), we display the convergence rate of the approximation of g 60,60 (60). Figure 4: Left: numerical and reference solutions of g 0,0 (t) for the semi-discretized heat equation. Right: Theoretical and numerical slopes of the convergence rates for the approximation of g 60,60 (60).

Numerical discretization

In this section, we design the corresponding discretization of the reduced ordinary differential problem with the accurate artificial boundary condition. For conciseness, we only give the discretized finite difference scheme for the Schrödinger equation. The discretization for the reduced heat equation can easily be transposed. Following subsection 2.1, we adopt the Crank-Nicolson scheme [START_REF] Antoine | Computational methods for the dynamics of the nonlinear Schrödinger/Gross-Pitaevskii equations[END_REF] for the time integration in (4)

u k+1 m,n - i∆t 2 (u k+1 m+1,n + u k+1 m-1,n + u k+1 m,n+1 + u k+1 m,n-1 -4u k+1 m,n ) = u k m,n + i∆t 2 (u k m+1,n + u k m-1,n + u k m,n+1 + u k m,n-1 -4u k m,n ). ( 56
)
After the computation of g m,n (t), the convolutions in [START_REF] Li | Local artificial boundary conditions for Schrödinger and heat equations by using high-order azimuth derivatives on circular artificial boundary[END_REF] are performed through the trapezoidal rule as follows

u k+1 I p = 4L-4 =1 ( ∆t 2 g 0 I p -I λ k+1 I + ∆t k α=1 g k+1-α I p -I λ α I + ∆t 2 g k+1 I p -I λ 0 I ), (57)
for p = 1, ..., 4L -4 and where g α m,n denotes the fourth-order Runge-Kutta approximation of g m,n (α∆t) for the time step ∆t. Rewriting the above equation into the matrix form leads to

                  u k+1 I 1 u k+1 I 2       g 0 I 1 -I 1 2 g 0 I 1 -I 2 2 • • • g 0 I 1 -I 4L-4 2 g 0 I 2 -I 1 2 g 0 I 2 -I 2 2 • • • g 0 I 2 -I 4L-4 2 . . . . . . . . . . . . g 0 I 4L-4 -I 1 2 g 0 I 4L-4 -I 2       ∆t 2 4L-4 =1 g k+1 I 1 -I λ 0 I + ∆t 4L-4 =1 k α=1 g k+1-α I 1 -I λ α I ∆t 2 4L-4 =1 g k+1 I 2 -I λ 0 I + ∆t 4L-4 =1 k α=1 g k+1-α I 2 -I λ α I . . . ∆t 2 4L-4 =1 g k+1 I 4L-4 -I λ 0 I + ∆t 4L-4 =1 k α=1 g k+1-α I 4L-4 -I λ α g 0 J 2 -I 4L-4 2 . . . . . . . . . . . . g 0 J 4L -I 1 2 g 0 J 4L -I 2 2 • • • g 0 J 4L -I 4L-4 2                                          λ k+1 I 1 λ k+1 I 2 . . . λ k+1 I 4L-4                   +                    ∆t 2 4L-4 =1 g k+1 J 1 -I λ 0 I + ∆t 4L-4 =1 k α=1 g k+1-α J 1 -I λ α I ∆t 2 4L-4 =1 g k+1 J 2 -I λ 0 L + ∆t 4L-4 =1 k α=1 g k+1-α J 2 -I λ α I . . . ∆t 2 4L-4 =1 g k+1 J 4L -I λ 0 I + ∆t 4L-4 =1 k α=1 g k+1-α J 4L -I λ α I                    . (59)
In addition, since g m,n (0) = 0 for |m| + |n| > 0, one deduces

                  u k+1 J 1 u k+1 J 2 . . . u k+1 J 4L                   =                    ∆t 2 4L-4 =1 g k+1 J 1 -I λ 0 I + ∆t 4L-4 =1 k α=1 g k+1-α J 1 -I λ α I ∆t 2 4L-4 =1 g k+1 J 2 -I λ 0 L + ∆t 4L-4 =1 k α=1 g k+1-α J 2 -I λ α I . . . ∆t 2 4L-4 =1 g k+1 J 4L -I λ 0 I + ∆t 4L-4 =1 k α=1 g k+1-α J 4L -I λ α I                    . (60)
We denote the vector [

u k I 1 , u k I 2 ..u k I 4L-4 ] T by u k I , [λ k I 1 , λ k I 2 ...λ k I 4L-4 ] T by λ k I and [u k J 1 , u k J 2 ..u k J 4L ]
T by u k J . Therefore, we can recast systems (58) and (60) through a matrix A 1 , and two vectors R k+1 1 , and R k+1 2 as

u k+1 I = A 1 λ k+1 I + R k+1 1 , u k+1 J = R k+1 2 . (61)
Now one can use (61) with (56). All the points (m, n), for 1 ≤ m, n ≤ L, are indexed following

N k = ( k -mod(k, L) L + 1, mod(k, L)), for 1 ≤ k ≤ L 2 . If I k = (m, n), we define I -1 (m, n) = k (and similarly for J k = (m, n) and N k = (m, n)). The vector [u k N 1 , u k N 2 ..u k N L 2 ] T is denoted by u k N .
If we have all the data until the k-th time step,

λ α I for (1 ≤ ≤ 4L-4 and 0 ≤ α ≤ k) and u k m,n for 0 ≤ m, n ≤ L, then we can compute u k+1 J = R k+1 2 .
We can also rewrite (56) as Bu k+1 N = F k+1 , where B is a L 2 × L 2 matrix defined for computing the values u k+1 m,n for 1 ≤ m, n ≤ L and F k+1 is the right-hand side. More precisely, for (m, n), with 2 ≤ m, n ≤ L -1, and from (56) we have

B N -1 (m,n),N -1 (m,n) = 1 + 2i∆t, B N -1 (m,n),N -1 (m-1,n) = - i∆t 2 , B N -1 (m,n),N -1 (m+1,n) = - i∆t 2 , B N -1 (m,n),N -1 (m,n+1) = - i∆t 2 , B N -1 (m,n),N -1 (m,n-1) = - i∆t 2 , F k+1 N -1 (m,n) = u k m,n + i∆t 2 (u k m+1,n + u k m-1,n + u k m,n+1 + u k m,n-1 -4u k m,n ). (62)
For (m, n) = (1, 1), we also can write

(1 + 2i∆t)u k+1 1,1 - i∆t 2 (u k+1 2,1 + u k+1 1,2 ) = u k 1,1 + i∆t 2 u k+1 1,0 + i∆t 2 u k+1 0,1 + i∆t 2 (u k 2,1 + u k 0,1 + u k 1,2 + u k 1,0 -4u k 1,1 ), (63)
and by (61) we have

B N -1 (1,1),N -1 (1,1) = 1 + 2i∆t, B N -1 (1,1),N -1 (1,2) = - i∆t 2 , B N -1 (1,1),N -1 (2,1) = - i∆t 2 , F k+1 N -1 (1,1) = u k 1,1 + i∆t 2 u k+1 1,0 + i∆t 2 u k+1 0,1 + i∆t 2 (u k 2,1 + u k 0,1 + u k 1,2 + u k 1,0 -4u k 1,1 ) = i∆t 2 (R k+1 2 ) J -1 (1,0) + i∆t 2 (R k+1 2 ) J -1 (0,1) + u k 1,1 + i∆t 2 (u k 2,1 + u k 0,1 + u k 1,2 + u k 1,0 -4u k 1,1 ). ( 64
)
We can deal with the points (m, n) = (1, L), (L, 1) and (L, L) similarly to the point [START_REF] Schrödinger | An undulatory theory of the mechanics of atoms and molecules[END_REF][START_REF] Schrödinger | An undulatory theory of the mechanics of atoms and molecules[END_REF]. For (m, n) = (1, ), for 2 ≤ ≤ L -1, one gets

(1 + 2i∆t)u k+1 1, - i∆t 2 (u k+1 1, -1 + u k+1 1, +1 + u k+1 2, ) = u k 1, + i∆t 2 u k+1 0, + i∆t 2 (u k 1, +1 + u k 1, -1 + u k 2, + u k 0, -4u k 1, ), (65)
which means

B N -1 (1, ),N -1 (1, ) = 1 + 2i∆t, B N -1 (1, ),N -1 (1, -1) = - i∆t 2 , B N -1 (1, ),N -1 (1, +1) = - i∆t 2 , B N -1 (1, ),N -1 (2, ) = - i∆t 2 , F k+1 N -1 (1, ) = i∆t 2 (R k+1 2 ) J -1 (0, ) + u k 1, + i∆t 2 (u k 1, +1 + u k 1, -1 + u k 2, + u k 0, -4u k 1, ). (66)
The points (m, n) = ( , 1), (L, ) and ( , L), for 2 ≤ ≤ L -1, can be treated as (1, ). The other coefficients of B are equal to 0. Now we get u k+1 N and all u k+1 m,n can be computed for 1 ≤ m, n ≤ L. To update λ k+1 I , one uses the relation

λ k+1 I = (A 1 ) -1 (u k+1 I -R k+1 1 ), (67)
and the λ k+1 I are next used to build R k+2 2 .

Numerical results

This wave strikes the boundary allowing therefore to analyze the quality of the boundary condition to avoid some spurious reflection that would come back into the computational domain. We fix the spatial grid to L = 46. The numerical reference solution u ref is calculated for a very small time step and on a large domain so that there is no interacting wave with the boundary. We report at time t the absolute 2 -and ∞ -errors between the numerical solution u ∆t (t) := (u ∆t m,n ) 1≤m,n≤L for a time step ∆t and the reference solution u ref (t), and defined by

E 2 ∆t (t) := L m=1 L n=1 |u ref m,n (t) -u ∆t m,n (t)| 2 , E ∞ ∆t (t) := max 1≤m,n≤L |u ref m,n (t) -u ∆t m,n (t)|. ( 69
)
The evolution of the numerical solution with the proposed boundary condition is reported in Fig. 5 for the times t = 0, 12, 14 and 36. We do not observe any reflection at the corners. The errors |u ref m,n (14)u ∆t m,n (14)| (with ∆t = 2 × 10 -2 ) on the grid are displayed in Fig. 6 (left) and with a maximum of the order of 7 × 10 -5 . For completeness, we report u ref 1,46 (t)u ∆t 1,46 (t) (with ∆t = 2 × 10 -2 ) on Fig. 6 (right). In Tables 1 and2, we provide respectively the errors E 2 ∆t (t) and E ∞ ∆t (t), for t = 2, 12, 14, 36, and for various time steps ∆t. We clearly observe that the scheme is at least second-order in time which is coherent with the expected second-order accuracy of the Crank-Nicolson scheme (the numerical rate of convergence are obtained as the least-square approximations of the pointwise rates).

The semi-discretized heat equation

We consider now the case of the semi-discretized heat equation (48) with the initial data uniformly set equal to 10 on a 8 × 8 sub square lattice, and zero otherwise. We choose a 10 × 10 computational domain, which is small enough to test the accuracy of the boundary condition. We report on Fig. 7 3 and4, respectively. Again, a second-order convergence is observed.

Conclusion

In this paper, we investigated the accurate numerical solution of the semi-discretized linear Schrödinger equation by introducing an artificial boundary condition on a rectangular domain to bound the computational domain. The accurate evaluation of the boundary condition is realized by means of some recurrence relations involving the discrete Green's functions. As a result, the corner reflection is fully eliminated. We provide some numerical results that illustrate the second-order accuracy of the proposed scheme with some details about the implementation. Closely related artificial boundary conditions are applied to the heat equation and some numerical simulations show that they are still reflection-free.

Concerning the perspectives, some questions regarding the theoretical and numerical aspects of these boundary conditions remain to be answered. For example, the rigorous proof of the numerical stability and estimates of the convergence rate must still be obtained. Finally, fast algorithms for the boundary conditions also need to be explored in details to accelerate the evaluation of the transparent operator. These questions will be addressed in a forthcoming paper.

Figure 3 :

 3 Figure 3: Theoretical and numerical slopes of the convergence rates for the approximation of g 60,60 (12) (left) and g 60,60 (36) (right).

  or der theor etical s lope l og (|g r ef 6 0 , 6 0 (t = 60) -g ∆ t 6 0 , 6 0 (t = 60)|)

4. 1 .

 1 The semi-discretized Schrödinger equation Let us consider the semi-discretized Schrödinger equation (4) with the discrete initial data u m,n (0) = exp -m -

Figure 5 :

 5 Figure 5: Amplitude of the numerical solution |u ∆t | (with ∆t = 2 × 10 -2 ) of the Schrödinger equation on the 46 × 46 square computational domain at different times t (upper-left for t = 0; upper-right for t = 12; lower-left for t = 14; lower-right for t = 36).

ℜ(u r e f 1 ,

 1 4 6 (t) -u ∆ t 1 , 4 6 (t)) ℑ(u r e f 1 , 4 6 (t) -u ∆ t 1 , 4 6 (t))

Figure 6 :

 6 Figure 6: Left: absolute error |u ref m,n (14)u ∆t m,n (14)| (with ∆t = 2 × 10 -2 ) for the Schrödinger equation on the 46 × 46 square computational domain. Right: error on the real and imaginary parts at the upper left corner.

  and for ∆t = 2 × 10 -2 the numerical solution in the truncated domain. We can see that there is no corner reflection. The errors |u ref m,n (40)u ∆t m,n (40)| on the grid E 2 ∆t (t) t = 2 t = 12 t = 14 t = 36 ∆t = 8 × 10 -2 2.98 × 10 -2 1.31 × 10 -1 8.37 × 10 -2 9.00 × 10 -3 ∆t = 5 × 10 -2 1.16 × 10 -2 5.11 × 10 -2 3.25 × 10 -2 3.58 × 10 -3 ∆t = 4 × 10 -2 7.44 × 10 -3 3.26 × 10 -2 2.07 × 10 -2 2.30 × 10 -3 ∆t = 2 × 10 -2 1.84 × 10 -3 8.05 × 10 -3 5.11 × 10 -3 5.79 × 10 -4 ∆t = 1 × 10 -2 4.38 × 10 -4 1.92 × 10 -3 1.22 × 10 -3 1.45 × 10 -4 ∆t = 5 × 10 -3 8.76 × 10 -5

  ∆t = 8 × 10 -2 2.61 × 10 -2 3.15 × 10 -2 3.38 × 10 -2 1.61 × 10 -2 ∆t = 5 × 10 -2 1.01 × 10 -2 1.23 × 10 -2 1.31 × 10 -2 6.25 × 10 -3 ∆t = 4 × 10 -2 6.49 × 10 -3 7.84 × 10 -3 8.40 × 10 -3 3.99 × 10 -3 ∆t = 2 × 10 -2 1.62 × 10 -3 1.96 × 10 -3 2.10 × 10 -3 7.97 × 10 -4 ∆t = 1 × 10 -2 4.05 × 10 -4 4.89 × 10 -4 5.24 × 10 -4 2.49 × 10 -4 ∆t = 5 × 10 -3 1.01 × 10 -4 1.22 × 10 -4

  ∆t = 8 × 10 -2 3.99 × 10 -3 4.53 × 10 -3 3.68 × 10 -3 1.75 × 10 -3 ∆t = 5 × 10 -2 1.55 × 10 -3 1.76 × 10 -3 1.43 × 10 -3 6.81 × 10 -4 ∆t = 4 × 10 -2 9.93 × 10 -4 1.13 × 10 -3 9.15 × 10 -4 4.35 × 10 -4 ∆t = 2 × 10 -2 2.48 × 10 -4 2.81 × 10 -4 2.28 × 10 -4 1.09 × 10 -4 ∆t = 1 × 10 -2 6.24 × 10 -5 7.03 × 10 -5 5.71 × 10 -5 2.72 × 10 -5 ∆t = 5 × 10 -3 1.59 × 10 -5 1.76 × 10 -5 1.43 × 10 -5

Figure 7 :

 7 Figure 7: Amplitude of the numerical solution |u ∆t | (with ∆t = 2 × 10 -2 ) of the heat equation on the 10 × 10 square computational domain at different times t (upper-left for t = 0; upper-right for t = 1; lower-left for t = 10; lower-right for t = 40).

Figure 8 :

 8 Figure 8: Absolute error |u ref m,n (4)u ∆t m,n (4)| (with ∆t = 2 × 10 -2 ) for the heat equation on the 10 × 10 square computational domain.

Table 1 :

 1 [START_REF] Hazeltine | The Framework of Plasma Physics[END_REF].84 × 10 -4 2.44 × 10 -4 3.63 × 10 -5 Error E 2 ∆t (t) at different times for various time steps ∆t (Schrödinger equation). × 10 -2 3.86 × 10 -3 1.94 × 10 -2 2.08 × 10 -2 3.47 × 10 -4 ∆t = 5 × 10 -2 1.51 × 10 -3 7.56 × 10 -3 8.10 × 10 -3 1.38 × 10 -4 ∆t = 4 × 10 -2 9.66 × 10 -4 4.83 × 10 -3 5.17 × 10 -3 8.88 × 10 -5 ∆t = 2 × 10 -2 2.39 × 10 -4 1.19 × 10 -3 1.28 × 10 -3 2.24 × 10 -5 ∆t = 1 × 10 -2 5.69 × 10 -5 2.84 × 10 -4 3.05 × 10 -4 5.61 × 10 -6 ∆t = 5 × 10 -3 1.14 × 10 -5 5.71 × 10 -5 6.22 × 10 -5 1.40 × 10 -6

	rate	2.09	2.09	2.09	2.00
	E ∞ ∆t (t) ∆t = 8 rate	t = 2 2.09	t = 12 2.09	t = 14 2.08	t = 36 1.99

Table 2 :

 2 E ∞ ∆t (t) at different times for various time steps ∆t (Schrödinger equation).

are displayed in Fig 8

.

The reference solution is calculated over a much larger domain and again with a very small time step. The errors in 2 -and ∞ -norms with respect to ∆t are reported in Tables

Table 3 :

 3 1.31 × 10 -4 6.23 × 10 -5 E 2 ∆t (t) at different times for various time steps ∆t (heat equation).

	rate	2.00	2.00	2.00	2.00

Table 4 :

 4 E ∞ ∆t (t) at different times for various time steps ∆t (heat equation).

	6.79 × 10 -6

Acknowledgements

This study is partially supported by the National Natural Science Foundation of China under contract numbers 11272009, 11521202 and 11502028. X. Antoine thanks the support of the LIASFMA (funding from the University of Lorraine) and the French ANR grant Bond (ANR-13-BS01-0009-01).