Supplementary Material

Supplementary Material

for the article Diphthong Dynamics under Lexical Tone: Cross-Dialectal Evidence for Category-Specific f0--F1 Coupling

Authors: Chenyu Li and Jalal Al-Tamimi from Laboratoire de linguistique formelle, CNRS, Paris 75013, France

S0. Full analysis pipeline and reading guide

The supplementary material follows the main stages of the analysis. S1 documents data screening and f0 trackability. S2 reports the GAMM analyses and the GAMM-based cascade reconstruction. S3 reports the primary pooled 2 × 2 z-scale fPCA analyses, including bridge and relative-configuration results. S4 provides the speaker-centered non-z robustness reruns. S5 reports the full-curve FoF analyses and residualized tone-recovery metrics. S6 reports the supplementary seven-group fPCA analyses and their non-z reruns, including the multivariate check in S6g. S7 gathers reconstruction-oriented visual summaries across methods.

Full GAMM model summaries, diagnostics, and complete compareML outputs are retained in the accompanying analysis notebook: SuppPub2_GAMM_summaries.html.

S0a. Participants, variety assessment, experimental procedure, and item inventory

Forty-eight native speakers participated: 24 from Xi'an (12 female, 12 male; 18–65 years; mean 42.875, SD 13.19) and 24 from Tianjin (12 female, 12 male; 20–65 years; mean 44.08, SD 14.10). All were native speakers of their local Mandarin variety and reported regular daily use of Standard Mandarin (SM). Local-variety background and SM proficiency were assessed using a questionnaire, self-report, and experimenter evaluation.

Each participant produced both SM and their local variety. Within each location, half completed the SM block first and half the local-variety block first. Before each block, participants read several naturalistic sentences in the intended variety to establish the corresponding speech mode. Target items were randomized within each variety block. Each item was repeated twice within the two carrier sentences described below. Misproduced items were repeated; persistently inaccurate tokens were excluded. After recording, two native SM speakers conducted a blind check of tone-category accuracy across diphthongs.

Recordings took place in a quiet sound-insulated environment using a Zoom H6 Handy Recorder and Shure WH20 head-mounted dynamic microphone, at 48 kHz and 24-bit resolution. The onset inventory was restricted to non-aspirated bilabial and alveolar consonants to limit anticipatory coarticulation and aspiration-related glottal effects. The two carrier sentences provided contrasting preceding tonal contexts within each variety and minimized systematic influences from the following syllable. Acoustic alignment and tracking procedures are detailed in S1f.

The stimulus materials were elicited in two carrier sentences: “Please read X twice” (请将X读两遍) and “Please play X eight times” (请播放X八遍) corresponding to IPA forms tɕʰjəŋT3 tɕjaŋT1 X tuT2 ljaŋT3 pjɛnT4 and tɕʰjəŋT3 poT1 faŋT4 X paT1 pjɛnT4. The item space combined two falling diphthongs (/ai/, /au/), four non-aspirated bilabial or alveolar onsets (/m/, /p/, /n/, /t/), and the four tonal categories (T1–T4), yielding 32 theoretically possible onset-by-diphthong-by-tone slots. Because several slots lacked suitable lexical items with attested characters, the realized stimulus set contained 28 items.

The realized items were:

S1. Data summary and exclusion profile

The functional analyses (fPCA analyses and function-on-function mapping) and the GAMMs were based on overlapping but not identical datasets. Because the functional analyses treated each token as a full trajectory, they required sufficiently continuous f0 curves: tokens with fewer than five observed f0 points were excluded, and any remaining missing values were filled by linear interpolation with edge filling. The GAMMs, by contrast, were fitted to broader observationwise datasets and therefore did not require the same degree of token-level f0 continuity.

Table S1a. Overall sample sizes by analysis type

Table S1b. Exclusion summary by variety

Table S1c. Exclusion breakdown by tone, variety, and diphthong

Table S1d. Highest-exclusion cells

Because functional exclusion is triggered by insufficient observable f0 support for interpolation, Table S1c also serves as the operational summary of f0 trackability loss by condition. The strongest exclusion concentrations occur in low-tone cells, especially SM T3 in both datasets.

S1e. Pointwise f0 missingness models

To complement the token-level exclusion audit above, we also modelled pointwise f0 missingness directly with binomial GAMMs in the Tianjin-speaker and Xi'an-speaker datasets. These models estimate the population-level probability that f0 is missing at each normalized time point for each ToneVarDiphInt condition. Together with Table S1c, they provide a more direct picture of f0 trackability loss.

Figure S1-1. Predicted probability of missing f0 in the Tianjin-speaker dataset

Population-level missingness curves from the Tianjin-speaker GAMM, with the strongest probabilities concentrated in low-tone cases (SM T3; TJ T1 and T3).

Table S1e-1. Tianjin-speaker missing-f0 summary by condition

Table S1e-2. Tianjin-speaker missing-f0 summary by condition and time region

The time regions were defined over the 11-point normalized trajectory as early (points 0–3), middle (4–7), and late (8–10), and the reported values are mean predicted missing-f0 probabilities within each region.

Figure S1-2. Predicted probability of missing f0 in the Xi'an-speaker dataset

Population-level missingness curves from the Xi'an-speaker GAMM. The most severe missingness again concentrates in low-tone cases (SM T3 and XA T1).

Table S1e-3. Xi'an-speaker missing-f0 summary by condition

Table S1e-4. Xi'an-speaker missing-f0 summary by condition and time region

In both datasets, low-tone conditions, especially SM T3, have the highest predicted probability of missing f0, consistent with the token-level exclusion patterns in Table S1c.

S1f. Segment boundaries, f0 extraction, and formant tracking

Segmental alignment was obtained with the Montreal Forced Aligner and manually corrected. Vowel onset was placed at the emergence of clear formant structure and offset at its disappearance. Boundaries were primarily guided by F2 and cross-checked against F1–F4. Spectrograms were generally inspected with Praat's default 70-dB dynamic range, adjusted to 50 dB when necessary for low-amplitude signals.

Acoustic measurements were extracted in Praat 6.3.02. f0 was estimated with a two-pass autocorrelation procedure: the first pass used a broad 75–600 Hz range, and the first and third quartiles of the resulting track defined an adaptive range for the second pass. The extracted f0 track was smoothed with Praat's built-in 10-Hz bandwidth smoothing function. F1 was extracted with the Burg method followed by Praat's tracking procedure. Five formants were estimated; the maximum formant frequency depended on speaker sex and diphthong as follows.

Table S1f-1. Maximum formant settings

f0 and F1 were sampled at 11 equidistant points across each diphthong. The primary representation used within-speaker z-scores. The non-z reruns used f0 in semitones relative to the speaker's mean and F1 in speaker-centered Bark. Praat returned no valid f0 estimate at some points; GAMMs retained partially observed tokens and omitted observations with missing model-required values. Functional analyses excluded tokens with fewer than five observed f0 points and filled remaining gaps by linear interpolation with edge filling. These trackability-related exclusions are quantified in S1a–S1e; the FoF check without 10-Hz smoothing is reported in S5h.

S2. GAMM analyses

S2a. Baseline and mechanistic GAMM fit summaries

The baseline GAMMs model tone-conditioned F1 and f0 trajectories over the normalized time axis. The mechanistic GAMMs replace the explicit tone term with time-varying f0, duration, and interaction structure. Categorical predictors (including ToneVarDiphInt and covariates) in these GAMMs were entered as ordered factors using treatment coding, and random smooths were included for speaker and item (character).

For full GAMM summary and gam.check results, please refer to the document "Full GAMM Summaries and gam.check results".

Baseline:    F1 ~ s(time) + s(time, by = ToneVarDiphInt) + covariates + random smooths + AR(1)
Mechanistic: F1 ~ s(time) + s(f0) + s(duration) + ti(time, f0, duration) + covariates + random smooths + AR(1)

Table S2a. Main GAMM fit indices under z and speaker-centered non-z scaling

Figure S2-1. Baseline F1 curves in the Tianjin-speaker dataset

Baseline F1 trajectory panels from the z-scale analysis. Solid lines show the predicted trajectories, and shaded ribbons show the corresponding confidence intervals.

Figure S2-2. Baseline F1 curves in the Xi'an-speaker dataset

Baseline F1 trajectory panels from the z-scale analysis. Solid lines show the predicted trajectories, and shaded ribbons show the corresponding confidence intervals.

Figure S2-3. Baseline f0 curves in the Tianjin-speaker dataset

Baseline f0 curves by tone for the Tianjin-speaker dataset. Solid lines show the predicted trajectories, and shaded ribbons show the corresponding confidence intervals.

Figure S2-4. Baseline f0 curves in the Xi'an-speaker dataset

Baseline f0 curves by tone for the Xi'an-speaker dataset. Solid lines show the predicted trajectories, and shaded ribbons show the corresponding confidence intervals.

Figure S2-5a. Pairwise F1 difference smooths in the Tianjin-speaker dataset

Combined pairwise F1 difference smooths across sex, variety, diphthong, and the six tone contrasts in the Tianjin-speaker dataset. Shaded regions indicate intervals of significant difference.

Figure S2-5b. Pairwise F1 difference smooths in the Xi'an-speaker dataset

Combined pairwise F1 difference smooths across sex, variety, diphthong, and the six tone contrasts in the Xi'an-speaker dataset. Shaded regions indicate intervals of significant difference.

Table S2b. Mechanistic tone-vs-no-tone comparison by compareML

The two models in Table S2b differ only in whether tone is allowed to structure the remaining time-varying F1 trajectory after the shared f0- and duration-based terms have been specified. In the no-tone model, diphthong-by-variety (DiphVar.ord) structures the factor-specific time smooth, whereas in the with-tone model this role is taken by the fuller tone-by-variety-by-diphthong factor (ToneVarDiphInt.ord). The f0- and duration-related smooth terms themselves were intentionally kept the same across the two models, because the goal of this comparison was to ask whether tone still explains residual F1 trajectory shape after controlling for shared nonlinear f0 and duration effects, rather than to build a tone-specific mechanistic mapping for those smooth terms.

Figure S2-6a. Mechanistic GAMM heatmap for /au/ in the Tianjin-speaker dataset

Mechanistic /au/ surface from the z-scale analysis.

Figure S2-6b. Mechanistic GAMM heatmap for /au/ in the Xi'an-speaker dataset

Mechanistic /au/ surface from the z-scale analysis.

Figure S2-6c. Mechanistic GAMM heatmap for /ai/ in the Tianjin-speaker dataset

Mechanistic /ai/ surface from the z-scale analysis.

Figure S2-6d. Mechanistic GAMM heatmap for /ai/ in the Xi'an-speaker dataset

Mechanistic /ai/ surface from the z-scale analysis.

Table S2c. Cascade summary comparison across scales

Tianjin shows lower mean RMSE and narrower mean non-overlap than Xi'an.

Note: In the GAMM cascade summaries, predicted f0 curves were passed through the mechanistic F1 model together with duration values matched at the full condition-cell level (variety × sex × diphthong × tone × left segment × carrier-type × reading-order); error was then averaged equally across the valid matched cells within each condition.

Table S2c-1. Tone-specific GAMM cascade RMSE averaged across sex and diphthong

RMSEs (z units) were averaged equally over realized onset, carrier-sentence, and reading-order combinations within each sex-by-diphthong condition, then equally across the two sexes and two diphthongs. Each value is a mean of condition-specific RMSEs, not an RMSE calculated from averaged trajectories.

Figure S2-7. Direct-versus-cascade trajectories for /ai/ in the Tianjin-speaker dataset (z scale)

Representative trajectory-level comparison between the direct F1 fit and the cascade prediction for /ai/. These panels are shown across sex and variety, with carrier type fixed at 1, left segment fixed at m, and reading order set to the panel-specific reference level (S for SM, T for TJ, X for XA). Colored lines distinguish the direct F1 fit from the f0-based cascade, and translucent ribbons in matching colors show their confidence intervals.

Figure S2-8. Direct-versus-cascade trajectories for /ai/ in the Xi'an-speaker dataset (z scale)

Representative trajectory-level comparison between the direct F1 fit and the cascade prediction for /ai/. These panels are shown across sex and variety, with carrier type fixed at 1, left segment fixed at m, and reading order set to the panel-specific reference level (S for SM, T for TJ, X for XA). Colored lines distinguish the direct F1 fit from the f0-based cascade, and translucent ribbons in matching colors show their confidence intervals.

Figure S2-9. Direct-versus-cascade trajectories for /au/ in the Tianjin-speaker dataset (z scale)

Representative trajectory-level comparison between the direct F1 fit and the cascade prediction for /au/. These panels are shown across sex and variety, with carrier type fixed at 1, left segment fixed at m, and reading order set to the panel-specific reference level (S for SM, T for TJ, X for XA). Colored lines distinguish the direct F1 fit from the f0-based cascade, and translucent ribbons in matching colors show their confidence intervals.

Figure S2-10. Direct-versus-cascade trajectories for /au/ in the Xi'an-speaker dataset (z scale)

Representative trajectory-level comparison between the direct F1 fit and the cascade prediction for /au/. These panels are shown across sex and variety, with carrier type fixed at 1, left segment fixed at m, and reading order set to the panel-specific reference level (S for SM, T for TJ, X for XA). Colored lines distinguish the direct F1 fit from the f0-based cascade, and translucent ribbons in matching colors show their confidence intervals.

Figure S2-11. Direct-versus-cascade divergence for /ai/ in the Tianjin-speaker dataset (z scale)

Representative band non-overlap diagnostic for /ai/ under the z-scale cascade comparison, using the same fixed reference setting as Figures S2-7 and S2-8. Shaded regions indicate non-overlap intervals.

Figure S2-12. Direct-versus-cascade divergence for /ai/ in the Xi'an-speaker dataset (z scale)

Representative band non-overlap diagnostic for /ai/ under the z-scale cascade comparison, using the same fixed reference setting as Figures S2-7 and S2-8. Shaded regions indicate non-overlap intervals.

Figure S2-13. Direct-versus-cascade divergence for /au/ in the Tianjin-speaker dataset (z scale)

Representative band non-overlap diagnostic for /au/ under the z-scale cascade comparison, using the same fixed reference setting as Figures S2-9 and S2-10. Shaded regions indicate non-overlap intervals.

Figure S2-14. Direct-versus-cascade divergence for /au/ in the Xi'an-speaker dataset (z scale)

Representative band non-overlap diagnostic for /au/ under the z-scale cascade comparison, using the same fixed reference setting as Figures S2-9 and S2-10. Shaded regions indicate non-overlap intervals.

S2b. Residual structure after the no-tone mechanistic GAMM

Residual-tone diagnostics quantify the structured variance that remains after fitting the no-tone f0-duration backbone.

Table S2d. Overall residual summary by dataset

Table S2e. Mixed-model residual ANOVA summary

Residual F1 variation shows tone- and variety-dependent structure, with significant tone-by-variety interactions in both datasets.

Figure S2-15. Residual curves by tone in the Xi'an-speaker dataset

Residual curves after the no-tone mechanistic GAMM.

Figure S2-16. Residual curves by tone in the Tianjin-speaker dataset

Residual curves after the no-tone mechanistic GAMM.

Figure S2-17. Residual RMSE heatmap in the Xi'an-speaker dataset

Residual RMSE by condition after the no-tone mechanistic GAMM.

Figure S2-18. Residual RMSE heatmap in the Tianjin-speaker dataset

Residual RMSE by condition after the no-tone mechanistic GAMM.

S3. Primary pooled 2 × 2 fPCA analyses under z scaling

S3a. Common functional spaces, analysis blocks, and explained variance

Trajectories from both speaker origins were pooled to estimate one f0 basis and one F1 basis for the primary /ai/ + /au/ analysis. Speaker origin (Tianjin or Xi'an) and variety status (SM or local) define four cells: Tianjin SM, Tianjin TJ, Xi'an SM, and Xi'an XA. Local status therefore denotes TJ for Tianjin speakers and XA for Xi'an speakers. The diphthong-specific reruns each re-estimate common bases across the same four cells. Origin-specific summaries below use these common bases; the separate seven-group decompositions are reported in S6.

Table S3a-1. Analysis blocks and retained sample sizes

Table S3a-2. Explained variance proportions

In the primary pooled analysis, the first two components account for 96.05% of f0 variance and 77.44% of F1 variance. PC1 primarily captures overall level and PC2 dynamic tilt; PC3 captures finer trajectory shape. Retention counts in Table S3a-1 refer to the decomposition, while the score and bridge models use the first three predictor PCs and the relative-configuration reconstruction uses PC1–PC2.

Figure S3-1. Pooled eigenfunctions (z scale)

Domain-specific eigenfunctions estimated across both origins and both variety statuses.

Figure S3-2. Pooled trajectory reconstructions (z scale)

Tone-wise trajectories reconstructed from mean PC scores in the four origin × variety cells. These are descriptive, unadjusted reconstructions; covariate-adjusted time shapes are shown in Figure S3-3.

S3b. Score-space models, duration, and condition-block comparison

The primary mixed models predict F1-PC1 or F1-PC2 from f0-PC1–PC3, standardized duration, tone × origin × variety status, diphthong, sex, reading order, onset, and carrier type, with a speaker random intercept. The duration R² comparison is a companion ordinary least-squares calculation with the same fixed predictors; its R² values are not mixed-model marginal or conditional R². A likelihood-ratio comparison was also conducted between each full mixed model and a reduced mixed model containing f0-PC1–PC3, duration, the non-tonal covariates, and the speaker random intercept, but omitting the entire tone × origin × variety-status block. Both models were fitted by maximum likelihood. Because the reduced model omits origin and variety status as well as tone and their interactions, this comparison tests the joint contribution of the complete 15-parameter condition block rather than a tone-only effect.

Table S3b-1. Incremental contribution of duration

Table S3b-2. f0-PC and duration coefficients in the primary mixed models

Duration contributes chiefly to F1-PC2 (ΔR² = 0.06498), with negligible change for F1-PC1 (ΔR² = 0.0000331). The conditional F1-PC1 coefficient for f0-PC1 is positive in the mixed model; this conditional coefficient should be distinguished from the negative unadjusted score correlations and the inverse relative tone configurations.

Table S3b-3. Unadjusted cross-domain score correlations in the common bases

Table S3b-4. Likelihood-ratio comparison of the full and reduced mixed models

For both F1-PC1 and F1-PC2, the full model fitted substantially better than the corresponding reduced model. Thus, tone, origin, variety status, and their interactions jointly accounted for F1 score variation beyond the low-dimensional f0 scores, duration, and the non-tonal covariates.

Figure S3-3. Covariate-adjusted pooled PC time shapes (z scale)

Tone-wise PC1 and PC2 contributions reconstructed from mixed-model estimated marginal mean scores across the four origin × variety cells.

S3c. Origin-specific bridge reconstruction summaries

Within the common pooled bases, separate linear bridges are fitted for each origin: each F1-PC1–PC3 score is predicted from f0-PC1–PC3. R² describes token-level in-sample fit. Absolute and Euclidean discrepancies are computed between direct and predicted tone × variety mean scores and averaged equally across the eight conditions within each origin; they are not token-level prediction errors.

Table S3c-1. Origin-specific bridge fit and condition-mean discrepancy

Xi'an has higher bridge R² for PC1 and PC2 (0.174 and 0.034) than Tianjin (0.015 and 0.007), but a larger PC1–PC2 condition-mean discrepancy (0.559 versus 0.297). Predictability and discrepancy therefore describe different aspects of recoverability.

S3d. Covariate-adjusted relative tonal configuration

Separate mixed models for each domain's PC1 and PC2 estimate tone × origin × variety marginal means, adjusting for duration and the non-tonal covariates with a speaker random intercept. Reconstructed trajectories are centered across the four tones at each time point within each origin × variety cell and domain. Spearman correlation measures agreement in tone rank; Pearson correlation measures agreement in centered spacing. The summary is the mean of the 11 pointwise correlations. Negative correlations indicate inverse configurations.

Table S3d-1. Adjusted relative-configuration summary

Both Xi'an cells show stronger inverse alignment than the corresponding Tianjin cells.

Figure S3-4. Adjusted relative f0–F1 configuration (z scale)

Centered trajectories reconstructed from adjusted PC1–PC2 scores; each row identifies one origin × variety cell.

Figure S3-5. Adjusted relative similarity over time (z; /ai/ + /au/)

Correlations across the four tones at each of the 11 time points, calculated from covariate-adjusted PC1–PC2 reconstructions centered within origin × variety cell and domain. Negative values indicate inverse ordering or spacing.

S3e. Covariate-adjusted direct-versus-bridged comparisons

For this supplementary paired-score evaluation, one briged is fitted across all four cells, then compares direct and bridged token scores using method × ToneOriginVar, onset, reading order, carrier type, and diphthong, with random intercepts for speaker and token. This pooled bridge differs from the origin-specific bridges summarized in S3c. Estimated marginal means and within-condition contrasts describe the pooled bridge's remaining condition-structured discrepancies. Each condition has one direct-minus-bridged contrast; the reported p values have no adjustment across the complete collection of conditions and PCs.

Table S3e-1. Direct-minus-bridged contrasts for all three F1 PCs

Figure S3-6. Pooled bridge marginal means (z scale)

Covariate-adjusted direct and bridged F1-PC1–PC3 scores with 95% confidence intervals. Bridged scores come from the all-origin pooled bridge described in S3e.

S3f. Split-by-diphthong pooled 2 × 2 checks

Each diphthong is analyzed in its own common f0 and F1 spaces across both origins and variety statuses. The models follow S3b–S3e, omitting diphthong as a covariate. Origin-specific bridge summaries and adjusted relative configurations are reported separately.

Table S3f-1. Incremental contribution of duration

Table S3f-2. /ai/ origin-specific bridge summary

Table S3f-3. /ai/ adjusted relative configuration

Figure S3-7. /ai/ pooled bridge marginal means (z scale)

Direct and bridged estimated marginal means with 95% confidence intervals, using the all-origin bridge within this diphthong.

Figure S3-8. /ai/ adjusted relative configuration (z scale)

Adjusted and centered PC1–PC2 reconstructions across the four cells.

Figure S3-9. Adjusted relative similarity over time (z; ai)

Correlations across the four tones at each of the 11 time points, calculated from covariate-adjusted PC1–PC2 reconstructions centered within origin × variety cell and domain. Negative values indicate inverse ordering or spacing.

Table S3f-4. /au/ origin-specific bridge summary

Table S3f-5. /au/ adjusted relative configuration

Figure S3-10. /au/ pooled bridge marginal means (z scale)

Direct and bridged estimated marginal means with 95% confidence intervals, using the all-origin bridge within this diphthong.

Figure S3-11. /au/ adjusted relative configuration (z scale)

Adjusted and centered PC1–PC2 reconstructions across the four cells.

Figure S3-12. Adjusted relative similarity over time (z; au)

Correlations across the four tones at each of the 11 time points, calculated from covariate-adjusted PC1–PC2 reconstructions centered within origin × variety cell and domain. Negative values indicate inverse ordering or spacing.

The supplementary multivariate f0 + F1 analysis belongs to the separate-group checks and is reported in S6g.

S4. Speaker-centered non-z robustness analyses

The pooled 2 × 2 fPCA analysis is repeated with f0 in semitones relative to each speaker's mean and F1 in speaker-centered Bark. The primary pooled and diphthong-specific runs retain the same 4,760, 2,174, and 2,586 tokens, respectively. Models, bridge definitions, and adjusted relative-configuration calculations follow S3. Duration remains standardized. The earlier separate-group non-z analyses are retained in S6h–S6i.

S4a. Pooled non-z score, bridge, and relative-configuration summaries

Table S4a-1. Incremental contribution of duration

Table S4a-2. Non-z origin-specific bridge summary

Table S4a-3. Non-z adjusted relative configuration

Table S4a-4. Non-z pooled direct-minus-bridged contrasts

As in S3e, the marginal-mean contrasts evaluate the all-origin pooled bridge, whereas Table S4a-2 summarizes separately fitted origin-specific bridges. Contrast p values are not adjusted across all conditions and PCs. Duration again contributes mainly to PC2 (ΔR² = 0.04640), with negligible change for PC1 (0.00000423). The origin-specific PC1–PC2 discrepancy remains smaller for Tianjin (0.293) than Xi'an (0.455), and the Xi'an relative configurations remain more strongly inverse.

S4b. Split-by-diphthong pooled non-z summaries

Table S4b-1. Incremental contribution of duration

Table S4b-2. /ai/ non-z origin-specific bridge summary

Table S4b-3. /ai/ non-z adjusted relative configuration

Table S4b-4. /au/ non-z origin-specific bridge summary

Table S4b-5. /au/ non-z adjusted relative configuration

S4c. Cross-method reconstruction summary across scales (z vs not-z)

The fPCA rows use the pooled 2 × 2 results in S3c–S3d and S4a. GAMM and FoF entries retain their respective analyses in S2 and S5. Discrepancy and RMSE values have representation-dependent units; compare the cross-origin pattern within each scale.

Taken together, the non-z analyses support the same broad interpretation as the main text: the GAMM cascade asymmetry remains, the fPCA bridge asymmetry remains, the stronger Xi'an inverse relative configuration remains, and the FoF direct full-trajectory contrast is partly scale-sensitive across representations.

Figure S4-1A. Non-z cascade reconstruction for /ai/ in the Tianjin-speaker dataset

Representative non-z direct-versus-cascade trajectory comparison for /ai/. These panels are shown across sex and variety, with carrier type fixed at 1, left segment fixed at m, and reading order set to the panel-specific reference level (S for SM, T for TJ, X for XA).

Figure S4-1B. Non-z cascade reconstruction for /au/ in the Tianjin-speaker dataset

Representative non-z direct-versus-cascade trajectory comparison for /au/. These panels are shown across sex and variety, with carrier type fixed at 1, left segment fixed at m, and reading order set to the panel-specific reference level (S for SM, T for TJ, X for XA).

Figure S4-2A. Non-z cascade reconstruction for /ai/ in the Xi'an-speaker dataset

Representative non-z direct-versus-cascade trajectory comparison for /ai/. These panels are shown across sex and variety, with carrier type fixed at 1, left segment fixed at m, and reading order set to the panel-specific reference level (S for SM, T for TJ, X for XA).

Figure S4-2B. Non-z cascade reconstruction for /au/ in the Xi'an-speaker dataset

Representative non-z direct-versus-cascade trajectory comparison for /au/. These panels are shown across sex and variety, with carrier type fixed at 1, left segment fixed at m, and reading order set to the panel-specific reference level (S for SM, T for TJ, X for XA).

Figure S4-3. Pooled bridge marginal means (non-z)

Adjusted direct and bridged F1-PC scores with 95% confidence intervals from the all-origin pooled bridge.

Figure S4-4. Adjusted relative configuration (non-z)

Centered trajectories from covariate-adjusted PC1–PC2 marginal means.

Figure S4-5. Adjusted relative similarity over time (non-z; /ai/ + /au/)

Correlations across the four tones at each of the 11 time points, calculated from covariate-adjusted PC1–PC2 reconstructions centered within origin × variety cell and domain. Negative values indicate inverse ordering or spacing.

Figure S4-6. Pooled adjusted PC time shapes (non-z)

Tone-wise adjusted PC1 and PC2 contributions in the four origin × variety cells.

Figure S4-7. /ai/ pooled bridge marginal means (non-z)

Adjusted direct and bridged scores with 95% confidence intervals.

Figure S4-8. /ai/ adjusted relative configuration (non-z)

Centered reconstructions from adjusted PC1–PC2 scores.

Figure S4-9. Adjusted relative similarity over time (non-z; ai)

Correlations across the four tones at each of the 11 time points, calculated from covariate-adjusted PC1–PC2 reconstructions centered within origin × variety cell and domain. Negative values indicate inverse ordering or spacing.

Figure S4-10. /au/ pooled bridge marginal means (non-z)

Adjusted direct and bridged scores with 95% confidence intervals.

Figure S4-11. /au/ adjusted relative configuration (non-z)

Centered reconstructions from adjusted PC1–PC2 scores.

Figure S4-12. Adjusted relative similarity over time (non-z; au)

Correlations across the four tones at each of the 11 time points, calculated from covariate-adjusted PC1–PC2 reconstructions centered within origin × variety cell and domain. Negative values indicate inverse ordering or spacing.

S5. Function-on-function (FoF) mapping

The FoF analyses treat the normalized f0 contour as a functional predictor and ask whether the whole f0 trajectory predicts the whole F1 trajectory, rather than only the concurrent relation between f0(t) and F1(t) at the same time point. The target model can be written as

where s indexes time on the input f0 curve, t indexes time on the output F1 curve, Z_i contains token-level covariates, and \(\beta(s,t)\) is a smooth coefficient surface describing how input-time deviations in f0 are associated with output-time variation in F1.

In practice, each token is converted from the original long format into an 11-point f0 curve and an 11-point F1 curve. Missing internal points are linearly interpolated, and leading or trailing gaps are filled by nearest-neighbor carry-forward/carry-backward so that each retained token has a complete grid on 0:10. The f0 matrix is then centered at each input time point, so that the FoF term is interpreted as the effect of token-level deviations from the average normalized f0 contour rather than the effect of the grand mean trajectory itself. The integral is approximated by trapezoidal quadrature, so the discretized FoF contribution becomes

The implementation uses the matrix covariates S, Tmat, and L. For each output row F1_i(t_j), S stores the full input-time grid 0:10, Tmat stores the current output time repeated across columns, and L stores the centered and trapezoidally weighted full f0 curve for that token. The FoF term is then fitted in mgcv::bam() as te(S, Tmat, by = L), together with parametric condition effects, smooths of output time and duration, factor-specific output-time smooths, a speaker-specific factor-smooth term, and an AR(1) residual structure. Cross-validation is speaker-blocked, so held-out speakers are never seen during training; predictions exclude the speaker-specific smooth so that evaluation reflects transferable whole-curve structure rather than speaker memorization.

The FoF coefficient surface shown later in this section is a visualization of \(\beta(s,t)\) on the discrete 11-by-11 grid. It is extracted from the full-data refit by probing the te(S,Tmat):L term one input time point at a time and then plotting the resulting grid as a heatmap. We verified this extraction in two formally equivalent ways: a predict(..., type = "terms") probe implementation and an lpmatrix implementation that isolates the same one-hot input location. The term-prediction and lpmatrix methods yielded identical coefficient estimates and standard errors at all 121 grid points in both datasets.

Table S5a. FoF overall prediction summary

Table S5b-1. Diagonal concentration statistics for the coefficient surfaces

Table S5b-2. Hardest and easiest FoF conditions by dataset and scale

Table S5c. FoF tensor-surface compareML summary

This table compares the reduced FoF model without the tensor-product surface to the corresponding model including te(S, Tmat, by = L). In both datasets, the tensor-product term significantly improved model fit and reduced AIC, indicating that a distributed input-output surface captures structure beyond a simpler non-surface specification.

Residualized FoF evaluation of tone-related structure

The main text supplements the direct full-trajectory FoF comparison with a residualized evaluation that targets tone-related structure more directly. The goal is to remove the shared diphthong trajectory within each diphthong-variety-by-sex-by-time cell before comparing observed and predicted curves.

The implementation proceeds cell by cell. Let g index a diphthong-variety-by-sex-by-time cell, and let tau index the four tones. For each curve type c in {Observed, Predicted}, a tone-equal shared baseline is computed as:

This baseline is tone-equal rather than token-count-weighted, so the common trajectory is not dominated by tone imbalance in the sample. Residualized values are then defined as:

The residualized tone-structure coefficient of determination is:

This statistic asks how much of the tone-structured F1 variation is recovered after the shared diphthong trajectory has been removed. The associated \(NRMSE_{\text{tone}}\) rescales the residualized RMSE by the observed residual standard deviation:

A complementary tone-separation recovery analysis summarizes how much observed tonal spacing is retained in the predicted residual curves. Within each cell g, the mean residual is first computed for each tone. Pairwise tone separation is then defined as the average of the six absolute pairwise differences among those four tone means, and the pairwise recovery ratio is:

An RMS-based separation ratio is also retained, using the root-mean-square of the four tone-specific residual means within each cell. Higher values in these recovery ratios indicate that predicted curves preserve more of the observed tone separation after removal of the common diphthong trajectory.

Under these residualized FoF measures, the Xi'an-speaker dataset shows stronger recovery of tone-related structure than the Tianjin-speaker dataset. The overall residualized \(R^2_{\text{tone}}\) is 0.023 in Tianjin and 0.124 in Xi'an, and the overall pairwise tone-separation recovery ratio is 0.365 in Tianjin and 0.616 in Xi'an.

Table S5d. Residualized FoF tone-structure summary (overall)

Table S5e. Residualized FoF tone-structure summary by condition

Table S5f. Residualized FoF tone-separation recovery (overall)

Here, the recovery-ratio columns indicate how fully observed tone separation is preserved (values closer to 1 are better), whereas the remaining columns report the absolute observed and predicted separation magnitudes.

Table S5g. Residualized FoF tone-separation recovery by condition

As in Table S5f, values closer to 1 in the recovery-ratio columns indicate better preservation of observed tone separation within each condition, while the remaining columns report the absolute observed and predicted separation magnitudes.

S5h. Robustness to omission of the 10Hz smoothing step

As a supplementary robustness check, we repeated the FoF analyses using the unsmoothed f0 input (f0_ori) instead of the 10Hz-smoothed f0_10 trajectory before within-speaker z normalization. Across both datasets, the main predictive results changed only minimally. In Tianjin, the cross-validated overall RMSE/R2 changed from 0.54423/0.69818 to 0.54405/0.69839; in Xi'an, they changed from 0.64581/0.56604 to 0.64571/0.56616. The residualized tone-structure and tone-separation summaries also showed only very small shifts (Tianjin R^2_tone: 0.02338 -> 0.02422; Xi'an R^2_tone: 0.12422 -> 0.12459), indicating that the main whole-curve conclusions are not materially dependent on the 10Hz smoothing step.

Differences in the estimated coefficient surfaces were most apparent near the temporal boundaries, despite the small changes in prediction metrics.

Figure S5h-1. FoF coefficient surface without the 10Hz smoothing step (Tianjin-speaker dataset)

Coefficient surface from the full-data FoF refit using unsmoothed f0_ori rather than the 10Hz-smoothed f0_10 input.

Figure S5h-2. FoF coefficient surface without the 10Hz smoothing step (Xi'an-speaker dataset)

Coefficient surface from the full-data FoF refit using unsmoothed f0_ori rather than the 10Hz-smoothed f0_10 input.

Figure S5-1. FoF cross-validated mean curves by tone (z scale, Tianjin-speaker dataset)

Observed-versus-predicted mean curves under speaker-blocked FoF cross-validation.

Figure S5-2. FoF cross-validated mean curves by tone (z scale, Xi'an-speaker dataset)

Observed-versus-predicted mean curves under speaker-blocked FoF cross-validation.

Figure S5-3. FoF coefficient surface (z scale, Tianjin-speaker dataset)

Coefficient surface extracted from the full-data FoF refit after the cross-validation procedure.

Figure S5-4. FoF coefficient surface (z scale, Xi'an-speaker dataset)

Coefficient surface extracted from the full-data FoF refit after the cross-validation procedure.

Figure S5-5. FoF cross-validated mean curves by tone (non-z, Tianjin-speaker dataset)

Observed-versus-predicted mean curves from the non-z FoF rerun.

Figure S5-6. FoF cross-validated mean curves by tone (non-z, Xi'an-speaker dataset)

Observed-versus-predicted mean curves from the non-z FoF rerun.

Figure S5-7. FoF coefficient surface (non-z, Tianjin-speaker dataset)

Coefficient surface extracted from the full-data non-z FoF refit after the cross-validation procedure.

Figure S5-8. FoF coefficient surface (non-z, Xi'an-speaker dataset)

Coefficient surface extracted from the full-data non-z FoF refit after the cross-validation procedure.

S6. Supplementary seven-group fPCA analyses and non-z reruns

This section retains the separate-group analyses: two speaker-group datasets, four origin-specific variety subgroups, and a combined SM dataset. These use separately estimated functional bases and supplement the primary common-basis pooled 2 × 2 analysis in S3. The z-scale results appear in S6a–S6g, followed by their non-z counterparts in S6h–S6i. Condition-specific bridge reconstructions from these separate-group analyses are collected in S7b.

S6a. Analysis blocks and token counts

S6b. Change in model fit after adding duration to score-space models

This separate-group sensitivity table supports the qualitative observation that duration contributes mainly to the dynamic F1 dimension (PC2), with little effect on the global level dimension (PC1).

S6c. Bridge reconstruction summaries

These summaries quantify how well low-dimensional f0 structure recovers low-dimensional F1 structure in the bridge analysis. The component-wise bridge \(R^2\) values indicate how much variance in F1-PC1 to F1-PC3 is explained by f0-PC1 to f0-PC3, while the mean absolute differences show how far the bridged scores remain from the directly observed F1 scores in each component. The Euclidean discrepancy measures summarize this direct-versus-bridged gap jointly in the PC1--PC2 and PC1--PC3 score spaces.

S6d. Covariate-adjusted relative-configuration summary

Figure S6-1. Group-level eigenfunctions (z scale, Tianjin-speaker dataset)

Pooled eigenfunctions for the Tianjin-speaker z-scale fPCA analysis.

Figure S6-2. Group-level eigenfunctions (z scale, Xi'an-speaker dataset)

Pooled eigenfunctions for the Xi'an-speaker z-scale fPCA analysis.

Figure S6-3. Paired f0-F1 PC effects (z scale, Tianjin-speaker dataset)

Pooled paired f0-F1 PC effect plots for the Tianjin-speaker z-scale analysis.

Figure S6-4. Paired f0-F1 PC effects (z scale, Xi'an-speaker dataset)

Pooled paired f0-F1 PC effect plots for the Xi'an-speaker z-scale analysis.

Figure S6-5. Group-level adjusted PC time-shape curves (z scale, Tianjin-speaker dataset)

Covariate-adjusted mixed-model PC reconstructions for the Tianjin-speaker pooled analysis.

Figure S6-6. Group-level adjusted PC time-shape curves (z scale, Xi'an-speaker dataset)

Covariate-adjusted mixed-model PC reconstructions for the Xi'an-speaker pooled analysis.

Figure S6-7. Group-level bridge marginal means (z scale, Tianjin-speaker dataset)

Direct-versus-bridged F1-PC marginal means for the Tianjin-speaker pooled analysis.

Figure S6-8. Group-level bridge marginal means (z scale, Xi'an-speaker dataset)

Direct-versus-bridged F1-PC marginal means for the Xi'an-speaker pooled analysis.

Figure S6-9. Adjusted relative configuration (z scale, Tianjin-speaker dataset)

Model-adjusted relative f0-F1 configuration trajectories.

Figure S6-10. Adjusted relative configuration (z scale, Xi'an-speaker dataset)

Model-adjusted relative f0-F1 configuration trajectories.

Figure S6-11. Adjusted relative similarity by time (z scale, Tianjin-speaker dataset)

By-time rank and level similarity for the adjusted relative configuration analysis.

Figure S6-12. Adjusted relative similarity by time (z scale, Xi'an-speaker dataset)

By-time rank and level similarity for the adjusted relative configuration analysis.

Table S6e. Method-context sensitivity checks for pooled bridge mixed models

This table shows that the direct-versus-bridged difference varies significantly across pooled group-level conditions for all three F1 principal components, indicating that bridge reconstruction only partially preserves the observed condition-structured score pattern.

S6f. Split-by-diphthong summaries (z scale)

Table S6f-1. Group-level split bridge summaries

Table S6f-2. Group-level split relative-configuration summaries

Figure S6-13. Split-by-diphthong bridge comparison (z scale, Tianjin-speaker /ai/)

Direct-versus-bridged F1-PC summaries from the fully split /ai/ subset of the Tianjin-speaker group.

Figure S6-14. Split-by-diphthong bridge comparison (z scale, Tianjin-speaker /au/)

Direct-versus-bridged F1-PC summaries from the fully split /au/ subset of the Tianjin-speaker group.

Figure S6-15. Split-by-diphthong bridge comparison (z scale, Xi'an-speaker /ai/)

Direct-versus-bridged F1-PC summaries from the fully split /ai/ subset of the Xi'an-speaker group.

Figure S6-16. Split-by-diphthong bridge comparison (z scale, Xi'an-speaker /au/)

Direct-versus-bridged F1-PC summaries from the fully split /au/ subset of the Xi'an-speaker group.

Figure S6-17. Split-by-diphthong relative similarity (z scale, Tianjin-speaker /ai/)

By-time spearman rank and pearson level similarity from the fully split /ai/ subset of the Tianjin-speaker group.

Figure S6-18. Split-by-diphthong relative similarity (z scale, Tianjin-speaker /au/)

By-time spearman rank and pearson level similarity from the fully split /au/ subset of the Tianjin-speaker group.

Figure S6-19. Split-by-diphthong relative similarity (z scale, Xi'an-speaker /ai/)

By-time spearman rank and pearson level similarity from the fully split /ai/ subset of the Xi'an-speaker group.

Figure S6-20. Split-by-diphthong relative similarity (z scale, Xi'an-speaker /au/)

By-time spearman rank and pearson level similarity from the fully split /au/ subset of the Xi'an-speaker group.

S6g. Multivariate f0+F1 robustness analysis

As a convergence check on the separate f0 and F1 decompositions used in the main text, a multivariate FDA was also run on paired f0 and F1 trajectories within each diphthong. This multivariate analysis uses the same 11-point normalized trajectories and the same split-by-diphthong logic as the supporting analyses above, but treats f0 and F1 as two dimensions of a single functional object rather than concatenating them or decomposing them separately.

The multivariate rerun uses f0_10Z_interp and f1Z and retains only tokens with at least five observed points in each signal before interpolation.

Table S6g-1. Split-by-diphthong explained variance in the multivariate f0+F1 analysis

Across all four split analyses, the first multivariate component accounts for roughly one half of the total joint f0+F1 variance, and the first three components together account for approximately 82--84% of the total variance. These results indicated that the shared f0+F1 structure is low-dimensional even when the two signals are modeled jointly.

Table S6g-2. Incremental contribution of duration in the split multivariate f0+F1 score-space models

The duration follow-up shows a selective pattern: duration contributes very little to the dominant shared component in most cells, but can contribute more strongly to higher multivariate components (PC3), especially for the Tianjin /ai/ split. This is compatible with the main-text interpretation that the strongest shared f0-F1 structure is not reducible to duration alone.

Figure S6-21A. Multivariate f0+F1 score space for the Tianjin-speaker /ai/ split

Shared multivariate score space for the split-by-diphthong Tianjin /ai/ analysis.

Figure S6-21B. Multivariate f0+F1 score space for the Tianjin-speaker /au/ split

Shared multivariate score space for the split-by-diphthong Tianjin /au/ analysis.

Figure S6-21C. Multivariate f0+F1 score space for the Xi'an-speaker /ai/ split

Shared multivariate score space for the split-by-diphthong Xi'an /ai/ analysis.

Figure S6-21D. Multivariate f0+F1 score space for the Xi'an-speaker /au/ split

Shared multivariate score space for the split-by-diphthong Xi'an /au/ analysis.

Figure S6-22A. Multivariate f0+F1 reconstruction for the Tianjin-speaker /ai/ split

Joint reconstruction of paired f0 and F1 trajectories from the first three multivariate components in the Tianjin /ai/ split.

Figure S6-22B. Multivariate f0+F1 reconstruction for the Tianjin-speaker /au/ split

Joint reconstruction of paired f0 and F1 trajectories from the first three multivariate components in the Tianjin /au/ split.

Figure S6-22C. Multivariate f0+F1 reconstruction for the Xi'an-speaker /ai/ split

Joint reconstruction of paired f0 and F1 trajectories from the first three multivariate components in the Xi'an /ai/ split.

Figure S6-22D. Multivariate f0+F1 reconstruction for the Xi'an-speaker /au/ split

Joint reconstruction of paired f0 and F1 trajectories from the first three multivariate components in the Xi'an /au/ split.

S6h. Seven-group non-z duration, bridge, and relative-configuration summaries

Table S6h-1. Additional contribution of duration in the non-z score-space models

Table S6h-2. Non-z bridge reconstruction summaries

Table S6h-3. Non-z covariate-adjusted relative-configuration summary

Table S6h-4. Non-z method-context sensitivity checks

S6i. Split-by-diphthong non-z summaries

Table S6i-1. Group-level split bridge summaries under non-z scaling

Table S6i-2. Group-level split relative summaries under non-z scaling

Figure S6-23. Group-level bridge marginal means (non-z, Tianjin-speaker dataset)

Direct-versus-bridged F1-PC marginal means from the non-z analysis.

Figure S6-24. Group-level bridge marginal means (non-z, Xi'an-speaker dataset)

Direct-versus-bridged F1-PC marginal means from the non-z analysis.

Figure S6-25. Adjusted relative configuration (non-z, Tianjin-speaker dataset)

Model-adjusted relative configuration from the non-z analysis.

Figure S6-26. Adjusted relative configuration (non-z, Xi'an-speaker dataset)

Model-adjusted relative configuration from the non-z analysis.

Figure S6-27. Adjusted relative similarity by time (non-z, Tianjin-speaker dataset)

By-time rank and level similarity from the non-z analysis.

Figure S6-28. Adjusted relative similarity by time (non-z, Xi'an-speaker dataset)

By-time rank and level similarity from the non-z analysis.

S7. Cross-method reconstructed F1 curves by condition

This section presents condition-specific F1 reconstructions from the three modelling frameworks. Each panel shows one variety × diphthong condition, with the four tones distinguished by color. The plots provide descriptive comparisons of reconstructed trajectory shapes.

S7a. GAMM cascade reconstructions

These GAMM panels use the cascade prediction only, under the reference-style setting used in the follow-up script (leftSegment = m, type = 1, with order held at the script's reference value for each variety).

Figure S7-1. Cascade reconstructed F1 curves by condition in the Tianjin-speaker dataset

Group-level cascade reconstructions from the z-scale analysis, shown by variety × diphthong condition with tones overlaid.

Figure S7-2. Cascade reconstructed F1 curves by condition in the Xi'an-speaker dataset

Group-level cascade reconstructions from the z-scale analysis, shown by variety × diphthong condition with tones overlaid.

S7b. fPCA bridge reconstructions

These figures are reconstructed from the separate-group, split-by-diphthong bridge outputs summarized in S6f; the primary pooled 2 × 2 bridge results are in S3c and S3e. Each panel corresponds to one variety × diphthong condition, so the two varieties within each group analysis are separated.

Figure S7-3. fPCA bridge reconstruction by tone

Bridge-reconstructed F1 curves for the SM /ai/ condition in the Tianjin-speaker group analysis.

Figure S7-4. fPCA bridge reconstruction by tone

Bridge-reconstructed F1 curves for the TJ /ai/ condition in the Tianjin-speaker group analysis.

Figure S7-5. fPCA bridge reconstruction by tone

Bridge-reconstructed F1 curves for the SM /au/ condition in the Tianjin-speaker group analysis.

Figure S7-6. fPCA bridge reconstruction by tone

Bridge-reconstructed F1 curves for the TJ /au/ condition in the Tianjin-speaker group analysis.

Figure S7-7. fPCA bridge reconstruction by tone

Bridge-reconstructed F1 curves for the SM /ai/ condition in the Xi'an-speaker group analysis.

Figure S7-8. fPCA bridge reconstruction by tone

Bridge-reconstructed F1 curves for the XA /ai/ condition in the Xi'an-speaker group analysis.

Figure S7-9. fPCA bridge reconstruction by tone

Bridge-reconstructed F1 curves for the SM /au/ condition in the Xi'an-speaker group analysis.

Figure S7-10. fPCA bridge reconstruction by tone

Bridge-reconstructed F1 curves for the XA /au/ condition in the Xi'an-speaker group analysis.

S7c. FoF predicted-only reconstructions

These FoF panels are based on the cross-validated predicted curves only.

Figure S7-11. FoF predicted-only reconstruction by tone

FoF predicted-only F1 curves for the SM /ai/ condition in the Tianjin-speaker dataset.

Figure S7-12. FoF predicted-only reconstruction by tone

FoF predicted-only F1 curves for the TJ /ai/ condition in the Tianjin-speaker dataset.

Figure S7-13. FoF predicted-only reconstruction by tone

FoF predicted-only F1 curves for the SM /au/ condition in the Tianjin-speaker dataset.

Figure S7-14. FoF predicted-only reconstruction by tone

FoF predicted-only F1 curves for the TJ /au/ condition in the Tianjin-speaker dataset.

Figure S7-15. FoF predicted-only reconstruction by tone

FoF predicted-only F1 curves for the SM /ai/ condition in the Xi'an-speaker dataset.

Figure S7-16. FoF predicted-only reconstruction by tone

FoF predicted-only F1 curves for the XA /ai/ condition in the Xi'an-speaker dataset.

Figure S7-17. FoF predicted-only reconstruction by tone

FoF predicted-only F1 curves for the SM /au/ condition in the Xi'an-speaker dataset.

Figure S7-18. FoF predicted-only reconstruction by tone

FoF predicted-only F1 curves for the XA /au/ condition in the Xi'an-speaker dataset.