E D I T O R I A L

Increasing our knowledge on grapevines physiology to increase yield, quality and sustainably Scandinavia is not (yet) a famous wine-producing area, but the Scandinavian Plant Physiology Society offered me the opportunity to highlight some of the recent scientific contributions on grapevines published in Physiologia Plantarum. I first would like to apologize to many colleagues who published important contributions in the field and I did not cite. The aim of this short communication was to take some examples in the literature to present to scientific wine lovers, but also the ones who like grape juice and table grapes or raisins, which topics grapevine scientists are studying to increase our knowledge in grapevine physiology that could at the end help winegrowers to increase yield, quality and sustainability.

Even if the storytelling around domestication is so interesting [START_REF] Badouin | The wild grape genome sequence provides insights into the transition from dioecy to hermaphroditism during grape domestication[END_REF][START_REF] Terral | Evolution and history of grapevine (Vitis vinifera) under domestication: new morphometric perspectives to understand seed domestication syndrome and reveal origins of ancient European cultivars[END_REF], for the purpose of this article, we should remember that Vitis vinifera is the most cultivated one.

Several thousand varieties have been listed [START_REF] Galet | Dictionnaire encyclopédique des cépages et de leurs synonymes[END_REF].

Grape is an important crop with high economic value for several countries, from the old European producing ones to the modern ones (North and South America, South Africa, Australia, New Zealand) and the coming one (China). Some specificities should be remembered to understand the complexity of grape production and research; within these, we can cite (1) a wide range of yields (from 15 to 200 hl/Ha), (2) a not-fully understood, and even defined, quality linked to the notion of terroir,

(3) growth at the same place since 200 years precluding crop rotations, (4) in most cultivated vineyards, the presence of grafted plants that were not introduced for yield or quality but for disease resistance (phyloxera), (5) a strong diversity in agropedoclimatic growth conditions, (6) an important development of organic and esoteric biodynamic culture, (7) an important spreading of wood diseases such as eutypiosis and esca, and (8) slow and low plant breeding outputs because of a reactionary attitude of both producers and consumers.

Like most crops, vineyards are facing biotic and biotic stresses.

Although drought has been studied for a long time, climate change has induced an increase in the studies dealing with temperature, light quality and intensity and CO 2 . Unfortunately, Vitis vinifera is very sensitive to several pathogens, the two most important being downy and powdery mildew, but other ones are also problematic in more specific areas such as trunk disease, black rot, Xylella fastidiosa. The combinatorial effects of these stresses cannot be predicted from the study of each stress studied alone [START_REF] Bouain | Getting to the root of plant mineral nutrition: combinatorial nutrient stresses reveal emergent properties[END_REF]. A recent and counterintuitive demonstration is the absence of drought effects on esca development [START_REF] Bortolami | Grapevines under drought do not express esca leaf symptoms[END_REF].

Although it was claimed above that new grapevine varieties are rare and infrequent, two exceptions should not be forgotten: disease (mildew)-resistant varieties and the use of genetic as a research tool.

Besides, some historical varieties from interspecific crosses are resistant to downy and powdery mildew (Jacquez, Clinton, Maréchal Foch), new breeding programs are developed to create new resistant varieties with properties closer to the international varieties. Several new varieties are available for the producer, even in conservative countries. Other breeding programs take into account climate changes adaptation [START_REF] Bigard | Vitis vinifera L. diversity for cations and acidity is suitable for breeding fruits coping with climate warming[END_REF].

| WATER USE EFFICIENCY AND DROUGHT

Grapevine is often cultivated in drought areas. The efficacity to use water to produce a defined biomass can be defined at different levels and take different names: Water Use Efficiency (WUE) or Transpiration Efficiency. The importance of such parameters lies in the link between water loss and CO 2 fixation through the stomatal pores. The regulation of stomatal aperture has been studied in different conditions and varieties. Daytime regulation is well-documented, whereas the role of nighttime transpiration has been demonstrated only recently [START_REF] Coupel-Ledru | Reduced nighttime transpiration is a relevant breeding target for high water-use efficiency in grapevine[END_REF].

Several studies were performed to compare the variability in WUE between different cultivars, for example, by comparing (near) isohydric (Grenache) and (near)anisohydric cultivars (Syrah) [START_REF] Coupel-Ledru | Reduced nighttime transpiration is a relevant breeding target for high water-use efficiency in grapevine[END_REF]. [START_REF] Buesa | Intracultivar genetic diversity in grapevine: water use efficiency variability within cv. Grenache[END_REF] studied the intracultivar genetic diversity using 13 Grenache genotypes (both red and white) under field conditions over three seasons. They demonstrate that three parameters (intrinsic WUE, crop WUE, and δ 13 C) are reliable indicators to select drought-resistant genotypes.

While it is expected that the growing conditions are important to study the effect of water availability of the different (eco)physiological parameters, it is also important to quantify the drawbacks that can occur. The effect of the container's volume on the plant response to drought has been studied by comparing several parameters on plants grown in 7 or 20 L pots [START_REF] Herrera | Container volume affects drought experiments in grapevines: insights on xylem anatomy and time of dehydration[END_REF]. Before drought

| PRIMARY AND SECONDARY METABOLISM

One important berry component is the sugar content. Sugar is made from CO 2 in the source (the mature leaves) and transported to the sink: berries, young leaves, roots/trunks. Combinatorial stress application is necessary because plant response cannot be inferred from their response to a single stress. The three parameters that will change in the coming years, namely, CO 2 , temperature and water availability, have been studied in two Tempranillo cultivars [START_REF] Kizildeniz | Future CO2, warming and water deficit impact white and red Tempranillo grapevine: photosynthetic acclimation to elevated CO2 and biomass allocation[END_REF]: CO 2 and temperature effect are cultivar-and combinatorial dependent. Photosynthetic acclimation is observed under long-term exposition to elevated CO 2 .

The relationships between photosynthesis, carbohydrate metabolism and expression of sugar-related genes in differentially light-exposed leaves (East/West/Sun/Shade) of Malbec cultivar were studied [START_REF] Dayer | Non-structural carbohydrates and sugar export in grapevine leaves exposed to different light regimes[END_REF]. East, West and Sun leaves exported sugars. Shade leaves accumulated starch and sugars, but the origin of these carbohydrates is unknown since the photosynthetic activity of these leaves is negligible.

| PATHOGENS

Grapevine trunk diseases (Esca, Black Dead Arm, Eutypa, …) are developing in Europe since the ban of sodium arsenite in 2003. Yield loss can be as much as 20%. Leaf metabolite fingerprints have been performed on two different clones of Trousseau, an esca-susceptible cultivar [START_REF] Moret | Esca-affected grapevine leaf metabolome is clone-and vintage-dependent[END_REF]. A clone-and year-dependent metabolic response is demonstrated. Furthermore, comparison with data obtained on Chardonnay [START_REF] Moret | Clone-dependent expression of esca disease revealed by leaf metabolite analysis[END_REF] did not allow identifying esca-specific metabolic fingerprints.

The nonproteinaceous amino-acid β-aminobutyric acid (BABA) is a priming defense elicitor known to induce broad-spectrum disease resistance. It activates the expression of several defense-related genes as PR-proteins and phytoalexin biosynthesis genes. A transcription factor (VvWRKY18) is identified as a BABA-priming regulator [START_REF] Wang | Dual function of VvWRKY18 transcription factor in the β-aminobutyric acid-activated priming defense in grapes[END_REF]

| MINERAL NUTRITION

Climate change is also affecting grapevine mineral nutrition. If we take into account only two major nutrients, potassium and nitrate, a strong effect on berry quality can be expected. For potassium, its increase in berries is correlated with a decrease in acidity, often a source of lowquality wines. Potassium flux in the plant is affected by climate changes and leads to an increase in berry potassium concentration [START_REF] Monder | Adjustment of K+ fluxes and grapevine defense in the face of climate change[END_REF]. Nitrogen is involved in the growth, development and synthesis of aromatic secondary metabolites. Furthermore, berry nitrogen (YAN, Yeast assimilable nitrogen) content is important for the fermentation process. Nitrogen availability, uptake, and metabolism affect berry nitrogen composition. The effects of agropedoclimatic condition on N nutrition have been recently reviewed [START_REF] Verdenal | Understanding and managing nitrogen nutrition in grapevine: a review[END_REF]. Mineral nutrition should be studied as an inclusive process taking into account interaction between nutrients [START_REF] Bouain | Getting to the root of plant mineral nutrition: combinatorial nutrient stresses reveal emergent properties[END_REF] as well as all the other abiotic and biotic factors.

| TECHNICAL ADVANCES

A technic that was first developed for research purposes is isotopic discrimination. Several physiological processes discriminate between the two stable carbon isotopes, 12 C and 13 C [START_REF] Farquhar | Carbon isotope discrimination and photosynthesis[END_REF].

The ratio between the two isotopes can be represented as the δ 13 C. The use of this δ 13 C on grapevine samples was reported in this journal in 1977 [START_REF] Marco | Carbon isotope ratios (13C/12C) in fractions of field-grown grape[END_REF] and was later demonstrated to be a good proxy of the summer grapevine water status [START_REF] Gaudillère | Carbon isotope composition of sugars in grapevine, an integrated indicator of vineyard water status[END_REF]. This parameter is now quantified in routine by the oenology laboratory for the winemakers.

OMICS approaches have been used in a lot of studies to understand plant responses to the environment and plant development. The ones who have eaten berries from the same cluster have observed that they do not all have the same taste. In spite of that, "representative" sampling is still performed. Recently, a thousand berries were collected and berry weight, sugar import and malate breakdown were quantified [START_REF] Shahood | First quantitative assessment of growth, sugar accumulation and malate breakdown in a single ripening berry[END_REF]. Kinetic analysis allowed us to build an improved model of berry growth and ripening. Should we go smaller in the analysis using the trendy single-cell approach?

| CONCLUSION

A number of challenges are faced by producers and scientists to cope with climate change. Being already cultivated in semi-arid areas, an increase in drought will not be fully counteracted by practice management such as irrigation or architectural forms [START_REF] Deloire | What could be the architectural forms of future vines adapted to climate change: a new challenge! Let's discuss the Gobelet (bush vine)[END_REF].

Exploiting natural diversity [START_REF] Wolkovich | From pinot to Xinomavro in the world's future wine-growing regions[END_REF] and developing breeding strategies [START_REF] Bigard | Vitis vinifera L. diversity for cations and acidity is suitable for breeding fruits coping with climate warming[END_REF] are required. Our understanding of abiotic stress effects due to climate change on different grape varieties was recently reviewed [START_REF] Carvalho | Cutting the Gordian knot of abiotic stress in grapevine: from the test tube to climate change adaptation[END_REF]. The cultivar-dependent resilience to heatwaves (up to 47 C in France in June 2019) should be studied because it cannot be predicted from response to other stress. For example, Carignan, a cultivar known to be drought-resistant or drought-adapted, is the one that was the more heat-sensitive. Whatever the stresses, recurrence is increasingly Email: benoit.lacombe@cnrs.fr

(

  46 days), the number of leaves and average leaf area were the same between plants grown in containers of different volume, but differences in xylem vessel anatomy have been identified. Furthermore, the plant response to drought is different. But what should be clarified is the origin of the difference. Is the stress sensed by the plant different or is the plant anatomy at the beginning of the stress different? Indeed, differences in nutrient availability could induce anatomy differences. Besides the importance of the development stage where the drought stress occurs, the effect of the speed of the stress should also be taken into account. Morabito et al. (2022) compared fastdeveloping drought to slow-developing drought on stress response and recovery. They demonstrate that Grenache, usually described as a near-isohydric cultivar, can behave as a near-anisohydric cultivar when the stress is slowly applied.

  Carbon balance calculation needs to take into account carbon inputs and outputs. These parameters have been quantified in leaves in a plethora of conditions, but fruit respiration has been far less considered.[START_REF] Hernández-Montes | Plant water status and genotype affect fruit respiration in grapevines[END_REF] studied this in fields conditions over two seasons at different stages in irrigated and nonirrigated Tempranillo and Grenache vineyards. Cluster respiration decrease during development. Tempranillo loss more carbon than Grenache. What is so (the most?) important in the berry intended to wine production? The amount and characteristics of the secondary metabolites. The biosynthetic gene expression regulation, characterization and quantification of the different metabolites are under study. Del-Castillo-Alonso et al. (2021) studied the effect of UV radiation on different physiological parameters and could demonstrate that ambient UVs increase flavonol contents and expression of flavonol synthase and chalcone synthase, with UV-B having a stronger effect than UV-A.

  . VvWRYK18 interacts with VvNPR1 and activates the transcription of Stilbene Synthase genes. Using specific organisms to negatively control phytopathogens is a biocontrol strategy developed in vineyards. Besides different seaweed, Trichoderma spp., soil fungi, are used and their mechanisms of action are partially understood. Lazazzara et al. (2021) studied the Volatile Organic Compounds (VOCs) profile of three Trichoderma strains and identified 31 compounds. They tested on leaf discs five of the identified VOCs against downy mildew, callose deposition and hypersensitive response and two of them on gene expression. All VOCs tested act as biofungicides against downy mildew.

Figueiredo

  [START_REF] Figueiredo | An apoplastic fluid extraction method for the characterization of grapevine leaves proteome and metabolome from a single sample[END_REF] developed a protocol to prepare samples of leaves apoplast, which were used to perform proteomic and metabolomic analyses. nanoLC-MS/MS and FT-ICR-MS identified 700 proteins and 514 metabolites.

  occurring and effects on plant development and yield have been observed. The cumulative effects of stress over the years are a challenge that should be tackled. Benoit Lacombe IPSiM, Université de Montpellier, CNRS, INRAE, Institut Agro, Montpellier Cedex, France Correspondence Benoit Lacombe, IPSiM, Université de Montpellier, CNRS, INRAE, Institut Agro, Place Viala, 34060 Montpellier Cedex, France.

Physiologia Plantarum


Pedro Lopez-Merino

Juliette Rouchier

An agent-based model of (food) consumption: Accounting for the Intention-Behaviour-Gap on three dimensions of characteristics with limited knowledge

Keywords: Intention-Behaviour-Gap, Agent-based modelling, Innovation diffusion, Lancaster goods

We propose an agent-based model to study purchase diffusion. We conceive goods with 3 dimensions of characteristics, price premiums and limited knowledge on the part of consumers. Despite purchase intention for each characteristic being high, only a reduced number of consumers are initially able to acquire them. The model thus reproduces the Intention-Behaviour-Gap often identified in sustainable food consumption, by exploring two of its known sources: price premiums on the presence of extra characteristics, and lacks in consumers' knowledge as to which goods contain them. We analyse the extent to which knowledge of characteristics diffuses throughout the population and purchases of them are adopted. By testing how different parameters affect these evolutions, we offer insights as to how wider adoption of desirable purchase behaviours can be encouraged. Results show that all parameters have significant effects on knowledge availability and purchase behaviour, and that the ensuing increased knowledge particularly affects purchases of combined characteristics-an interesting and unexpected result. Modifying network parameters (average network degree and knowledge spillovers) produces effects comparable to those of external ones 1 (initial knowledge availability and price premiums), an interesting feature in terms of policy recommendations since the former can arguably imply less costly interventions than the latter.

Introduction 1.Motivation for the model

In many respects, the health, environmental and social impacts of the contemporary globalised food system are deemed negative [START_REF] Tscharntke | Global food security, biodiversity conservation and the future of agricultural intensification[END_REF][START_REF] Lamine | Sustainability and Resilience in Agrifood Systems: Reconnecting Agriculture, Food and the Environment: Sustainability and resilience in agrifood systems[END_REF]. Recommendations at all levels [United Nations General Assembly, 2015, IPES-Food, 2017] promote a higher consideration of health, environment and society across the food industry. While the production and distribution sides of the issue have been the subject of a relatively large body of research and policy recommendations, consumption remains somewhat understudied [START_REF] Carroll | Locating the locale of local food: The importance of context, space and social relations[END_REF]Fahy, 2015, Keller et al., 2016].

Consumers' general awareness and interest in these issues have increased over the past few years, although actual purchases of goods that can qualify as healthy, environmentally sustainable and socially optimal remain lower than declared intentions. Studies find price premiums and lack of knowledge on the part of consumers as to which goods contain the said characteristics to be amongst the reasons for this intention-behaviour gap (IBG) [START_REF] Joshi | Alimentation durable : les Français de plus en plus attentifs à ce qu'ils mangent[END_REF].

Theoretical approach

We conceive goods as (binarily) satisfying or not the characteristics of being healthy, good for the environment and good for society. In this, we follow the seminal work of [START_REF] Lancaster | A New Approach to Consumer Theory[END_REF], who introduced the notion of goods as collections of characteristics that consumers seek to satisfy.

Evidence shows that each of our three dimensions of analysis (i.e. health, environment, society) have in the past few years been increasingly adopted by [French] consumers [START_REF] Joshi | Alimentation durable : les Français de plus en plus attentifs à ce qu'ils mangent[END_REF][START_REF] Joshi | Alimentation durable : les Français de plus en plus attentifs à ce qu'ils mangent[END_REF]. The innovation-diffusion approach thus appears as an appropriate framework to describe current trends and study possible evolutions, since it was conceived to study how groups of consumers move from low to high adoption. For this, we draw from an agent-based model (ABM) presented by [START_REF] Leite | Innovation diffusion with heterogeneous networked agents: a computational model[END_REF].

Our work, however, expands theirs in two different ways: first, by using Lancaster-type goods we allow for multiple dimensions of innovations to run in parallel and-to some extent-in competition (as consumers may have to make trade-offs between them). Second, we study intention-behaviour gaps generated by price differences and lack of knowledge on the part of consumers, thus ac-counting for a typical finding of sustainability studies (and of the aforementioned surveys): that of discrepancies between intention and action.

Two dimensions of diffusion are present in our analysis: knowledge and purchase. The former is a direct output of social interaction, whereby individual consumers can gain knowledge from one another as to what goods contain which characteristic. The second is indirect, since the acquisition of knowledge will permit consumers to assess whether they are able or not to buy the good they want, and so purchase is not an immediate result of interaction.

Working hypotheses, objectives and main findings

In order to proceed with our modelling work, we establish two main working hypotheses, thus far impossible to verify in the literature: H1: The level of interest on a given characteristic determines whether or not a consumer will initially have the knowledge of which goods contain it. H2: Individuals with knowledge of a characteristic and the ability to purchase it can influence others' knowledge and thus drive diffusion.

H1 is derived from the intuition that the higher the interest in a given dimension (e.g. environmental sustainability), the more a consumer will seek information as to how to satisfy it, thus increasing his or her chances of gaining knowledge. H2 is derived from a common element in innovation diffusion studies, whereby experience of earlier adopters trickles down and can influence others' adoption, as well as from the fact that individuals largely seek information on their social networks [START_REF] Everett | Diffusion of preventive innovations[END_REF]. These hypotheses are worth discussing further, and we'll get back them at the end of the article.

Other than these, our model is built under a specific economic assumption: goods objectively have the quality of containing or not a given characteristic, and there is no subjectivity in identifying them, nor an evolution over time of what constitutes healthy, sustainable or socially fair food.

With regards to our subject-matter (food), there is reason to think that the process of knowledge diffusion for each of the characteristics is in the early stages of diffusion, with some 80% of respondents declaring they lack information as to how to identify them on products they purchase [START_REF] Joshi | Alimentation durable : les Français de plus en plus attentifs à ce qu'ils mangent[END_REF][START_REF] Joshi | Alimentation durable : les Français de plus en plus attentifs à ce qu'ils mangent[END_REF]. We thus calibrate the model so that about 20% of consumers have knowledge for each of the characteristics from the beginning. Conversely, preferences are preset so that 60% of consumers have the intention of buying each characteristic.

The model recreates a situation in which consumers with experience buying each dimension disseminate knowledge of them and in turn can influence adoption over time. Although goods containing extra dimensions are pricier than their more basic counterparts 1 , sufficiently interested consumers will pay the price premium on condition of knowing where they can be found. Those who do not adopt a given dimension are either uninterested in it, do not have the knowledge of where to find it, or are unable to pay for it-these last two being behind the aforementioned IBG.

When studying this model, we seek to observe the relationship between purchases, knowledge and price premiums for networked consumers. We look at the extent to which knowledge is a source of IBG and how its diffusion throughout the simulation can reduce the said gap. More particularly, we seek to find emerging properties in the relationship between knowledge diffusion and IBG for combined characteristics.

Model description

Consumers

We initialise 400 agents ("consumers") each of whom is randomly linked to s i others, s i ∼ U (1, 2S). The average number of links per consumer is thus S.

Each consumer is endowed with an individual budget w i ∼ U ($1, W ).

They have a randomly-assigned individual preference for each of the dimensions p i,a ∼ U (0, 1), a = {1, 2, 3}. With regards to these preferences:

-A common intention-threshold H (0 < H < 1) determines the point above which preference becomes intention:

p i,a > H → h i,a = 1
-A common knowledge-threshold K (H < K < 1) determines the point above which preference is associated with knowledge of a dimension:

p i,a > K → k i,a = 1

Goods

We create a space of 8 types of goods, depending on whether each of the dimensions is present or not (2 3 ). Types of goods thus range from [0 0 0] (a basic good) to [1 1 1] (a complete one).

The presence of a dimension in a good comes with a price premium. We set x, x+10% and x+15% as the premium for the presence of one, two and three characteristics on a given good. The price of a basic good being set at $1, that of a complete one will then be $(1+x + 10% + 5%). The supply for each type of good is fully elastic (there is no limit to the quantity that can be supplied at the price of a good, and consumers are price-takers).

extra 10% on the second, and a further 5% on the third, resulting in a 48% premium on a good containing all three characteristics.

Time-dynamics

At each time-step consumers have to purchase one unit of a type of good (which can be pictured as a weekly basket). They will aim for their most preferred good according to p i,a and k i,a .

A consumer unable to purchase their preferred good due to an insufficient budget will settle for a sub-optimal one. For this, they will drop the dimension(s) for which they have the lowest preference and aim to buy the next relevant good.

Consumers who are able to purchase a given dimension accumulate experience following the rule of [START_REF] Leite | Innovation diffusion with heterogeneous networked agents: a computational model[END_REF] :

E t i,a = 1 -(t -t i,a ) -1 , if t > t i,a 0, otherwise (1)
Where E t i,a is the amount of experience on characteristic a accumulated by consumer i up to the t-th time-step, and t i,a is the time-step at which he or she first knowingly bought it.

Consumers interested in a given dimension but without knowledge of it (H < p i,a < K) are unable to knowingly buy it2 . They thus seek it from other consumers in their immediate network who have been able to purchase it, following:

P r(k t i,a = 1) = l jn j1 E t j,a , if l jn j1 E t j,a < 1 1, otherwise (2)
Where j 1 , ..., j n represent each of the nodes in i's immediate network, and l is an experience spillover parameter. Equations 1 and 2 imply that the more experience i 's network gains on a given characteristic, the likelier it will be for i to gain knowledge of which goods contain it.

The reader will notice that preferences are fixed throughout the simulation. This means that all through our runs, only consumers whose preference on a given dimension was initially sufficiently high (p i,a > H) will have the intention to buy, and that whether they do or not will depend on budget and knowledge considerations only. In this, the model can be seen as a snapshot in which neither preferences nor prices evolve, and diffusion takes place only through knowledge spreading3 .

Initial setup

We initially parametrise our model in the following way4 .

Fixed parameters

-W = $2 -H = 0.4
(Consumers are endowed with between 1 and 2 dollars for their food purchases per unit of time) (For any given dimension, 60% of consumers have the intention of purchasing it, whereas 40% don't)

Variable parameters

-l = 0.42 -S = 6 -K = 0.81 -x = +33%
(The knowledge spillover per time-step is of 42% of the experience accumulated by agents on a consumers' network (Consumers are connected on average to 6 others) (81% of consumers do not have spontaneous knowledge of the goods containing each dimension. Only 19%-those with a very high preference of a dimension-do) (Price premiums are of 33, 43 and 48% for the presence of 1, 2 and 3 characteristics in a good)

We run our model under this baseline setup in order to study the evolution of a number of indicators described in section 2.3.

Parameter modification

After testing our model under the baseline setup described above, we look at how variations in our chosen parameters affect its evolution. We modify each parameter individually by ±1/3 rd (±33.3%), starting from the baseline situation. The direction chosen for each parameter is that which increases knowledge circulation and adoption: +33.3% for l and S, -33.3% for K and x.

The parameters we test allow us to assess how our model reacts to higher or lower knowledge circulation (l), larger or smaller networks (S), modifications in the price premiums (x), and a higher or lower knowledge threshold (K).

In doing so, we study how network-related parameters (l and S) affect the model's dynamics, as opposed to external ones (x and K). With this, we seek to gain insight into how an increased social interaction impacts purchases and the intention-behaviour gap, and compare it to the related effects of having external interventions such as subsidies, better labelling of characteristics or rich literature uses ABM to study the evolution of opinion (preference), and has the potential of producing interesting results when combined with our approach. See [START_REF] Deffuant | An Individual-Based Model of Innovation Diffusion Mixing Social Value and Individual Benefit[END_REF], [START_REF] Huet | Resisting hostility generated by terror: An agent-based study[END_REF], [START_REF] Butler | An opinion diffusion model with deliberation[END_REF], [START_REF] Barbet | Tension Between Stability and Representativeness in a Democratic Setting[END_REF]. communication campaigns.

Observed Indicators

When studying the evolution of the model, we pay attention to a number of indicators, described in Table 2 below.

Indicator Description %know1D the average percentage of consumers with knowledge of one, two or three of the dimensions %know2D %know3D %buy1D

the average percentage of consumers actually purchasing one, two or three of the dimensions %buy2D %buy3D %IKG the average intention-knowledge gap per dimension %IBG the average intention-behaviour gap per dimension %KBG the average knowledge-behaviour gap per dimension IBG refers to the consumers (as a proportion of the total population) who have the intention to buy a dimension, yet fail to do so because of knowledge or budget limitations. IKG refers to how many have intention but do not have knowledge of a dimension, and KBG to how many have knowledge (and thus intention) but fail to buy because of budget limitations only. Other than IBG, these concepts are not drawn from the existing literature and can thus be seen as an addition. IKG, IBG and KBG are normalised to 100% to account for the fact that there are three dimensions.

There is a trivial relationship between the last three indicators, since the IBG is caused either by knowledge or by budget limitations (or both). Note however that KBG + IKG does not necessarily equate IBG because of the unknowing pickups we spoke about earlier.

For each configuration of parameters, we run the model 25 times during 80 time-steps (the time it normally takes for evolving indicators to stabilise). We look at indicators at time-steps 1 and 80.

3 Simulation results We first run the model on the basic setup we propose. Figure 1 above and Table 3 below synthesise the evolution of the main indicators of our virtual world. Left to itself, each indicator evolves towards a certain equilibrium, the main observations being:

Baseline setup

-The evolution of indicators for knowledge, IBG, IKG and KBG, as well as (albeit less visibly) purchases shows an evolution producing S-shaped curves-a usual finding in innovation diffusion studies. This is not entirely unexpected since we built our model on the basis of a preexisting one that already showed these features. Still, since the majority of indicators (IBG, IKG, KBG) are new constructions, this is a relevant finding.

-Both knowledge and purchases of one (blue), two (green) and three (red) dimensions increase. Since neither prices nor preferences vary, this means that word-of-mouth-type information is circulated and helps consumers identify the dimensions they look for. Increases in actual purchases are purely the result of knowledge gains.

-Knowledge of at least one dimension moves from being available to less than 46.51% of the population initially, to 66.55% by t = 80, an increase of 20.04%. This means that two-thirds of our consumers are able to identify at least one characteristic in a product containing it by the time the simulations stabilise. Knowledge of two and three dimensions also increase, the proportion of consumers able to identify all three dimensions reaching 4,39% by the end of baseline simulations.

-Purchases also evolve, as can be gathered from inspecting Figure 1 and Table 3. Since neither prices nor preferences evolve, all increases in purchases can be attributed to knowledge gains, and they are more pronounced for consumers buying 2 or 3 dimensions than for those buying one alone.

-The IKG drops from 40.60% to 27.44% over the course of the simulations.

-The IBG also drops, though its fall is more moderate, moving from 35.28% to 31.07%.

-The KBG nearly doubles, to 11.52%. This, because gains in knowledge place an increasing number of consumers in a situation of incapacity due to price premiums.

These initial results open up a number of issues. As stated earlier, we know that consumers currently declare lack of knowledge as an important limiting factor in their food purchases. In the reality of our model, we expect knowledge to disseminate by itself to an important proportion of the population (a result already observed, under certain circumstances in earlier works by [START_REF] Rouchier | Opinion transmission in organizations: an agent-based modeling approach[END_REF]). We'll see however, that certain parameter modifications produce stronger knowledge gains.

The observed increase in knowledge translates as changes in purchases in what is perhaps the most interesting finding of the model: the proportion of consumers buying one of the dimensions shows little change, whereas that of consumers buying two or three increases more perceptibly. The increase in the percetage of consumers buying 1 characteristic is of 1.99%, compared to about 5% for 2 and 3 characteristics. This, despite the fact that increases in knowledge are decreasing in the number of dimensions (20.04% for 1, 16.90% for 2 and 3.67% for 3). This is possibly explained by both probabilities and by a trickling effect as described in Section 4. At any rate, it is not something we expected ex ante when building the model, and provides with an interesting emergent property of the model. We get back to this in the next sections.

The last three indicators provide with more aggregate information. What we see is an important shift from high to low IKG, and a corresponding inverse for KBG. This means that our model has an important proportion of consumers moving from intention-without-knowledge to knowledge-without-budget situations. Nonetheless, the reduction seen for the IBG indicates that a relevant part of the gains in knowledge translate as gains in consumption.

Parameter modification

We look in some detail at the values of our different indicators. Table 6 in the Appendix 5, shows that:

• All parameters modified (l, S, K, x) have an impact on knowledge indicators (%know1D, %know2D, %know3D) when compared to initial baseline results.

• Parameters' influence on purchase indicators (%buy1D, %buy2D, %buy3D) is also important, particularly when purchases of more than one characteristic are considered, thus providing a confirmation of the result described above.

• All gap indicators (IKG, IBG, KBG) are substantially modified for all parameters, implying that:

=⇒ Knowledge of characteristics effectively circulates reaching a wider proportion of individuals who seek it (thus reducing the intentionknowledge gap).

=⇒ Higher knowledge availability increases both the number of people with knowledge who are able to buy (thus reducing the intentionbehaviour gap) and of those who are unable to (increasing the knowledgebehaviour gap). This, since the extra knowledge will favour purchases only of those with sufficient budget.

• Variations in l and S do not imply a significant difference at time-step 1 when compared with the baseline, meaning that differences in final results are the result of time-based dynamics. Variations in K and x -which create important differences in initial knowledge and purchase capacitysignificantly alter indicators at t = 1 making comparison of results at t = 80 more complex. In the next subsections we look at these in more detail.

3.3 Network-relevant parameters (l and S ). The modification of each l and S -which we have called network-relevantproduce significant results, although knowledge of dimensions by time-step 80 is higher for changes in l, thus giving it a lower intention-knowledge gap. This increased knowledge when l is modified translates also as better results for each of the different indicators, including purchases, intention-behaviour and knowledge-behaviour gap. For purposes of illustration, Figure 2 shows the difference between baseline and a modified l.The results for l and S at t = 80 are shown on Table 4.

Table 4 confirms what we had discussed above: gains in purchases increase as more dimensions are considered, overproportionately to the knowledge gains.

Knowledge threshold and price premiums.

Giving easier initial access to knowledge of the dimensions as well as reducing the price premium distorts the departure point of simulations. This is summarised on the left side of Table 5, where we see that a modified K greatly increases initial knowledge and purchases, and that x increases initial purchases. results at t = 80 are reinforced. The shift is also visible in Figure 3. Also of interest is that the IKG and KBG reach their most extreme values when modifying K, and that purchases of three dimensions seem to depend more on knowledge than on price premiums-meaning that complete goods will be favoured by information more than prices. This is much in line with results for l and S on Subsection 3.3. Lower price premiums also affect the baseline situation from the start, by facilitating purchases at t = 1, although their effect on initial knowledge is minimal. Nonetheless, by t = 80 the higher purchases have strongly encouraged knowledge diffusion-and in turn purchase-thus improving all indicators when compared to baseline.

Concluding remarks

Our work is an effort to use the innovation diffusion framework to study sustainable consumption of (food) products. For this, we have expanded an existing model to include a measure of intention-behaviour gap (caused by price premiums and lacks of knowledge) and the multi-dimensionality that characterises food consumption. We work on three dimensions of characteristics, which can be interpreted as being how healthy, environmentally sustainable and socially fair food purchases are. We first observe that knowledge of characteristics diffuse throughout the simulations regardless of the parameter setup chosen, and that this causes an increase in purchases of sustainable goods and a reduction in the intentionbehaviour gap faced by individuals. In other words, the accumulation of experience and knowledge of sustainable consumption by certain agents creates a shift towards a larger adoption of such practices in the society, and social interaction can be behind reductions in the observed intention-behaviour gap. As it has been noted that "cultural factors are perhaps the most powerful determinants of which food we consume" [Prescott & Graham-Bell, 1995, p. 201], cited in [START_REF] Béné | When food systems meet sustainability -Current narratives and implications for actions[END_REF])], this feature of the model is interesting in terms of its relationship to a real-world phenomenon.

t = 1 t = 80 (∆ baseline ) (∆ baseline ) base K x base K x %know1D
We also see that the modification of certain parameters that affect the dynamics of the network (knowledge spillovers and network density) further reinforce the diffusion of knowledge, the increase in purchases and the reduction of the intention-behaviour gap. This means that societies where interaction on the topic of sustainable consumption is more prevalent, should move quicker towards a widespread adoption of sustainable practices. Conversely, external The most interesting result in our view is that the increase in purchases of additional characteristics is directly proportional to the number of characteristics considered, whereas the increase in knowledge is inversely proportional. In other words, gains in knowledge of characteristics translates into higher purchases when more dimensions are considered. The implications of this are substantial: social interaction and the exchange of knowledge is a more strong determinant of purchases for increasingly complex goods-such as is the case with sustainable food. Given the inherent complexity of sustainable food products, their adoption by a larger proportion of the population will be strongly dependent on how social interaction is exploited in that sense. As Carroll and Fahy has put it, contemporary shopping practices are socially situated activities, to which we add that the more complex the practice, the more socially dependent it will be. Reductions in the intention-behaviour gap will benefit from interaction when an increasing number of dimensions are considered.

Although we did not ex ante expect this result, we have identified one reasons why it may appear, which we attribute to probabilities: the chances of unknowingly purchasing a good will be lower for consumers seeking a larger number of characteristics. Purchases of several dimensions of sustainability will thus be more dependent on the information consumers can gather from their social networks than for those seeking to buy only one (or none).

When modifying parameters, we get further insights into these dynamics. What we see is, first and foremost, a process in which parameters' variations can have important effects on the resulting values of our indicators. We observe that the effect of higher knowledge spillovers and of network degree is, in some respects, as important as that of initial knowledge availability and price premiums. This is an encouraging result, since it can be argued that network-relevant parameters relate to cheaper interventions (creation of discussion forums, other forms of gatherings around food issues) than their knowledge and price premium counterparts (which can be conceived of as messages in the sense proposed by Rogers [2002] or production subsidies). The policy implications of this are important.

In terms of future work, the first avenue we identify is the need to validate our working hypotheses (H1-H2, Section 1.3), to better understand how knowledge of health, environmental and social dimensions in food is held by different actors, and how to more realistically describe its dissemination. It would be welcome to see surveys and field work to put these hypotheses to test and understand how interest is related to knowledge of dimensions, how this knowledge circulates within social groups influencing purchases, and the extent to which we can claim that mental representations of each of these characteristics is common to all individuals. Moreover, empirical work on our model's results, in particular with respect to social interaction being a stronger determinant of purchases of increasingly complex goods would be an interesting endeavour. From a modelling point of view, it would be interesting to explore these results analytically, as well as to generalise to n dimensions of characteristics. Our emerging property of a disproportional influence of knowledge dissemination for a higher number of dimensions is worth further exploration and a generalisation to more dimensions. The model also contains a number of assumptions that could be relaxed in order to increase its descriptive realism, such as fixed preferences throughout the simulations, absence of contrarians or negative influencers, fixed price premiums, etc. Checking how it responds to modifications in these would also be a welcome effort.

The main conclusion our model seems to bring forward is that creating instances of reflection and sharing around the topic of food can be per se a driving force towards a consumption that is respectful of environment, society and health. This, we believe, is a notion worth exploring further.

Appendix: Full table of results

Figure 1 :

 1 Figure 1: Evolution of main indicators (25 runs). The evolution of all indicators is S-shaped.

Figure 2 :

 2 Figure 2: Evolution of main indicators (25 runs). Dashed lines correspond to baseline setup, continuous ones to the parameter l modified by +33%.

Figure 3 :

 3 Figure 3: Evolution of main indicators (25 runs). Dashed lines correspond to baseline setup, continuous ones to the parameter K modified by +33%. The process is shifted right from the beginning.

Table 1 :

 1 Table 1 below sums up the different modifications tested. The full detail of results is shown on Table 6 in the Appendix 5. Different parameters' modifications (±33.3% variations from baseline listed). Each parameter is modified either positively or negatively in the direction that favours knowledge circulation and purchases.

			Parameters modified
		Baseline	l	S	K	x
	l	0.42	0.56 -	-	-
	S	6	-	8	-	-
	K	0.81	-	-0.54	-
	x	0.33	-	-	-	0.22

Table 2 :

 2 Indicators used to analyse simulations.

Table 3 :

 3 Indicators on baseline setup, 25 runs (in brackets, change from t = 1). Knowledge availability increases for 1, 2 and 3 dimensions, as do purchases. IKG and IBG go down, KBG goes up.

Table 4 :

 4 Indicators at t = 80 for +33% modifications of parameters l and S (in brackets, changes from baseline). Compared to baseline, both parameters have a significant impact on all indicators studied.

Table 5 :

 5 Indicators at t = 1 & t = 80 for -33% modifications of parameters K and s (in brackets, changes from baseline). Indicators at t = 1 are strongly modified, making comparisons at t = 80 more complex.

Table 6 :

 6 Baseline results and different parameters modified by 33.3% at t = 1 & t = 80, 25 runs per setup (in brackets, change from baseline).

Data on price premiums for healthy, sustainable and socially fair goods is scarce. As something of an educated guess we add a 33% premium on the presence of one dimension, an

The model allows for the possibility of unknowing pickups, whereby a sufficiently wealthy consumer buys a characteristic he or she has no knowledge of, in which case no experience is accumulated. In this, we do not restrict purchases of goods only to consumers with intention and knowledge of them.

A possible extension of our proposed framework would be to relax this assumption. A

The reader should be informed that the choice of values, which may seem strange at first sight, was done in order to facilitate manipulation. Therefore, we choose $1 as a base price instead of the more common $0 since it allows for a better grasping of percentual variations of price premiums. Likewise, choosing 0.42 or 0.81 for two of the parameters gives round values for our 1/3 rd variations.


S Omid Sayedaghaee

Sergey Prosandeev

Sergei Prokhorenko

Yousra Nahas

Charles Paillard

Bin Xu

L Bellaiche

Domain-wall-induced electromagnons in multiferroics

come

I. INTRODUCTION

Recent observations of magnetoelectric (ME) responses in multiferroics, which characterize cross couplings of electric and magnetic properties, has reawakened the interest of scientists and engineers in this territory of research, particularly in the last two decades [1][2][3][4][5]. Such cross couplings open opportunities to, e.g., control electrical polarization by magnetic field or magnetization by electric field. Switching of polarization by magnetic field can be employed technologically in the design of memory devices and magnetic field sensors whereas switching of magnetization by electric field is useful in designing low energy consumption storage devices in which data is written electrically and read magnetically [5][6][7][8][9][10][11].

The existence of excitations other than phonons, such as magnons, in these systems provides more functionalities such as controlling the spin waves by electric field or by optical phonons [6,12,13]. The dynamical ME couplings of magnons and phonons, known as electromagnons, have also attracted considerable attention due to their role in magnetocapacitance and magnetodielectric effects [14][15][16][17]. Starting in 2006, experimental reports such as exciting spin waves by light in TbMnO 3 [18] and activating magnetic excitations by applying electric fields in YMn 2 O 5 and TbMn 2 O 5 [19] kicked off investigations on the origin and mechanism of electromagnons. In particular, it was also experimentally shown that adjusting Terahertz pulses with electromagnons of 0.7 and 1.8 THz in TbMnO 3 can affect its magnetic structure [20].

In 2008, electromagnons were proposed to be observed for the first time in single crystals [21,22] and thin films [23] of bismuth ferrite BiFeO 3 (BFO). Atomistic methods also predicted electromagnons in antiferromagnetic ferroelectrics [24]. A novel dynamical coupling of magnons with acoustic phonons and optical phonons, coined as electroacoustic magnons, was also suggested in 2019 [25,26]. Electromagnons can have other benefits. For instance, switching the orientation of domains by applying electric fields on multiferroics is typically of the order of milliseconds [20,27], although the speed of ferromagnetic or ferroelectric domain reorientation (which is a key in the design of data storage devices) is typically of the order of nanoseconds [20]. In order to improve the operation speed of devices, one can apply high frequency magnetic fields or optical pulses which increase the functionality speed to the order of femtoseconds to picoseconds [28][29][30] or coherently excite electromagnons [31].

Based on the latter possibility, and also to advance fundamental knowledge, it is timely to address a presently unknown issue regarding electromagnons, namely if they can also be preferentially localized within the walls of ferroelectric domains (or near the middle of the domains) rather than "only" equally vibrate in the entire system. If such preferential localizations exist, novel effects can emerge. Examples include new resonances in the magnetoelectric responses and/or a local control of the dynamical conversion between electric and magnetic properties. In order to resolve the aforementioned issues, we performed molecular dynamics (MD) simulations within the GHz-THz regime of applied magnetic fields in BFO materials possessing ferroelectric domain walls. Fourier analysis of temporal evolution of dipole and magnetic moments along with the calculations of quadratic magnetoelectric response of BFO reveals that the existence of domain walls does result in the emergence of domain-wall-induced electromagnons, which in turn induce large enhancement of ME response around their resonances. These novel quasi-particles were further found to either be preferentially localized near the domain walls or near the center of the domains (depending on their frequency and direction of propagation of their optical phonon and magnon), as a result of scattering by the periodic array made by the domain walls. Such discoveries have the potential to pave the way for designing original devices. For instance, although the capabilities of magnons as the carriers of information in the field of magnonics is quite promising [32], one of the main challenges is to find a way to control the two-dimensional propagation of the magnons in an energy-efficient fashion, especially in small length scales.

A previously proposed approach, which was also realized in experiment recently [33], is to guide the spin-waves in a reconfigurable domain wall as magnonic waveguides. Here, the present work revealing the existence of domain-wall-induced electromagnons strongly suggests that magnonic circuitry can be utilized by activating, via electric fields, the magnons associated with such preferentially localized electromagnons. This should allow to design novel configurable electronic circuits with much less energy consumption.

The article is organized as follows. Section II details the method used here, while results are presented and discussed in Section III. Finally, Section IV provides a summary of this work.

II. METHODS

An effective Hamiltonian framework developed for bismuth ferrite and detailed in Ref. [34] is employed to perform molecular dynamics simulations for the investigation of the temporal evolution of polarization and magnetization as well as the quadratic magnetoelectric response of BFO (which is one of the few materials with magnetoelectric functionalities at room temperature, making it a prime candidate for designing novel devices [10,[35][36][37]) under the application of magnetic fields with various frequencies from 10 to 1300 GHz. Note that the applied time-dependent magnetic field can generate an inhomogeneous electric field which we neglected here in order to merely focus on the ME response rather than the polarization caused by the dielectric response of BFO. In addition, modeling such an electric field was not feasible in the simulations because of using periodic boundary conditions in all three Cartesian directions. Note also that our molecular dynamics simulations are performed at the temperature of 10 K to reduce numerical noises. However, we numerically found that conducting them at room temperature leads to similar results albeit with lower numerical resolutions for the magnetoelectric coefficients [25]. This presently used effective Hamiltonian scheme encompasses three main energy contributions given by Equation 1 [13,34,38]:

E tot = E F E ({u}, {η l }) + E AF D ({u}, {η l } , {ω}) + E M AG ({m}, {u}, {η l } , {ω}) . (1)
In Equation 1, local modes {u} are directly related to local electric dipoles. Strain tensor {η l } includes both homogeneous {η H } and inhomogeneous {η I } strains. The tiltings of oxygen octahedra (also known as antiferrodistortive or AFD motions) are denoted by {ω}, and the magnetic moments associated with iron ions are denoted by {m}. Also, in Equation 1, the first term includes the energy of local modes and elastic deformations that can be written as [39,40]

E F E =E self ({u}) + E dpl ({u}) + E dpl ({u}) + E short ({u}) + E elas ({η l }) + E int ({u}, {η l }) , (2)
in which E self is the self-energy of local-modes, E dpl denotes the long-range dipole-dipole interaction, E short is the short-range interaction of soft modes, E elas denotes the contribution of elastic energy, and E int is the interaction of local modes and strain.An extensive overview and explanation on the energy terms of Equation 2 is provided in Ref. [40].

The second term of Equation 1, whose energy is only related to AFD motions and their couplings with local modes and strains, is given by

E AFD = i κ A ω 2 i + α A ω 4 i + γ A ω 2 ix ω 2 iy + ω 2 iy ω 2 iz + ω 2 ix ω 2 iz + ij αβ D ijαβ ω iα ω jβ + i αβ C lαβ η l (i)ω iα ω iβ + ij αβ K ij,αβ u j,α ω i,α ω j,β + i αβγδ E αβγδ ω iα ω iβ u iγ u iδ , (3)
where indices i and j numerate the lattice sites and α, γ, ν, and δ are the Cartesian components. In Equation 3, the first term shows the contribution of self-energy of the tilting of the oxygen octahedra while the second term shows the contribution of the short range interactions of AFD motions. The third term shows the coupling of the AFD motions with strain. The fourth and fifth terms show the contributions of trilinear and biquadratic energies coupling the local modes and AFD distortions. More details about AFD energy contributions can be found in Ref. [41]. The third term of Equation 1shows the magnetic moments and the coupling of local modes, strain, and antiferrodistortive motions with magnetic moments as given in the following

E M AG = ijαγ Q ijαγ m iα m jγ + ijαγ S ijαγ m iα m jγ + ij,αγνδ E ij,αγνδ m iα m jγ u iv u iδ + ij,αγνδ F ij,αγvδ m iα m jγ ω iv ω iδ + ijl,αγ G ijl,αγ η l (i)m iα m jγ + ij L ij (ω i -ω j ) • (m i × m j ) . (4)
In Equation 4, the first term indicates the dipolar interactions of the magnetic moments.

The second term shows the short-range magnetic exchange coupling. The third, fourth, and fifth terms indicate the contribution of local modes, AFD motions, and strain in magnetic exchange interaction, respectively. The sixth term shows the contribution of the tilting of oxygen octahedra through Dzyaloshinskii-Moriya interaction [42][43][44]. In addition to the energy contributions discussed above, if external magnetic and/or electric fields are applied, their contribution to the energy asi m i • H for magnetic field andi P i • E for electric field should be added to the total energy terms. The electrical polarization, P , can be computed from its proportionality to the local modes via P = Z * u V 0 in which Z * is Born effective charge obtained from first-principle calculations and its value for BFO is

Z * = 5.68 C [45]
, while V 0 is the volume of the unit cell.

To predict the dynamical properties of BiF eO 3 within the Effective Hamiltonian framework used here, Newtonian equations of motion are applied for theionic degrees of freedom (i.e. local modes, AFD, strains) discussed in Equation 1where the forces were computed by taking partial differential of the total energy terms similar to the techniques used in Refs.

[ [46][47][48]. Moreover, as demonstrated in Ref. [13], the dynamics of the magnetic moments are considered through the Landau-Lifshitz-Gilbert (LLG) equation [49] as indicated in the following:

dm i dt = -Γm i × B i eff (t) + b i fl (t) -Γ λ |m i | (1 + λ 2 ) m i × m i × B i eff (t) + b i fl (t) . (5)
Here, B i eff is the effective magnetic field which acts on the ith magnetic moment and is given by B i eff = -∂E tot /∂m i , while Γ and λ are gyromagnetic and damping coefficients, respectively. b i fl is is a fluctuation field which is also acting on the ith magnetic moment.

More details on the dynamics of magnetic properties calculations are given in Ref [13]. that is applied along the [ 211] direction and which is a combination of a dc field with an ac field of ω frequency. We chose H dc = 245 T and h ac = 61.2 T in order to get good statistics on the quadratic ME coefficients since, for such fields, linear magnetoelectricity has a small effect with respect to the quadratic magnetoelectric coefficients in BFO (see Refs. [36,38] and the Supplementary Information of Ref. [26]). However, and as shown in the Supplementary Material, the simulations were performed at lower magnetic fields to make sure that the observed results are still valid at magnetic fields that can be applied in practical applications.

Note that, in the entire supercell, the initial magnetization is along the [ 211] direction, i.e., along which the magnetic field is applied, while the G-type antiferromagnetic vector is along the perpendicuar [01 1] direction. Furthermore, ω is fixed during each simulation, but is allowed to vary from 10 GHz to 1300 GHz between different computations. The upper limit of investigated frequencies is intentionally chosen such as to avoid resonating with the A 1 phonon frequencies (known to occur at much higher frequencies of the order of ∼4300 and ∼7000 GHz [25,52]), but rather to focus on the possibility of having electromagnons associated with domain walls and thus occurring at lower frequencies (note that we do see these A 1 phonon frequencies in our simulations but do not report them here).

Note also that the homogeneous strain has been recently predicted to combine with optical phonons and magnons to create a new quasi-particle, coined electroacoustic magnon, that can result in the emergence of resonance in the magnetoelectric response [26]. Since we are especially interested here in revealing domain-wall-induced electromagnons (i.e., mixing of optical vibrations with magnons in a multidomain structure), we perform two kinds of MD simulations: one for which the homogeneous strain is clamped (implying that no electroacoustic magnon can exist) versus another for which such strain is fully relaxed. The inhomogeneous strain, on the other hand, is considered in the effective Hamiltonian and can vary during the course of all types of simulations conducted here.

Note also that, in multiferroics, the electrical polarization under the influence of applied electric and magnetic fields can be obtained by differentiating the expansion of the free energy with respect to electric field [5]. In the absence of external applied electric fields and when neglecting the linear magnetoelectric coefficient, the temporal evolution of the polarization as a function of our applied magnetic fields is given by

P (t) = P (t=0) + β(0, ω)H dc h ac sin(ωt) + β(ω, ω) 1 2 h 2 ac sin 2 (ωt). (6)
where β(0, ω), on which we primarily focus, is the quadratic ME coefficient coupling the ac and dc fields, while β(ω, ω) couples the ac field with itself and generates a second harmonic at 2ω.

The real and imaginary parts of β(0, ω) can thus be obtained by

β ′ (0, ω) = 2 τ τ 0 (P -P (t=0) )sin(ωt)dt H dc h ac , (7)
and

β ′′ (0, ω) = 2 τ τ 0 (P -P (t=0) ) cos(ωt)dt H dc h ac , (8)
in which τ is the overall time of the simulation. More details about the derivation of equations ( 7) and ( 8) are provided within the Supplementary Material.

B. Frequency Analysis

Let us first consider the case of the aforementioned multidomain structure when the homogeneous strain is clamped. During and at the end of the simulation, the overall polarization is along the [100] direction, since the y-and z-components of the dipoles in the left and right domains cancel each other. Furthermore, the direction of the magnetic moments during and at the end of the simulation is found to be along the applied dc magnetic field, i.e., [ 211], in both domains.

We then perform Fourier analysis of the temporal evolution of this electrical polarization and resulting magnetization along different directions (in order to, e.g., reveal the directions of vibration of phonons and magnons that are dynamically coupled to form electromagnons).

Note that here we are dealing with two types of frequencies, one being the frequency of applied magnetic field which is denoted by ω, while the other is the frequency spectrum that we get as a result of Fourier analyses of the dynamical evolution of the polarization or magnetization under the applied field. For instance, the frequency spectra of the electrical

[100]-component of the dipoles (direction of polarization, as well as normal to the ferroelectric domain wall) and the perpendicular, [011]-component, when the frequency of the applied magnetic field, ω, is 10 GHz are demonstrated by red solid lines in Figures 2a and2c, respectively, and up to about 1300 GHz. Apart from the frequency of applied field and its second harmonic at 10 and 20 GHz, two peaks are observed at ∼390 and ∼1150 GHz in the spectrum of the [100]-component (Figure 2a). Also a peak is observed at ∼700 GHz for the

[011]-component (Figure 2c).

In addition, the Fourier analyses of the magnetization were performed for the [ 211]component which is also the direction of the magnetic moments and the applied magnetic field (Figure 2b) as well as the [01 1]-component which is also a direction that is perpendicular to both polarization and magnetization (Figure 2d). Since the peaks at ∼390 GHz In other words, one can conclude that the 390, 1150, and 700 GHz modes are domain-wall-induced electromagnons (to be denoted as DWI electromagnons in the following). Such observations are also consistent with previous findings that vibrations of the order of 360, 700 and ∼ 1200-1800 GHz caused by domain wall motions were predicted for the case of BFO adopting domain structures [50,53] (note that other domain wall activities, but with lower frequencies, were also reported in Refs. [54,55] but the limitation of molecular dynamics methods below the range of few GHz prevents us from determining them).

Figures 2a, 2b, 2c, and, 2d also report the case of a multidomain but for unclamped homogeneous strain, which is represented by dashed blue lines. In this latter case, three additional peaks are observed at 70, 100 and 270 GHz with respect to the case of clamped multidomain, for both the Fourier transform of the temporal evolution of the polarization and magnetization. These three peaks, therefore, indicate that strain can mix with optical lattice vibrations and magnons to produce new particles -which is in-line with the recent prediction of the so-called electroacoustic magnons for unclamped monodomain in Ref. [25]. Note that these electroacoustic magnons at 70 and 270 GHz in our multidomain involve vibration of the polarization along the normal of the domain wall (see Fig. 2a) while that at 100 GHz is better seen in the oscillation of the electric dipoles perpendicularly to that normal (cf Fig. 2c). Moreover, the magnon associated with these electroacoustic magnons at 70, 100 and 270 GHz have vibrations along both [ 211] (that are stronger) and [01 1] (that are weaker), as demonstrated by Figures 2b and2d.

C. Magnetoelectric Response

Let us from now on focus on results obtained for the case of unclamped homogeneous strain in the multidomain structure, since we now know the origins of all the peaks appearing in Figure 2. The real and imaginary parts of the quadratic magnetoelectric coefficient, β(0, ω), for polarization lying along the [100] and [011] directions, are then computed using Equations ( 7) and ( 8) at each frequency of applied magnetic field from 10 to 1300 GHz in that case. They are shown in Figures 3a and3b, respectively. The calculated values are in agreement with experimental measurements of β of bulk BFO at small frequencies. For instance, the obtained value for the real part of β(0, ω) along the [100] direction by our MD simulations when the frequency of applied magnetic field is 20 GHz and the homogeneous strain is unclamped, is 2.457 × 10 -8 C/m 2 T 2 = 0.388 × 10 -19 s/A which remarkably agrees with the observation of β 311 = 0.3 × 10 -19 s/A reported in Reference [56] for the static case.

Also, in measurements that were performed in high magnetic fields [57], the quadratic ME coefficient was found to be

β 122 = β 212 = β 221 = -β 111 = -5 × 10 -9 C/m 2 T 2 which is in
reasonable agreement with the value of -1.004 × 10 -9 C/m 2 T 2 that we have computed along the [011] direction at ω=20 GHz. The ME coefficients exhibit resonances when the frequency of applied magnetic field, ω, approaches the frequency of the electroacoustic-magnon-like particles mentioned above, that are 70, 100 and 270 GHz here, as previously found for the case of a monodomain, for which electroacoustic magnons have resonant frequencies at 90 and 267 GHz [25,26].

Moreover, another important novelty of the present work is that resonances also happen when ω reaches the frequency of the three domain-wall-induced electromagnons. For instance, along the [011] direction, a nicely seen resonance takes place at ω=700 GHz, while, along the [100] direction, one resonance is rather sharp at 390 GHz while the other one around 1150 GHz is broader and even seems to be made of several peaks -which is reminiscent of the broad response found in Ref. [53] around these frequencies for the domain wall motions. These results are also consistent with the dielectric and domain-wall-related modes A and B reported in Ref. [50] for which mode A around 12 cm -1 ≃ 360 GHz is sharp while mode B around 50 cm -1 ≃ 1500 GHz is broader. Our DWI electromagnons therefore provide a path to engineer devices based on the functionality of the domain walls.

Although domain walls have already been employed in spintronics and for designing novel devices [58][59][60], their presently discovered role in optimizing ME responses make them also fascinating for magnetoelectric-based devices. As a matter of fact, our results show that one can tune the frequency of the applied magnetic field with the frequency of the domain wall vibrations (DWI electromagnons) to maximize the ME response.

D. Degree of localization of the domain-wall-induced electromagnons

In order to investigate the degree of localization that these DWI electromagnons possess, a layer-by-layer analysis is performed at different planes that are parallel to the domain wall.

As the dimension of the supercell is 20 × 10 × 10, each 10 × 10 cross section of the supercell along the [100] direction is defined as a layer for which the average of magnetization and polarization can be computed from the MD simulation results. Also, the quadratic ME coefficient of each layer can be obtained via equations ( 7) and ( 8) but for which the overall polarization is replaced by the average polarization of each layer. The layers are indicated by their plane numbers, with the planes 10 and 20 corresponding to the domain wall layers whereas planes 5 and 15 correspond to the center of each domain. Since the magnetoelectric coefficients β(0, ω) are calculated at each frequency of applied magnetic field for frequencies up to 1300 GHz, in each plane that is parallel to the domain wall, the absolute value of the real part of magnetoelectric coefficient, |Re(β(0, ω))| (which is related to real-space location of optical phonons) is used to investigate the preferential location of the optical phonons.

On the other hand, for each frequency of applied magnetic field, the Fourier analysis of magnetic moments provides a spectrum of the frequencies that are present in the vibration of the magnetic moments. Therefore, the Fourier intensity of the magnetization which is related to real-space location of magnons is used to investigate the preferential location of magnons. Here we used the Fourier analyses when the frequency of applied magnetic field is 10 GHz.

Figures 4a and4b provide an example of the results on such analysis for the optical phonon (Figure 4a) and magnon (Figure 4b) associated with the frequency of 700 GHz.

One can see that the maximum value of the resonance peak inside the domain (plane 5), which is 1.732 × 10 -8 C/m 2 T 2 , jumps to 4.393 × 10 -8 C/m 2 T 2 at the domain wall (plane 10)

showing that the ME coefficient is increased by 154%. Such comparison clarifies that not only the ME coefficient exhibits a resonance when the frequency of applied magnetic field is compares the intensity of Fourier analyses peaks around 700 GHz for magnetization along the [01 1] direction in each layer. The Fourier intensity of vibrations of magnetization at ω=700 GHz is stronger at the domain wall. Consequently, one can conclude that the DWI electromagnon at ω=700 GHz can be seen as waves that adopt maximum at the domain wall, that is both the optical phonon and magnon forming such electromagnon preferentially localize at this ferroelectric domain wall. Such domain-wall-induced electromagnon's wave behavior and preferential localization are also found for the frequency of 1150 GHz (but for optical phonon vibrating along [100] and magnons oscillating along [ 211]). As a matter of fact and as demonstrated in Figures 5(a) and (b), both the optical phonon and magnon of the domain-wall-induced electromagnon at 1150 GHz adopt wave-like behavior with maxima occurring in the vicinity of the domain walls (near plane 10). On the other hand, the DWI electromagnon having a resonant frequency of 390 GHz, while still possessing wave-like behavior, prefers to possess maxima occurring, both for its optical phonon and magnon, near the middle of the domains (near planes 5 and 15) rather than at the domain wall, as demonstrated in Figures 6a and6b. Note that the Supplementary Material also reports information about localization/delocalization of the electroacoustic magnons and some results on inhomogeneous strain analysis.

E. The origin of the domain-wall-induced electromagnons

To have a clear understanding of the origin of the 390 GHz, 700 GHz and 1150 GHz electromagnons induced by the domain structure, we have analyzed the spectral energy density (SED) associated with local dipoles (Φ P ) and inhomogeneous stain (Φ η ) when the homogeneous strain was unclamped and for a 40×10×10 supercell. In order to conserve the periodicity of the domain structure, the supercell is split into four alternating domains of equal width. Note that a larger supercell (i.e. 40×10×10 as compared to 20×10×10) is used to increase the wave-vector q resolution. Note also that Φ X , with X=P or η, describes the wave-vector q and frequency ν resolved power of the kinetic energy spectrum for the vector degrees of freedom X. It is computed [61] as

Φ X (q) = 1 N τ α=x,y,z Ẋα q,ν 2 , ( 9
)
where Ẋα q,ν denotes the Fast Fourier Transform coefficient of the α component of the velocity Ẋα . N and τ denote the total number of unit cells and the molecular dynamics trajectory time, respectively. The poles of Φ P provide the spectrum of excitations independently of the mode's polarization vector. Specifically, Φ P and Φ η can be used to obtain dispersions of optical and acoustic phonons, respectively.

Practically, we use a single 0.5 ns molecular dynamics trajectory computed for this N = 40 × 10 × 10 supercell. The integration time step is taken to be equal to 0.5 fs and the first 10,000 steps of the trajectory are not taken into consideration. The velocities are sampled each 10 MD steps for a subsequent Fast Fourier transforms to the wave-vector and frequency domain.

Figure 7a and 7c shows the computed low-frequency spectra of polar excitations along the [100] pseudo-cubic direction (i.e., different panels on the left/right side correspond to different q-points along that direction) for mono-and multi-domain states, respectively. The spectrum for the mono-domain configuration (panel a of Fig. 7) clearly shows the dispersion of three electromagnon branches labeled ( 1)-( 3). From the analysis of Φ η , we also find that these branches follow the dispersion curves of acoustic phonons (not shown here). Branches

(1) and ( 2) match the two transverse acoustic (TA) modes, while the third branch is mixed with the longitudinal acoustic (LA) mode.

Interestingly, the three described monodomain branches can be also clearly seen in the spectrum of the multidomain state (see panel c of Fig. 7). However, the multidomain spectrum additionally features multiple "satellite" peaks not present in the monodomain case (panel a of Fig. 7). These features most likely stem from the elastic scattering of the three mono-domain electromagnons by a periodic and static array of domain walls. In such hypothetical scenario, a single scattering event would conserve the mode frequency and lead to a shift of its momentum q x by ±∆q x . The possible shifts

∆q x = 2 20 π a 0 , 6 20
π a 0 , 10 20
π a 0 (10)
correspond to the dominant Fourier components of the local mode distribution, or, in other words, are determined by the periodicity of the domain structure. Here a 0 indicates the 5-atom lattice constant.

The simplest single scattering mechanism described above allows to explain the appearance of the predicted domain-wall-induced electromagnons at q x = 0. Namely, the emergent zone-center excitations at 390 GHz and 700 GHz can be matched with monodomain modes Such hypothesis about the origin of the domain-wall-induced electromagnons is further supported by visualizing the local mode vibrations at ν 0 =390 GHz, 700 GHz and 1150 GHz. To perform this analysis, we first note the deviation of the average local mode value δu(x, t) = u(x, t) -⟨u(x)⟩ within each x-layer of the supercell and at each MD step t from its time average ⟨u(x)⟩. Then, we perform the time-frequency Fourier transformation δu(x, t) → δu(x, ν) that allows us to look at how the layer-averaged local modes change along the [100] direction and only at the frequencies of interest (i.e., at these ν 0 ). We perform the described procedure using the generated MD trajectories for the multidomain case which technically corresponds to isolating the excitations pattern at a specific frequency in Figure 7c. Figure 9 shows how the resulting complex amplitudes δu α (x, ν 0 ) and phases θ α = arg(δu α (x, ν 0 )) change with x for each Cartesian component α of the local modes.

As expected from the scattering picture, the amplitudes of local modes displacements at ν 0 =390 GHz (Fig. 9a) and ν 0 =700 GHz (Fig. 9b) have a pronounced q x = 2 20 π a 0 character with half of the wave period falling within each of the domains (Figure 7b). At the same time, and as also consistent with Figure 8c, the amplitudes δu α (x, ν 0 ) at ν 0 =1150 GHz (see Figure 9c) have a different, three half-period, character rather described by q x = 6 20 π a 0 (Figure 7b). This latter periodicity is more clearly visible in the dependence of the phase on x at ν 0 =1150 GHz (see Figure 9f).

It is also interesting to note the markedly distinct behavior of the phase plots at ν 0 =390 GHz and ν 0 =1150 GHz from the complex phase distribution at ν 0 =700 GHz.

Indeed, the two former feature a constant (flat) phase θ α (x) dependencies on x within each domain (Figures 9d andf). Such distribution along with a spatially modulated amplitudes (Figures 9a andc) are unique signatures of standing waves. Notably, the standing-wave character (confinement) of these modes also goes in-line with the TA nature of the original monodomain branch (1) giving rise to these modes. Indeed, similar situations were previosuly reported in references [62,63] describing the confinement of TA acoustic phonons in semiconductor superlattices. In contrast, the ν 0 =700 GHz vibration pattern rather corresponds to a propagating wave as the latter typically yield a linear phase dependence on the coordinate (see Figure 9 e). This observation also goes in-line with our prediction of the common origin of the ν 0 =390 GHz and ν 0 =1150 GHz electromagnons distinct from that of the ν 0 =700 GHz domain-wall-induced mode.

y = A + C|sin(q(x -x 0 ))| (11)
where q = 2 20 π a 0 is the scattering vector for that frequency, A is the amplitude of the zonecenter mode and C can be seen as an x-component amplitude of the confined electromagnon at q. One can clearly see that such function fits well the data of Fig. 6a, which confirms our proposed concept of scattering. Similarly, the values of the Fourier intensity peaks of magnetization along the [ 211] direction at 390 GHz were fitted by a function of the form

y = A + C ′ |sin(q(x -x 0 ))| + D ′ sin(q(x -x 0 )) ( 12
)
where q is still equal to 2 20 π a 0 , C ′ is related to the x-component contribution of the electromagnon as in Eq.( 11), while D ′ originates from its y-and z-components. Once again, this function fits well the data of Fig. 6b, which is another confirmation that scattering is the main source of the domain-wall-induced electromagnons. It should be noted that the difference between the data and the fit at domain walls comes from the modes localized on the domain walls.

IV. SUMMARY

In summary, we report here the existence and explain the origin of domain-wall-induced electromagnons, and discover that they enhance the magnetoelectric response of the material to the application of external magnetic fields. Obtaining such results for bisumth ferrite not only implies that ME enhancement can be achieved at room temperature, but also infers that such enhancement can be generalized to other multiferroics possessing multidomains. Interestingly, the optical phonon and magnons forming such domain-wall-induced electromagnons behave as waves for which the maxima jointly occur in some specific parts of the multidomain, namely near the domain walls for the 700 and 1150 GHz frequencies versus near the middle of each domain for the 390 GHz frequency. The presently discovered domain-wall-induced electromagnons may thus be of some advantages for, e.g., the design of (i) novel ultra-fast energy-efficient data storage devices via local probes and (ii) and of electrically-driven configurable magnonic circuits taking advantage of the spatial preferential localization of both the phonons and magnons making the domain-wall-induced electromagnons. Each sub-panel in panels (a) and (c) shows the dependence of the logarithm of the spectral energy density Φ P on frequency at a single wavevectors q = (q x , 0, 0). The labels (1), ( 2

  III. RESULTS AND DISCUSSIONA. Polarization, Magnetization, and Magnetoelectric CoefficientHere and as shown in Figure1(a), we investigate BFO with two domains separated by a 109 • domain wall through a periodic 20 × 10 × 10 supercell (For comparison, Figure1(b)displays a monodomain). We have investigated 71 • and 180 • domain configurations as well, and the results of our analyses are consistent with 109 • configuration. However, since 109 • domain walls are the ones that are found most frequently in BFO bulk materials and the most seen ferroelastic domain walls in single crystal samples of BFO, the results for 109 • configuration are shown and discussed here[50,51]. In Figure1(a), the left half-side along the x-axis (which lies along the pseudo-cubic [100] direction), is made of electric dipoles being initially along[111] while the right half-side is formed by electric dipoles initially lying along[1 11 ]. This supercell is subjected to a magnetic field, H = H dc + h ac sin(ωt),

  and ∼1150 GHz are present in the Fourier analyses of both [100]-component of polarization and [ 211]-component of magnetization, these peaks represent electromagnons for which the optical phonon has a component of its vibration along the direction that is perpendicular to the domain wall and the magnon has a component of its vibration along the magnetization. (It should be noted that, we found such electromagnons for the 71 • domain configuration at ∼455 GHz and ∼1390 GHz. Also, for the 180 • domain configuration, such electromagnons are found to be at 240 GHz and 1430 GHz). Similarly, the peaks at ∼700 GHz seen in Figs 2c and 2d characterize an electromagnon made of an optical phonon that has a vibration direction with a component lying along the [011] direction and of a magnon for which the vibration has a finite component along [01 1] (note that [011] and [01 1] are both lying in the plane of the domain wall but are perpendicular to each other). Interestingly, and as verified by black solid lines in Figures 2a to 2d, these three peaks disappear when the simulations are performed for the case of monodomain structure (whether the homogeneous strain is clamped or unclamped).

( 1 ) 20 π a 0 20 π a 0 20 π a 0 by ∆q x = -6 20 π a 0

 120020020200 and (3) at q x = 2 that are scattered by ∆q x = -2 (see panels a and b of Fig.8), while the homogeneous 1150 GHz mode stems from scattering of mode (1) at q x = 6 (see panels a and c of Fig.8).

20 π a 0

 200 Finally, from Figures 9 it can be further noted that for all of the considered frequencies the vibration patterns are neither the classical standing nor classical propagating waves. Indeed, the amplitude plots in Figures9clearly indicate that all modes are mixed with Gamma point vibrations while the amplitude deviations at the domain walls indicate the presence of domain-wall localized modes. Moreover, the phase velocity of the propagating electromagnon at ν 0 =700 GHz clearly increases upon approaching the domain walls (Figure9e) while, at the same time, its amplitude experiences a drop. Such change of the phase velocity can be equivalently described as a mixing of the propagating mode with a domain-wall localized mode and a standing wave at q x = 2 20 π a 0 . All of the described deviations are due to the aforementioned scattering which couples phonons with wavevectors differing by the multiples of 2 due to which all of the three modes have a mixed character. To complete the picture, we summarize the information presented in Figure9by plotting the instantaneous snapshots of the calculated profiles of dipole vibrations at ν 0 =390 GHz, 700 GHz and 1150 GHz in Figure10.Based on the gathered insights, we now fit (as an example) the values of the Fourier intensity peaks of ME coefficient, β(0, ω), along the [100] direction and at 390 GHz by a function of the form

10 FIG. 2 :FIG. 3 :

 1023 FIG. 2: (Color online) The Fourier analysesThe Fourier analyses of the polarization (panel a is for polarization along the [100] direction while panel c is for polarization along the [011] direction) and magnetization (panel b is for magnetization along the [ 211] direction while panel d is for magnetization along the [01 1] direction) for the case of monodomain when the strain is unclamped (black) and for the case of multidomain when the homogeneous strain is unclamped (blue) and clamped (red). Here, the frequency of the applied magnetic field is 10 GHz. Note the change of vertical scale between Panels (b) and (d)

FIG. 4 :

 4 FIG. 4: (Color online) The layer-by-layer analysis for the domain-wall-induced electromagnons at 700 GHz Layer-by-layer comparison of the intensity of the (a) resonance peaks of the real part of ME coefficient, β, along the [011] direction when the frequency of applied magnetic field is 700 GHz; (b) Fourier peaks around 700 GHz for magnetization vibrations along the [01 1] direction; For panels (b) the frequency of applied magnetic field is 10 GHz.

FIG. 5 :

 5 FIG. 5: (Color online) The layer-by-layer analysis for the domain-wall-induced electromagnons at 1150 GHz Panel (a) shows the magnetoelectric coefficient, β, along the [100] direction when the frequency of applied magnetic field is 1150 GHz;. Panels (b) shows the Fourier spectrum of magnetization, along the [ 211] direction for 1150 GHz. For panels (b) the frequency of applied magnetic field is 10 GHz.

5 FIG. 6 :

 56 FIG. 6: (Color online) The layer-by-layer analysis for the domain-wall-induced electromagnons at 390 GHz Layer-by-layer comparison of the intensity of the (a) resonance peaks of the real part of ME coefficient, β, along the [100] direction when the frequency of applied magnetic field is 390 GHz; (b) Fourier peaks around 390 GHz for magnetization vibrations along the [ 211] direction. For panel(b) the frequency of applied magnetic field is 10 GHz.

FIG. 7 :

 7 FIG. 7: (Color online) Low-frequency spectrum of polar excitations for monodomain and multidomain states.

  ) and(3) indicate the relevant electromagnon branches. The value of q x for each row of sub-panels is indicated in panel (b). Panel (b) also contains a schematic representation of a 40 × 10 × 10 supercell with four alternating domains (blue and white regions) along with examples of the q wave modulations.

FIG. 8 : 20 π a 0 and q x = 6 20 π a 0

 820200 FIG. 8: (Color online) Origin of the zone-center electromagnons. Panel (a) shows the logarithm of the optical phonon SED at q x = 0 for the multidomain state. Panels (b) and (c) shows the logarithm of the optical phonon SED for the monodomain state at q x = 2 20

2

 2

FIG. 10 :

 10 FIG. 10: (Color online) Instantenous distribution of the local mode deviation vectors for the three considered domain-wall-induced electromagnons. The arrows indicate the local mode deviations with red (purple) color corresponding to vector with positive (negative) z components. Electromagnons at 390 GHz and 1150 GHz originate from the same monodomain electromagnon branch and have predominantly confined (standing-wave) character. The 700 GHz domain-wall-induced mode features a pronounced propagating wave nature.

Panels a, b and c (d, e and f) show the changes of the complex amplitude (phase) of the local mode vibrations along the [100] direction at 390 GHz, 700 GHz and 1150 GHz, respectively.

ACKNOWLEDGEMENTS The authors are thankful to the DARPA Grant No. HR0011-15-2-0038 (under the MATRIX program). S. Prosandeev also thanks ONR Grants No. N00014-17-1-2818 and N00014-21-1-2086. B.X. further acknowledges the financial support from National Natural Science Foundation of China under Grant No. 12074277, the startup fund from Soochow University, and the support from Priority Academic Program Development (PAPD) of Jiangsu Higher Education Institutions. S. Prokhorenko and Y.N. thank DARPA Grant No. HR0011727183-D18AP00010 (under the TEE Program). The authors are also thankful for the support of MRI Grant No. 0722625 from NSF, ONR Grant No. N00014-15-1-2881 (DURIP), as well as a Challenge grant from the Department of Defense, and also acknowl-