Intonation of Greek in contact with Turkish: a diachronic study

Abstract Asia Minor Greek (AMG) speakers cohabited with Turkish speakers for eight hundred years until the 1923 Lausanne Convention, which forced a two-way mass population exchange between Turkey and Greece and severed their everyday contact. We compare the intonation of the continuation rise tune in the speech of first-generation AMG speakers born in Turkey with three subsequent generations born in Greece. We examine how long contact effects in intonation persist after contact has ceased, through comparison of the f0 patterns in four generations of AMG speakers with those of their Athenian Greek- and Turkish-speaking contemporaries. The speech of the first-generation of AMG speakers exhibits two patterns in the f0 curve shape and time alignment of the continuation rises, one Athenian-like and one Turkish-like. Over subsequent generations use of the latter diminishes, while the Athenian pattern becomes more frequent, indicating intergenerational change.

Peninsular Spanish) while Cuzco speakers who are Spanish-Quechua bilinguals do not consistently differentiate between broad and contrastive focus, something that is attributed to the influence of Quechua (O'Rourke, 2012).
Sometimes contact gives rise to novel patterns attested in neither language. Lekeitio Spanish, a contact variety of Basque and Spanish, displays mixed properties of Basque and Spanish intonation, that is, accents in broad focus declaratives are not falls, as in the H*+L of Lekeitio Basque, but rises as in the L*+H or L+H* of Madrid Spanish. On the other hand, the peak alignment of the H tone in the L*+H rises does not appear after the stressed vowel as in Spanish but within the stressed vowel, as in Lekeitio Basque (Elordieta & Calleja, 2005).
In addition, bilingual speakers have been shown to codeswitch in prosody. Queen (2012) reported codeswitched productions in the polar question tunes of German-Turkish bilinguals: when the Turkish question particle /-mI/ was inserted into a German matrix utterance, the question particle was produced with a falling (H*L) Turkish f 0 contour.
Evidence is still lacking on whether and how long such effects persist after contact has ceased, as few studies exist on the diachronic change of intonation (Hualde, 2004). Colantoni and Gurlekian (2004) argue that the alignment patterns still found in Buenos Aires Spanish declaratives almost sixty years after the end of the immigration wave of Italians are due to the covert prestige of Italian. Barnes and Michnowicz (2015) report that the contact between Veneto Italian and Mexican Spanish in Chipilo, Mexico resulted in Italian-like patterns of prenuclear peak alignment in broad focus declaratives. They maintain that the influence of Veneto Italian has persisted because speakers identify with Veneto and are proud of their language and heritage.
Similarly lasting effects of contact have been documented in the direction from a dominant to a nondominant language. Bullock (2009) discusses prosodic patterns in French-English contact in a minority variety of French spoken in Frenchville, Pennsylvania since the 1830s. In the last two remaining fluent heritage speakers in the village, there was adoption of additional structures of the dominant language (i.e., pitch accents used in English but not in French) in the intonation of these French heritage speakers, which have distinct pragmatic uses.
In this paper, we trace the trajectory of intergenerational prosodic change in four generations of speakers of Asia Minor Greek (AMG). AMG patterns are compared with those of four generations of their Turkish-speaking contemporaries to bring to light any similarities that may have remained after the end of AMG-Turkish contact and four contemporaneous generations of Athenian speakers, since AMG speakers have now been living in Greece for almost a century (see Baltazani, Przedlacka, & Coleman, 2020).
Our approach differs from those cited above in several respects. First, we make a three-way comparison analyzing benchmark data from Athenian and Turkish in addition to AMG, whereas the papers cited above (with the exception of Baltazani et al., 2020) only analyze data in the contact variety, relying on previously published analyses of the two source languages. Second, we investigate data from four generations, whereas most of the papers cited above use data from just one or two. Third, whereas most studies use small corpora of speech elicited in a laboratory setting, we employ a large corpus of mostly spontaneous and semispontaneous speech because of the stronger influence of the source language on a contact variety in natural than in controlled laboratory speech (Barnes & Michnowicz, 2015).
Background on Athenian, Turkish, and AMG Athenian, the Modern Greek variety spoken in Athens, is the standard used for official purposes, in education, and the media. Athenian is thought not to be distinctly marked by any single traditional Greek dialect (Trudgill, 2003:48). Rather, extensive internal migration from various rural parts of the country to Athens between 1950 and 1980 resulted in today's amalgamation 1 of different varieties brought by migrants (Allen, 1986).
Despite reports on the leveling influence on Turkish of the Istanbul standard used in mass media and the education system since the 1930s (Johanson, 2001:16), there are descriptions of phonological, morphological, and syntactic dialectal features for several regions in Turkey, including the Black Sea (Johanson, 2001), Marmara, Aegean, and Central Anatolia regions (Karlık & Akbarov, 2015). To our knowledge, no descriptions exist about dialectal intonation differences. As we base the analysis of Turkish patterns on recordings from several geographic regions, we conducted a statistical comparison of the intonational parameters in these dialectal recordings (in 3.1), which revealed that they are similar, so for all subsequent analyses we grouped all Turkish data together.
AMG is a variety of Modern Greek previously spoken in Asia Minor by people whose ancestors had been part of the Byzantine Empire until the eleventh century conquest by Seljuk Turks. (For a detailed, well-documented linguistic history, see Karatsareas [2011:11-63] and sources cited there.) The Greek-speaking peoples in Asia Minor mixed culturally and linguistically with Turks for eight centuries in a society where Turkish was the dominant language. There is historical evidence of some linguistically Turkicized Greek communities and of others retaining Greek as their first language (Karatsareas, 2011:16-19). In 1923, the cohabitation was forcibly ended due to the Convention Concerning the Exchange of Greek and Turkish Populations. Approximately 1,300,000 Orthodox Christians of the Ottoman Empire were relocated to Greece, in exchange for roughly 360,000 Muslims from Greece to Turkey (Blanchard, 1925).
AMG is a heritage variety (Polinsky, 2018) only spoken at home and not dominant at the national level. It survives mostly in northern Greece (Katsapis, 2011:71) where currently there are second-, third-, and fourth-generation speakers of AMG who, unlike the first generation, did not grow up in a Turkish-speaking environment. The speakers are bi-or tridialectal, in AMG, local varieties, and Athenian Greek as part of their linguistic repertoire (Janse, 2009;Karatsareas, 2011). First-generation speakers in this study were born in Cappadocia prior to the 1923 Convention while speakers of subsequent generations are descendants of the Cappadocian refugees and were born in Greece. Refugees from Asia Minor made up more than a quarter of Greece's population in the 1928 census (Katsapis, 2011:126-29) and were viewed by many as an economic burden. The already ailing economy and the politically unstable system at the time resulted in the marginalization of these refugees (Gizeli, 1984), negative attitudes toward them (Kalafati, 2018), and stigmatization of their language.

Data, speakers, and annotation
We depart from the controlled laboratory methodology followed in most Autosegmental-Metrical studies of intonation (Ladd, 2008) and instead draw on natural speech corpora, because many of the sociolinguistic parameters affecting the behavior of bilingual speakers are not well understood and cannot be replicated in the laboratory. We extracted 2,977 continuation rise utterances from 111 speakers (seventy-one male, forty female) varying in length, complexity, lexical makeup, syntactic structure, and style from the sources listed in Appendix A. Table 1 presents the number of tokens and speakers analyzed per language variety and generation.
AMG speakers here were divided into four generations based on date of birth and historical circumstances. For consistency with the generational categorization in one of the AMG corpora we used (Janse, 2015), we classified speakers born in Turkey before 1923 and their Turkish and Athenian contemporaries as Generation 1. Subsequent generations were divided as follows: second-generation speakers born between 1923 and 1945 (the oldest was in fact born in 1927), third-generation from 1950 to 1975, and fourth-generation after 1978. 2 The oldest speaker in the recordings across all three varieties was born in 1879, the youngest in 1994, their birthdates spanning over a century.
The first-generation AMG speech was drawn from archival as well as contemporary field recordings, thus some members of this generation were recorded in their youth, others in their old age. For this generation, the shape of the modeled f 0 curves and the location of peaks and troughs in continuation rise tunes in archival recordings were compared to those in contemporary recordings (see Statistical analysis section and Baltazani et al., [2020] for details of this method). This comparison revealed no significant differences in f 0 shape or in the location of peaks and troughs due to recording date (Baltazani et al., 2020), so all the data on first-generation AMG speakers regardless of the recording date are pooled together in the subsequent analysis.
We defined continuation rises in all three varieties as phrases marked with an H tone on their right boundary (in Autosegmental-Metrical terms, H− or H%), indicating that the speaker has not finished speaking (see Analysis of f 0 contours section for details). In the Athenian continuation rise tune there is typically an L* nuclear pitch accent (i.e., an f 0 trough) aligned with the stressed vowel of the nuclear word, followed by an H− phrase-final accent (i.e., the phrase ends in high f 0 ; Arvaniti & Baltazani, 2005;Baltazani & Jun, 1999;Figure 1 top). In the Turkish continuation rise tune a H*+L nuclear pitch accent (i.e., an f 0 peak plus a fall) is followed by a H− phrase accent (Ipek & Jun, 2014;Özge & Bozsahin, 2010; Figure 1 bottom). The f 0 movement is a simple rise in Athenian but a rise-fall-rise in Turkish.

Planned comparisons
For all comparisons, we examined the modeled f 0 curves to find differences in their shape, including location of peaks and troughs (see Analysis of f 0 contours section for details). Two sets of comparisons were conducted: (1) Geographical region comparison (for Turkish only). The Turkish recordings were from several regions, including Istanbul, the Black Sea, the Aegean, and Central Anatolia. (Recordings from the Doegen corpus of Turkish speakers from the territory of the Ottoman Empire in what is now Romania are included as a fifth sample in order to examine whether or not this data is similar to the other varieties); and (2) Diachronic comparison. This was a comparison across the three varieties and four generations to determine the similarities of AMG to Turkish and Athenian, and whether or how their intonation may have changed over the generations.  Analysis of f 0 contours Our quantitative analysis of the shape and relative timing of f 0 contours of continuation rises involves the following steps: 1. Standardization of audio file formats. 2. Identification of utterances containing continuation rises.
3. Estimation of f 0 contours, quality control, and elimination of f 0 tracking errors. This step includes transformation of f 0 to semitones above or below the utterance mean. 4. Functional Data Analysis I: Modeling the f 0 contours of whole utterances, using tenth order polynomials, in order to precisely determine the timing of the trough that is used to compare alignment between the three languages. Alignment parameter τ is estimated in this step. 5. Functional Data Analysis II: Modeling the f 0 contour of the continuation rise portions, using fourth-order polynomials. Shape parameters of the continuation rises are estimated in this step.
For these steps we use a variety of standard signal processing tools: sox, ffmpeg, ESPS, and GNU Octave (Eaton et al., 2019), a public-domain package that is similar to Matlab in its syntax and functionality. Subsequent statistical analysis of the modeled contour features was carried out in R.
Standardization of audio file formats. All our data sources were digital, in a variety of formats (e.g., mp3, mp4, and .wav PCM; 2-channel or monophonic), bit rates or sampling rates (e.g., 44.1 kHz, 22.05 kHz, or 16 kHz). Additionally, some digital recordings made from ¼-inch tape, recorded at different tape speeds, required speeding up or slowing down by a factor of 2 or ½ to restore the correct original recording rate. A small number of such digitized tape recordings ran backwards on one channel as the tape spool had originally been turned over for the second half of a monophonic recording, but it had been digitized as if it were a two-track stereo recording. To permit for the subsequent functional data analysis steps to be performed as batch computations, we converted all the recordings to 16 kHz, monophonic, uncompressed PCM .wav audio files. 3

Identification of utterances containing continuation rises.
With the assistance of native speaker researchers, we employed pragmatic criteria for choosing the utterances we would analyze in all three language varieties to avoid the circularity that would ensue had we based our selection solely on intonational features. The included phrases were parts of broad focus declaratives in which no constituent carries narrow focus, which were nonfinal in the speaker's turn, and ended in a high boundary that typically (but not always) was followed by a short pause. In some cases, this high boundary tone was followed by a slight fall in pitch, but we took the closeness of the high tone to the end of the phrase as the only intonational criterion; the internal intonational details of the phrase were not used for the selection.
For each variety, the native speaker researchers identified the relevant utterances, orthographically transcribed them, translated them into English, located the nuclear word and manually annotated the beginning and the end of the stressed vowel ("v") using Praat (Boersma & Weenink, 2020) as in Figure 3 (top). The start of this stressed vowel served as the beginning of the "Region of Interest," that is, the part of the continuation rise tune compared across all three languages extending from the start of the vowel to the end of the utterance.
Estimation and verification of f 0 contours. For each utterance, f 0 was measured every 10 ms (1 cs) using the ESPS (Entropics Signal Processing System) get_f0 function (Talkin, 1995), to obtain an f 0 time series.
To identify f 0 tracking errors, especially octave doubling and halving, and to provide a degree of normalization of between-speaker differences in voice pitch, the f 0 contours of each utterance were converted to semitones above or below the whole utterance mean, f 0 , as in (1). This makes every token directly comparable to all the others and goes some way to normalizing for interspeaker differences. If the data had included better information about speaker identities, we might preferably have transformed f 0 to semitones above or below the speaker mean, but we lack that data.
Given this transformation, an octave error is a sudden jump of twelve semitones (or thereabouts; it can be twelve semitones + some small change in f 0 ), as shown in Figure 3. All such potential octave jumps were individually inspected to check that they were octave errors, not actually occurring large pitch changes; the error portions were then corrected by doubling or halving (as appropriate) the f 0 measurements throughout the mistracked interval. Other large or sudden jumps in f 0 were similarly inspected to determine their cause. Rapid pitch perturbations due to consonantal onsets or offsets were retained unaltered as they are a natural part of f 0 contours. In other cases, rapid jumps could be attributed to environmental noises or to other speakers; such portions were zeroed out in the measured f 0 tracks.
Modeling the f 0 contours of whole utterances.
To determine differences or changes in the shape of continuation rises, we employ techniques of Functional Data Analysis (Ramsay, Hooker, & Graves, 2009), a statistical approach to analyzing continuous data such as curves, signals, or surfaces used in a broad spectrum of disciplines such as physiology (growth curves), demographics (population variables), weather forecasting, and, more recently, speech (Andruski & Costello, 2004;Aston et al., 2010;Chen & Wang, 1990;Grabe et al., 2005Grabe et al., , 2007Gubian, Cangemi, & Boves, 2011;de Ruiter, 2011;Zellers et al., 2010). Functional data analysis begins with data smoothing and fitting to a selected "basis function," a process that transforms a time series of discrete, raw data points (in this case measured f 0 values) into a smoothly varying function. This emphasizes patterns in the data by minimizing short-term deviation due to measurement errors or inherent system noise. Functional data also allows information on the shapes of the curves to be obtained, which facilitates comparisons across languages, language varieties, genders, generations of speakers, and other sociolinguistic or phonetic variables. A significant advantage of functional data analysis is that it augments the highly abstract and impressionistic Autosegmental-Metrical analysis of intonation with quantitative, empirically testable models of tunes, allowing comparisons of large numbers of pitch curves of whole utterances or their parts and facilitating the study of melodic variability, co-occurrence of patterns of tone combinations in melodies, and their frequency of occurrence.
Tenth-order polynomials (2) were then fitted to the f 0 contours using the GNU Octave/Matlab polyfit function (Eaton et al., 2019).
Tenth-order polynomials were used to fit the f 0 contour of utterances that may be long, with up to nine peaks and troughs (e.g., Figure 2 top). That was sufficient for this dataset; in longer utterances with more than nine peaks and troughs, an even higher-order polynomial could be used. 4 Figure 2. An example of a tenth-order polynomial fitted to the f 0 contour of an utterance. Top: the smooth modeled curve is superimposed on the observed curve, which is characterized by rapid pitch perturbations, some due to consonantal onsets or offsets. Note the two voiceless stretches in the observed data at around 50 and just before 100 cs. Bottom: The modeled f 0 maxima and minima locations are indicated by stars at the relevant times.
Fitting a specific polynomial function to the data also permitted us to find the peaks and troughs in the f 0 contour analytically, by differentiating (2), since they lie at points where its first derivative (3) is zero. 5 df 0 /dt = 10a 10 t 9 + 9a 9 t 8 + . . . + a 1 (3) The roots (values of t at which df 0 /dt = 0) were found using the GNU Octave/Matlab function real(roots(polyder(a))) . These roots provide the times of all the maxima and minima (peaks and troughs) of the intonation contour (Figure 2 bottom). Putting these times back into equation (2) gives the f 0 of those peaks and troughs.
Since there are nine such turning points in a tenth-order polynomial, we manually selected the first local minimum after the start of the nuclear vowel in each continuation rise for the alignment comparison. As this local minimum is present in all the three varieties, this choice ensures crosslinguistic comparability. We determined the alignment (τ) of this local minimum to a segmental landmark, the end of the nuclear vowel; specifically, τ = t(min f 0 )t(V end). τ > 0 means that the trough follows the end of the nuclear vowel (Turkish pattern) and τ < 0 means that the trough precedes the end of the nuclear vowel (Athenian pattern).
Thus, polynomial modeling enables us to calculate very precisely the location of peaks and troughs in the f 0 track. Prior work on intonation has often found it difficult to estimate the location of such turning points due to the presence of voiceless stretches or microprosodic perturbations (e.g., the voiceless stretches in Figure 2), finding recourse instead in specially designed laboratory utterances containing mostly sonorant segments. Simply picking maxima and minima of observed f 0 is a rather questionable and inexact method (Kochanski, 2010). The location of a high or low f 0 "target" may coincide with a portion of voicelessness in the speech signal, and the observed f 0 contour will not have a measured value at such a point. In contrast, the method used here is able to infer an estimate of the peak or trough location and its modeled f 0 value even during intervals of voicelessness.
Modeling f 0 in the Region of Interest/continuation rise portions. The Region of Interest, as defined in Identification of utterances containing continuation rises section, contains from one to three syllables, depending on the variety and stress position (antepenultimate, penultimate, or final). To characterize similarities or differences in the overall shapes of the continuation rises in the Region of Interest, we visually inspected the shape of f 0 in this region. The most complex shape involved two peaks and a trough, and therefore the f 0 in this region was modeled in further detail using fourth-order polynomials, an order of polynomial that is necessary and sufficient for a model that has three extrema.
We then modeled the normalized f 0 data using polynomial basis functions, following Grabe et al. (2007). This process converts discrete data points ( f 0 values) into a smoothly varying, continuous function, that is, for the fourth-order case (Figure 4), one among various possibilities. The values of constants a 0 to a 4 were found using the GNU Octave/Matlab polyfit function. Disregarding the error term ε, the modeled f 0 function in (4) was then converted into the corresponding form (5), the best-fitting sum of Legendre polynomials.
This transformation has the benefit that the c n coefficients of each term are orthogonal, unlike the a n coefficients of an ordinary polynomial such as (4). Low-ranking terms of a polynomial pick out slowly varying properties, and higherranking terms pick out successively more rapidly varying properties. Legendre polynomial L 0 models the overall level or height (roughly, the average) of the delimited portion of the data, L 1 models the slope of that extract, L 2 fits a parabola to the data, L 3 models the data as an up-down-up (or down-up-down) wavy shape, and L 4 as a more complex wavy m-or w-shape. If the sign of the respective c n coefficient is inverted, those components of the pattern are flipped upside down about the horizontal axis (see Grabe et al. [2007] for further details).
We refer to the coefficients of the lowest four components as AVERAGE (i.e., c 0 ), SLOPE (c 1 ), PARABOLA (c 2 ) and WAVE (c 3 ) (see Grabe et al. [2007]), and coefficient c 4 as M-OR-W. The orthogonality of these components, coupled with the simplicity of their shapes (unlike, for example, Principal Components), makes the physical interpretation and particular contribution to the overall f 0 contour quite straightforward and intuitive. Other basis functions are of course possible, but we consider the independent prosodic transparency of these components to be an appealing feature of this methodology.
Legendre polynomials take values in the range [1, −1], so corrected f 0 [semitones] was scaled so that the whole f 0 contour lies within the interval [1, −1], with 0 as the utterance mean and the greatest positive or negative deviation from the mean being 1 or −1, as the case may be: As a consequence of this step, the normalized average f 0 of every utterance = 0; f 0 [normalized] > 0 means "higher than average f 0 ", f 0 [normalized] < 0 means "lower than average f 0 " (see Figure 4, lower panel, for an example). The c n coefficients are calculated from the a coefficients through a set of simple formulae, the details of which vary according to the order of the polynomial. 6 The normalization or scaling needed for the transformation to Legendre polynomials reduces interspeaker variation, to some extent. In particular, it reduces individual differences in mean f 0 , f 0 range, and speech rate. But we are not interested in those parameters in this study; the normalization and scaling do not affect interspeaker differences in the shape of the intonation contours, represented by the c n coefficients, and we examine the relative timing of the pitch peaks and troughs using the a n coefficients of the standard polynomial (4), prior to scaling. Figure 5 presents orthogonal polynomial analysis of prototypical examples of Athenian and Turkish continuation rise contours, showing each of the five polynomial terms independently. The same scale is used for both varieties in order to show their differences in magnitude. The Athenian example, at left, is a typical L* H− contour with a slight initial descent, a low plateau, and a large final rise. The magnitudes of the M-OR-W and WAVE components are therefore almost 0, and those components are almost flat. PARABOLA and SLOPE are much larger; the size of the PARABOLA coefficient c 2 gives a broad and rather shallow bowl shape. The strong overall rising profile of this intonation contour is mostly determined by the positive SLOPE component.
The Turkish example (Fig. 5 right) is quite different, a typical rise-fall-rise, with a slight deceleration/flattening off at the very end. Consequently, the WAVE component is large, more than ten times the magnitude of the WAVE component in the Athenian example. It is positive, meaning that it is a rise-fall-rise. The PARABOLA component is the largest and is over eight times larger than in the Athenian example, giving this component a narrower trough in Turkish than in Athenian. The M-OR-W component is negative, giving it two peaks rather than two troughs (M-like rather than W-like); the second of those peaks captures the slight flattening-off of the final rise. The AVERAGE component in the Turkish pattern is higher than in the Athenian pattern, because Turkish has two large peaks whereas Athenian has a long low plateau that is mostly a few semitones below the average.
To appreciate the individual contributions of the five components further, appendix B provides further illustrations of the same two examples in which each component is set to zero, one by one.

Hypotheses
Four hypotheses were tested. H1 (geographic region): based on impressionistic analysis, no differences are expected in the intonation patterns of the continuation rise tune between the regional varieties of Turkish. H2 (diachronic development): based on Baltazani et al. (2020), the intonation of AMG speakers is expected to display a mixture of Athenian-like and Turkish-like patterns. In addition, H3: firstgeneration AMG speakers are expected to make more use of Turkish-like patterns and less use of Athenian-like patterns than the subsequent generations. H4:

Language Variation and Change
similarities between AMG and Turkish are expected to diminish with each subsequent generation.
Based on the prior literature on continuation rises in Athenian and Turkish, we expect any influence of Turkish on AMG to be manifested in the shape of the modeled curves, through similarity in the polynomial coefficients. For example, because Turkish displays a rise-fall-rise f 0 movement while Athenian displays a simple rise from a trough, we expect Turkish-like continuation rises to have a positive WAVE coefficient that models the data as an up-down-up wavy shape and Athenian-like ones to have a WAVE coefficient around zero, that is, no such wavy shape, though we would expect Athenian-like continuation rises to have a positive SLOPE. In addition, we expect the alignment of the trough in the Region of Interest to be relevant: the nuclear accent in Athenian is typically realized as a trough aligned within the nuclear vowel (τ < 0), while the H*+L nuclear pitch accent in Turkish is realized with an H* tone within the nuclear vowel followed by the trough, a trailing L tone, which typically occurs after the end of the nuclear vowel (τ > 0).

Statistical analysis
Because there were relatively few archival recordings of first generation AMG speakers, we pooled data from archival recordings made in the 1920s and 1930s and contemporary ones made in the 2000s, having first checked that there were no significant differences between the two types. Baltazani et al. (2020) tested the hypothesis that the influence of Turkish would be found to be greater in the archival than in the contemporary recordings. Turkish influence was expected to have diminished in contemporary recordings due to longer contact with Greek (eighty years between 1930 and 2011) and absence of contact with Turkish. A Mann-Whitney U test revealed no significant difference in the curve shape (e.g., for c 3 , U = 13,452, p = .518; for c 2 , U = 12,656, p = .129) or in the alignment (U = 13,317, p = .427) between archival recordings (mean ranks c 3 = 173, c 2 = 166, alignment = 184) and contemporary recordings (mean ranks c 3 = 181, c 2 = 184, alignment = 175).
For the geographic region comparison, six Kruskal-Wallis one-way analyses of variance (Kruskal & Wallis, 1952) tested for differences between the Turkish regional varieties in the five polynomial coefficients (c 4 to c 0 ) and trough alignment (τ).
For the diachronic development hypothesis, we used a Gaussian mixture model (e.g., Marin et al., 2005), which formalized the assumption that the distributions of shape coefficients in the AMG data are either Athenian-like (with a probability of λ) or Turkish-like (with a probability of 1−λ). As ANOVAs and frequentist linear mixed effects models are more widely used in phonetics research, the choice of a (Bayesian) Gaussian mixture model may warrant justification. The main reason is that we were not only interested in whether the coefficients varied by generation or dialect (as an ANOVA would allow us to detect). To test hypothesis 4, we also wanted to determine the relative proportions of Turkish-like and Athenian-like utterances in AMG in each generation, since we had previously found that our first generation of AMG speakers used a mixture of Athenian-like and Turkish-like patterns (Baltazani et al., 2020).
In consequence, the distributions of coefficients for every AMG speaker were assumed to be a mixture of an Athenian-like component and a Turkish-like component, scaled by the mixture proportion λ, as illustrated in Figure 6 for different examples of mixture distributions with different mixture proportions λ (left, center, and right panels) and different distances between the modes of the two components that constitute the mixture (upper and lower panels). Figure 6 also illustrates that not all mixture distributions are evidently bimodal: when the two mixture components (dotted and dashed lines, respectively) are sufficiently far away from each other relative to their variance, the resulting mixture distribution does have two clearly visible modes regardless of the mixture proportion parameter λ (panels A-C). In contrast, when the mixture components are close, with significant overlap between them, there might be only one local maximum in the mixture distribution (panels D and F). The presence of distinct modes depends on the mixture proportion λ and the distance between the distributions.
While the distributions of some coefficients for the AMG data were bimodal (as shown in Figure 9 below), similar to those in panels A-C in Figure 6, some had a single local maximum. Even in such cases, mixture proportions can be estimated if the distributions of the mixture components are known. Thus, we modeled the distribution of AMG coefficients as a mixture of Athenian-like and Turkish-like distributions. To this end, we assumed that the coefficients for Athenian and Turkish tunes were normally distributed with means μ A , μ T and standard deviations σ A , σ T , while AMG followed a mixture of the distributions 7 N(μ A , σ A ) and N(μ T , σ T ) with different mixture proportions λ g for each generation g.
We further assumed that λ g represents the proportion of Athenian-like utterances in a given generation rather than the proportion of Athenian-like speakers. This is because most AMG speakers produced continuation rises with both Athenian-like and Turkish-like patterns, as illustrated in the two examples in Figure 7, which were produced by the same speaker.
A more detailed description of the model structure is presented in appendix C. 8 All free model parameters (μ A , μ T , σ A , σ T , λ 1 , λ 2 , λ 3 , λ 4 ) were estimated from the data using a Bayesian model implemented in brms (Bürkner, 2017) and rstan (Stan Development Team, 2022a, 2022b in R (R Core Team, 2021) with by-speaker random effects for the means of the speaker-specific Athenian-like and the Turkish-like coefficient distributions μ A and μ T, accounting for the fact that each speaker's scores could be systematically lower or higher than the group averages.
To establish the validity of the Gaussian mixture model, we compared it to a baseline model for each coefficient (c 0 to c 4 ), and τ, the trough alignment. As a baseline, we used Bayesian linear mixed-effects models (e.g., Gelman & Hill, 2006;Vasishth & Nicenboim, 2016) with generation, dialect, and their interaction as fixed effects, and by-speaker random intercepts. A detailed description of the model structure is presented in appendix D. We compared the models based on their expected log pointwise predictive density (elpd), calculated using PSIS-LOO CV in the R loo package (an implementation of Vehtari, Gelman, & Gabry, 2017). The elpd statistic provides an estimate of a model's out-of-sample performance and thus implicitly penalizes excessive model flexibility, thus putting models of different complexity on an equal footing.
Each group's mean and median values are also important for the interpretation of these results. As mentioned in the section Modeling f 0 in the Region of Interest, the sign of each coefficient is important because the patterns of coefficients with negative values are flipped upside down. The median values for each coefficient in the Turkish data (Fig. 8) have the same sign across the varieties, indicating no difference in shape for the parameters, despite the significant differences in their magnitudes. In view of this, in combination with the fact that all these tokens have the same meaning of continuation rise, we interpret the results of statistical significance as phonetic microvariation. That is, the same phonological pattern for the continuation rise tune is displayed for every coefficient in all varieties (corresponding to the H*+L H− tune), though there are fine phonetic differences, for example, in the breadth of the curvature or the steepness of the slope.
For the subsequent analysis, data from all five varieties were grouped together because no phonological differences were found between them.   Figure 9 shows the distribution and variation of the coefficients and trough alignments for each generation in AMG, Athenian, and Turkish (in each panel, upper, middle, and lower rows, respectively). Vertical dashed lines superimposed on the histograms correspond to the medians of the observations for Athenian and Turkish. The histograms for AMG are overlaid with both the Athenian and Turkish median lines to indicate the approximate position of the expected modes if the mixture were bimodal.

Diachronic comparison
In Athenian and Turkish, the coefficient medians, for all generations, confirm the impressionistic descriptions of the tunes in 2.1. The positive SLOPE coefficient c 1 for Athenian indicates a rise, but the 0 for Turkish is the average of a rise-fall-rise, the shape of which is captured by other parameters. The positive PARABOLA coefficient c 2 for Turkish models the data as an accelerating f 0 movement, and the 0 for Athenian indicates that this parameter does not contribute much to the contour, mostly a low plateau. The WAVE c 3 is positive for Turkish, modeling the data as a rise-fall-rise f 0 movement, but it is close to 0 for Athenian, indicating that it is not so wavy: it is a plateau plus a rise. The median of the M-OR-W coefficient c 4 is negative in Turkish, reflecting the fact that it has two peaks: an H accent on the stressed vowel followed by a trough and then another peak at the end of the utterance. The median M-OR-W is close to 0 in Athenian as it does not have two peaks. The median alignment parameter τ is negative for Athenian, showing that the L tone is aligned before the end of the stressed vowel (approximately twenty cs), and is positive for Turkish as the L tone is aligned approximately twenty cs after the end of the stressed vowel (Fig. 9).
The AMG histogram reveals a shift from greater use of Turkish-like patterns to greater use of Athenian-like patterns over the generations. The distribution of the SLOPE c 1 coefficient (top row in Fig. 9), has a mode in between the Turkish and Athenian medians. Figure 10 shows that the median SLOPE for AMG changes over the four generations, becoming progressively more Athenian-like. The next two rows in Figure 9 show that the PARABOLA c 2 and WAVE c 3 coefficients have a bimodal appearance in the first-generation. One peak of the AMG distribution is near to the Athenian median and the other is near to the Turkish median. Over the generations, the Turkish-like peak becomes smaller, so that these distributions have an increasingly unimodal appearance as the proportion of Turkish-like values seems to decrease, a pattern like the one modeled in Figure 6F. The M-OR-W coefficient c 4 also changes across the generations: the AMG mode always coincides with the Athenian mode, but the AMG distribution is substantially wider than the Athenian with a tail toward Turkish values, a pattern as in Figure 6D. The Turkish-like aspect of this tail in the first generation is reduced in subsequent generations.
The bottom row of Figure 9 gives the visually most striking finding: the distribution of the L tone alignment τ is clearly bimodal in the first three generations, the two peaks in the first generation almost coinciding with the medians of Athenian and Turkish, as in Figure 6B. For this parameter as well, the proportion of Turkish-like values decreases over the generations, giving rise to a more Athenian-like distribution, similar to Figure 6F.
In view of the above, we compared a Bayesian analysis of the Gaussian mixture model with a baseline linear mixed effects model for all five dependent variables as described in 2.6. The results of the model comparison using PSIS-LOO CV are  presented in Table 2. They show that the mixture model provided a better description of the data than the linear mixed effects model for all parameters. The mixture model only weakly outperformed the baseline model for c 1 but showed a substantial improvement over the baseline for the rest of the parameters.
In the Gaussian mixture model, estimates of the mixture parameter λ (proportion of Athenian-like values, expressed as a percentage) for the remaining parameters are shown in Table 3 and Figure 10. Several trends are revealed for all parameters except SLOPE c 1 . First, the median estimates of λ are very close to Turkish values (less than 7% Athenian-like) in generation 1. Second, in generation 2, the median estimates of λ are gradually approaching Athenian (11% -21% Athenian-like). Third, the median estimates of λ in generation 3 are over 68% Athenian-like. Among all the parameters we examined, only c 1 showed higher Athenian-like values than Turkish-like values. Yet, the rate of this similarity increased gradually from the first-generation to thefourth generation as well.
All in all, these results reveal a consistent pattern of a change over the generations in the realization of continuation rise tunes. The strong similarity to Turkish in AMG generation 1 weakens over the generations, with Athenian characteristics becoming predominant in generations 3 and 4. Table 2. Results of the model comparisons between (i) the baseline linear mixed effects model and (ii) the Gaussian mixture model for the coefficients c 1 −c 4 , and τ using PSIS-LOO CV estimates of elpd (Vehtari et al., 2017). Larger elpd LOO values indicate more parsimonious models. Positive Δelpd LOO values indicate that the Gaussian mixture model provides a more parsimonious account than the baseline linear mixed-effects model

Discussion
In this paper, we examined how long contact effects in intonation persist after contact has ceased. We compared the f 0 patterns in continuation rise tunes of Athenian Greek, Turkish, and Asia Minor Greek (AMG) across four generations. We modeled the sampled f 0 data using Legendre polynomial basis functions to characterize similarities or differences between these varieties in the final stretch of the continuation rises, comprising the nuclear pitch accent and the edge tones. This modeling allowed the analysis of a fairly large dataset of approximately three thousand continuation rise utterances from diverse corpora containing mostly natural speech in conversations and narratives. Using fitted polynomials for the analysis and comparison of the tunes in the three language varieties captures simple properties of the fitted curves and has an easily understood, cognitively plausible interpretation, such as the average f 0 , whether the contour is generally rising or falling, the size and location of intonation peaks and troughs, as well as the number of peaks and troughs in the intonation contour. Furthermore, the functional data analysis can be related to the well-established Autosegmental-Metrical analysis of intonation: the derivative of the quartic polynomial used here has three roots, corresponding to the time-points of the two peaks and the trough, which allow us to estimate their corresponding f 0 more reliably than through the observed peaks and troughs, which are susceptible to noise and mistracking, especially if they occur during voiceless intervals.
The results of the diachronic comparison across four generations of Turkish and Athenian speakers confirmed the realization of the continuation rise tunes reported in the literature for these two languages. The Turkish continuation rise tune is reported as a rise-fall-rise H*+L H−, the H* tone in the H*+L nuclear pitch accent occurring within the nuclear vowel. In our data the median alignment τ of the trailing L tone is 200 ms after the nuclear vowel offset, and the H− phrase accent lies at the right edge of the phrase. The negative median M-OR-W c 4 coefficient in our modeling reflects the fact that the Turkish tune has two peaks, the positive WAVE c 3 captures its rise-fall-rise f 0 movement, and the positive PARABOLA c 2 captures a general rising f 0 movement. The Athenian tune also ends in an H− phrase accent but has a L* nuclear pitch accent that is aligned on average 200 ms before the end of the nuclear vowel. The positive SLOPE coefficient c 1 for Athenian indicates that, in general, it rises while the low plateau from the nuclear vowel to the phrase edge is reflected in the PARABOLA c 2 , WAVE c 3 , and M-OR-W c 4 coefficients, all of which are close to 0, indicating that the Athenian contour lacks a wavy shape, and that it does not have two peaks.
As for AMG, a previous analysis of this tune in two generations of speakers showed a dual pattern of the L tone alignment, one Athenian-and one Turkish-like. Here we established a more complete picture by showing that the mixed pattern of the data is true for all five parameters through a Gaussian mixture model. The bimodality of the L tone alignment found in the speech of the firstgeneration AMG speakers decreased in each subsequent generation toward a more Athenian-like alignment. A similarly diminishing proportion of Turkish-like patterns was revealed for the polynomial coefficients c 1 to c 4 . While in generation 1 the mixture coefficient λ is less than 8% Athenian-like for four of the five parameters, in generation 2 there is a greater proportion of Athenian-like coefficient values, and in the two youngest generations the greatest proportion (almost 100%) is Athenian-like.
The Bayesian approach gave us the flexibility to define a mixture model that expresses our hypotheses while allowing us to account for speaker-specific deviations from the average in the form of by-speaker random effects (e.g., Gelman & Hill, 2006;Vasishth, Nicenboim, Beckman, Li, & Kong, 2018). This approach enabled us to encode the assumption that the distribution of the AMG coefficients is a mixture of an Athenian-like and a Turkish-like Gaussian distribution, and, at the same time, estimate the mixture parameters and the uncertainty associated with them in datasets with unequal numbers of observations per speaker. For AMG speakers, all coefficients were better described by a mixture of the respective Turkish-like and Athenian-like distributions than by simple linear mixed effects models.
Contact between Greek and Turkish resulted in phonological transfer of the pitch accent found in Turkish continuation rises into AMG. The AMG speakers alternated between the Turkish and the Greek tunes in a way similar to what is reported in Queen (2012), albeit in different proportions depending on the generation. No phonetic changes such as the delayed peaks reported in O'Rourke (2012) and Elordieta and Calleja (2005) were revealed in the tunes we examined. That is, when the tune was produced with the Athenian L pitch accent, the alignment of the trough was Athenian-like, that is, within the nuclear vowel. Similar results were found for the Turkish-like tune. It is not yet clear whether these two patterns of continuation rise differ in pragmatic interpretation or whether they are codeswitched equivalents. We leave such questions for future investigation.
This paper demonstrates intonational effects of language contact, a generally understudied area, since so far the knowledge about the results of contact has mainly relied on morphosyntactic or segmental phenomena. We have shown that prosodic characteristics of Turkish have persisted in AMG for close to a century after the cessation of contact, which we attribute to the strong sense of ethnolinguistic identity in the AMG-speaking community. However, usage of the Turkish-like intonation pattern has diminished over the generations due to a complex mixture of factors such as the prestige of Athenian Greek, the stigma that some attach to AMG, and the fact that the younger generations increasingly have less exposure to their AMG heritage variety, at least the Cappadocian variety examined here, than earlier generations did.
Competing interests. The authors declare none.

2.
For the majority of the sources, some metadata with sociolinguistic information were available (either directly or indirectly from public information available on the web) that included at least the speaker's age. Where information about the speaker's age was not available (marked with a star next to the token number in Appendix A) approximate age and classification to one of the four generations was inferred from characteristics of the speaker's appearance in videos or information about their career. 3. Although mp3 and mp4 files typically use audio compression with some loss of fidelity, Fuchs and Maxwell (2016:525) found that mp3 "compression rates between 56 and 320 kbps show relatively small mean errors of 2% or less, with median errors well below 0.5%. … mp3 compressed data is viable for the analysis of F0." They found that measurement errors in mp3-compressed data were greatest of all in consonants, for example, the perturbations to f 0 caused by obstruents, not relevant here. The variation in measured f 0 due to different estimation algorithms and analysis parameters can be far greater than the 0.5% to 2% error rate arising from mp3 compression. 4. Grabe, Kochanski, and Coleman (2005) also used polynomial equations to describe f 0 in complete intonation phrases, to investigate and compare global properties of statements and questions. The utterances in that study were much shorter, simpler, lab-elicited sentences, for which third-order polynomials were sufficient. Sherr-Ziarko (2019) used 25th-order polynomials (!) to first fit and then segment f 0 contours in long Japanese utterances containing many pitch accents. 5. Not all choices of basis functions support this kind of analytical approach to finding minima and maxima; it is an attractive feature of polynomial modeling that differentiation and root-finding is so easy. 6. For fourth-order polynomials, the transformations are: c 4 = 8/35 a 4 ; c 3 = 2/5 a 3 ; c 2 = ⅔ (a 2 -30/35 a 4 ); c 1 = a 1 + 3/5 a 3 ); c 0 = a 0 -3/35 a 4 + ⅔ (a 2 -30/35 a 4 ). See Grabe et al. (2007) appendix C for the transformations required to orthogonalize third-order polynomials. An explanation of the form of such transformations for polynomials of an arbitrary order requires quite a bit of algebra and lies well beyond the scope of this paper and this journal. 7. The notation N(μ, σ) stands for a Gaussian distribution with mean μ and standard deviation σ. 8. All data and code used in this analysis are available at https://osf.io/mus6v/. 9. We obtained many of these films from YouTube, but as they have since been deleted, we cannot provide their URLs here. Appendix B. Independent contributions of the five orthogonal polynomial components In the section Modeling f 0 in the Region of Interest we presented an orthogonal polynomial analysis of prototypical examples of Athenian and Turkish continuation rise contours. Figure 5 showed the contribution of each of the five polynomial terms independently. In order to appreciate the individual contributions of the five components further, Figure B1 shows how the overall polynomial fit in the same two examples is affected when each component is set to zero, one by one. In the Athenian example (Fig. B1, left) setting c 4 , c 3 , or c 0 to 0 has very little effect on the overall fit (those panels all have a low plateau with a final rise), because those coefficients were close to 0 in any case. Setting c 2 to 0 results in the fit losing its very broad PARABOLA component, so that the fitted contour is now rising rather than gently falling at the start. Since the SLOPE (c 1 ) component is doing a lot of the heavy lifting in this analysis of the Athenian contour, setting c 1 to 0 results in the fitted curve mostly losing its upward-rising profile at the end.
In the Turkish example (Fig. B1, right), setting the M-OR-W component c 4 to 0 eliminates the two distinct peaks seen in Figure 5. Setting c 3 to 0 eliminates the initial rise, as this WAVE component contributes strongly to the rise-fall-rise contour of the fitted f 0 pattern. In both of these cases, what remains is a falling-rising bowl-like pattern which does not model very well the much more dynamic rise-fall-rise (with final tail-off) of the Turkish pattern. Setting c 2 to 0 results in a fitted curve that does show those two peaks, but now the trough between them is no longer so deep, because of the removal of the PARABOLA component. Setting c 1 to 0 has very little effect on the overall fit as it was close to 0 (−0.07) anyway. Setting c 0 to 0 has no effect at all on the overall shape of the fitted curve, as this is just a constant term; the y-axis scales of the two lower-right panels differ by about 2.4 (c 0 = 1.6654), but the shapes are otherwise similar. The baseline linear mixed-effects models for all dependent variables included fixed effects for generation, dialect, as well as their interaction, and by-speaker intercept as stated in the brms model formula in D1. Generation and dialect were represented using successive differences contrast matrices given in Tables C2 and D1 (e.g., Venables & Ripley, 2002). Tables D2 and D3 list the prior distributions for the fixed effects and random effects of the model respectively. y 1 + (G (2−1) + G (3−2) + G (4−3) ) * (D (AMG−Ath.) + D (Tur−AMG) ) + (1|participant) (D1) a l , b λ,1 , b λ,2 , b λ,3 Normal(0, 2)