T Paulauskas

B Čechavičius

V Karpus

L Jočionis

S Tumėnas

J Devenson

V Pačebutas

S Stanionytė

V Strazdienė

A Geižutis

M Čaplovičová

V Vretenár

M Walls

M Kociak

A Krotkus

Spontaneous atomic ordering of dilute GaAsBi bismides: Structural and optical study

Keywords: GaAsBi, CuPt ordering, polarized photoluminescence, spectroscopic ellipsometry, birefringence, linear dichroism, atomic resolution STEM

The investigated dilute GaAs 1-x Bi x (x = 0.034-0.052) bismides have been MBE-grown on exact and low-angle offcut (001) GaAs as well as on (001) Ge substrates. Atomic-resolution scanning transmission electron microscopy and X-ray diffraction measurements show the presence of spontaneous CuPt B -type ordering in all the bismide samples. The ordering lifts the valence band degeneracy and induces the optical polarization anisotropy, which manifests in photoluminescence and transmittance spectra of the bismides samples, and is also confirmed by birefringence and linear dichroism measurements. The ordering-induced valence band splitting in dilute GaAs 1-x Bi x is estimated to be much larger than that in the conventional III-V semiconductor alloys (e.g. GaInP).

Introduction

The dilute GaAs 1x Bi x alloy is thus far the most extensively investigated material among Bicontaining III-V semiconductors [START_REF] Batool | Bismuth-containing III-V semiconductors: Epitaxial growth and physical properties[END_REF][START_REF] Wang | Bismuth-Containing Alloys and Nanostructures[END_REF]. The main interest in this alloy, here forth referred to as bismide, arises due to its potential and already implemented applications in the near-to mid-infrared range optoelectronics, lasers, photodetectors, solar cells (see, e.g., [START_REF] Hepp | MOVPE growth and device applications of ternary and quaternary dilute bismide alloys on GaAs Substrates[END_REF]). Despite wide-scale physical properties investigations of bismides, which started in the 1990s, only rather recently, in the early 2010s, first reports were made on the observation of spontaneous CuPt-type atomic ordering in GaAsBi [46].

The CuPt-type atomic ordering in conventional III-V alloys, most notably in Ga x In 1x P and GaAs 1x Sb x , is well known from its discovery in the late 1980s [START_REF] Mascarenhas | Spontaneous Ordering in Semiconductor Alloys[END_REF]. It is widely accepted that the CuPt-ordering results from surface reconstructions during epitaxial growth, and the (21) reconstruction, which is assumed to be a prerequisite for the ordering, is also observed in bismides [811].

A perfectly CuPt B -type ordered zinc-blende AB 0.5 C 0.5 alloy is a superlattice of alternating AB and AC monolayers on or planes, conventionally termed as B + and B  subvariants. Ordering on the other two planes, and , is termed CuPt A -type, but is rarely observed in practice. The CuPt B ordering occurs in domains with no-preference for B + and B  subvariants if an epitaxial layer is grown on exact (001) substrate. A single subvariant throughout the layer can be stabilized by employing (001) low-angle offcut wafers towards direction. This typically leads to larger ordered domains, and a single subvariant ordering has also been confirmed in GaAsBi [START_REF] Paulauskas | Atomic-resolution EDX , HAADF , and EELS study of GaAs 1-x Bi x alloys[END_REF][START_REF] Paulauskas | GaAs 1-x Bi x growth on Ge: anti-phase domains, ordering, and exciton localization[END_REF]. The ordering is never perfect and is quantified by the ordering parameter η. Partially ordered alloy AB 1x C x can be described as the AB 1(xη/2) C xη/2 and AB 1(x+η/2) C x+η/2 monolayers on or planes, with a maximal value of ordering parameter η max = min[2x, 2(1x)].

The ordered conventional III-V alloys are usually investigated at the commensurate compositions x  0.5 where the ordering parameter can approach unity. On the other hand, the dilute bismides have been mostly investigated at the atomic Bi concentrations close to

x  0.05, and an expected value of ordering parameter is one order of magnitude lower,  ~ 0.1, than in the conventional semiconductors. Another distinctive feature of bismides is an epitaxial strain of GaAsBi layers due to a lattice mismatch with GaAs substrates, while the conventional ordered III-V alloys are usually grown on the lattice-matched substrates, namely

GaInP on GaAs and GaAsSb on InP.

Previous studies of the atomic ordering of bismides focused primarily on their atomic and microstructure analysis [START_REF] Norman | Atomic ordering and phase separation in MBE GaAs 1-x B x[END_REF][START_REF] Reyes | Bismuth incorporation and the role of ordering in GaAsBi/GaAs structures[END_REF][START_REF] Bastiman | Bi incorporation in GaAs(100) -2  1 and 4  3 reconstructions investigated by RHEED and STM[END_REF][START_REF] Paulauskas | Atomic-resolution EDX , HAADF , and EELS study of GaAs 1-x Bi x alloys[END_REF][START_REF] Paulauskas | GaAs 1-x Bi x growth on Ge: anti-phase domains, ordering, and exciton localization[END_REF]. However, no relationships between the ordering and physical properties have been established in bismides as of this writing. The present article expands the studies by examining the influence of the ordering on GaAsBi optoelectronic properties.

This article is organized as follows. Firstly, in section 2, a comprehensive structural analysis of the investigated bismides samples is presented. Polarized photoluminescence and spectroscopic ellipsometry results are reported and discussed in section 3. The article ends with summarising conclusions and details of the experimental techniques employed in this study. Six GaAsBi samples, labelled S1S6, were grown for this study, all epitaxially deposited by molecular-beam epitaxy (MBE, see Methods for more details). Samples S1S3 were deposited on exact (001) GaAs substrates, while GaAs substrates with 2° and 6° offcuts towards direction were used for samples S4 and S5, respectively. Exact (001) Ge substrate was employed for sample S6, which also included a nominally 50 nm thick GaAs buffer layer deposited on Ge. The sample information is summarized in Table 1.

Atomic structure analysis

Sample characterization

X-ray diffraction rocking-curve (XRD-RC) measurements were used to examine the structural quality of bismide layers and determine their thickness [see Fig. 1(a)]. As can be seen from the sharp peaks and resolved thickness fringes, high quality epitaxial layers were achieved in samples S1S5. XRD reciprocal-space mapping (RSM) showed that samples S1S4 were fully biaxially compressively strained to GaAs substrates, while sample S5 showed a minor, ~ 2 %, lattice relaxation. GaAsBi is known to be resilient to lattice relaxation and layers well in excess of the critical thickness for dislocation nucleation are often observed [START_REF] Batool | Bismuth-containing III-V semiconductors: Epitaxial growth and physical properties[END_REF][START_REF] France | In situ strain relaxation comparison between GaAsBi and GaInAs grown by molecular-beam epitaxy[END_REF]. XRD-RSM was not performed on S6 since the sample was of polycrystalline nature. The bismuth concentrations, determined from XRD measurements, are presented in Table 1.

The band gap values E g listed in Table 1 are based on the room-temperature band-edge PL peak positions. Sample surfaces were examined using atomic-force microscopy (AFM) and scanning electron microscopy (SEM). Surface roughness values, obtained using AFM, are shown in Table 1 and indicate high surface quality of the S1S5 samples. Sample S6 surface was rough and showed high-density of Bi droplets. Apart from S6, only S2 had surface droplets with a size of about 50 nm, albeit with a low density of less than one particle per µm 2 . The droplets on S2 were found to be composed of Ga, as was determined using X-ray energy dispersive spectroscopy (EDX). Sample S5, on the other hand, displayed circular pits with a density of ~ 2 per µm 2 . Among samples S1S5, only S3 displayed phase-separated, elongated GaAs-like domains within the bismide layer as was revealed by the cross-sectional scanning transmission electron microscopy (STEM). The phase-separation can be associated with higher surface roughness in S3. No Bi segregation was observed in S1S5 samples.

Spontaneous atomic ordering

The so-called skew-symmetric XRD-RC measurements were performed to probe the presence of ½{111}-type superlattice reflections associated with CuPt-type ordering [START_REF] Liu | Determination of CuPt-type ordering in GaInP by means of x-ray diffraction in the skew, symmetric arrangement[END_REF]. Sample S4 grown on 2° offcut substrate manifested the largest intensity of superlattice reflections. We note that no MBE growth optimizations were undertaken for investigated samples to maximize the CuPt-type ordering. Our previous atomically-resolved STEM-EDX ordering parameter evaluation showed the η = 0.07 value for a sample with 5.5 % Bi [START_REF] Paulauskas | GaAs 1-x Bi x growth on Ge: anti-phase domains, ordering, and exciton localization[END_REF].

Although the ordering parameters for the present S1-S6 samples were not examined by EDX, one can expect them to be of the order of η ~ 0.1.

Cross-sectional structural investigations of S1-S6 samples were performed using aberration-

Optical spectroscopy

Scheme of optical transitions

The spontaneous CuPt-type atomic ordering in the conventional zinc-blende alloys and its influence on their electronic structure and optical properties has been extensively studied by Zunger, Wei, Mascarenhas, Zhang, and others (see, e.g., [7, 1719]). The CuPt-type ordering changes the cubic zinc-blende point symmetry T d to the trigonal C 3v one. The perturbative part of the Hamiltonian, describing an influence of the ordering, has the same form as that of the uniaxial strain along the ordering direction [START_REF] Wei | Optical properties of zinc-blende semiconductor alloys: Effects of epitaxial strain and atomic ordering[END_REF], namely the axis for the CuPt B+ subvariant [see Fig. Due to the ordering-induced changes in zinc-blende point symmetry, the PL intensity should be polarization-dependent even at the light propagation along the sample-surface normal

[001]. Indeed, the E 0,1 optical transitions between the conduction and hh bands can be induced by the light, electric vector of which is perpendicular to the ordering direction,

. The E 0,2 optical transitions between the conduction band and the lower, lh, valence band occur at the polarization only. Therefore, the PL polarization dependence at the surface-normal light propagation will be determined by a variation of the projection of E-vector onto the B + ordering plane . The PL intensity should be maximal at the polarization and minimal for , because the direction lies in the plane, while the direction projects onto the plane with the largest angle [see Fig. 4(a)]. Therefore, the polarized PL ratio theoretically should be . When a bismide layer is grown on GaAs substrate miscut towards direction, the ordering axis rotates further away from the sample surface normal, and the PL polarization ratio becomes , which is 2.7 and 2.3 for typical 2 and 6 offcuts, respectively [START_REF] Zhang | Effects of strain, substrate misorientation, and excitonic transition on the optical polarization of ordered zinc-blende semiconductor alloys[END_REF].

The normal-PL polarization dependence can be suppressed by the epitaxial strain, which also lifts the valence band degeneracy of zinc-blende semiconductors. The strained layer becomes optically uniaxial, with the optical axis perpendicular to the layer. Hence, in the presence of epitaxial strain (and in the absence of ordering), no polarization dependence is expected for the normal light propagation.

Polarized photoluminescence

The polarized PL measurements along a sample-surface normal were performed at room (RT, 300 K) and liquid nitrogen (80 K) temperatures. Representative PL spectra of samples S1, S4, and S6 are-+ shown in Fig. 5. The PL intensity of GaAsBi with increasing Bi concentration tends to be weaker and hence PL spectra become noisier. This is most evident in highly defective and non-uniform lattice of S6, which shows a wide and noisy RT PL peak.

Generally, GaAsBi alloys manifest pronounced exciton localization effects, which usually are traced in the temperature dependence of PL. Typical feature is the so-called PL "s-shift" due to exciton hopping and trapping at localized Bi-related density of states near the valence band maximum [START_REF] Imhof | Clustering effects in Ga(AsBi)[END_REF][START_REF] Wilson | Assessing the nature of the distribution of localised states in bulk GaAsBi[END_REF]. In our samples S2S5 we observed defect-related RT PL peak below the RT band-edge emission, at about 0.9 eV [see Fig. 5(b)], which disappears at 80 K.

As seen from Fig. 5, the samples grown on GaAs substrates consistently show higher intensity for the polarization as compared to the one, while no polarization dependence is observed for sample S6. Photoluminescence polarizations parallel to directions 45° and 135° from [110] were also tested, showing the same PL intensities for each of all S1S6 samples. The observed PL polarization dependence is in accord with that expected for the CuPt B -type ordered bismide layers under the assumption that PL is predominantly determined by the E 0,1 optical transitions. Since the S6-sample is composed of crystallites whose ordering B-axes are rotated in-plane by right angles, no polarization dependence is expected and none is observed.

The RT polarized PL ratio in S1S5 bismide samples is the range of about 1.62.4. At liquid nitrogen temperature, values remain in a similar range of about 1.41.9.

It should be noted that the 80 K spectra show a complex inner structure. This could be a manifestation of the "s-shift", however, detailed temperature-dependence investigations are needed for a definite answer.

The experimental values of R-parameter, ~ 2, are close to the theoretical one, R = 3.

However, the R = 3 value corresponds to unstrained epitaxial layers, which is obviously not the case here. The strains and epitaxial splitting Δ e in bismide layers can be easily calculated at given values of elastic constants C 12 /C 11 = 0.46, the deformation potential b = 2 eV (which for dilute GaAsBi bismides can be taken to be equal to those of GaAs [START_REF] Vurgaftman | Band parameters for III-V compound semiconductors and their alloys[END_REF]), and the a GaAs = 5.653 Å and a GaBi = 6.324 Å lattice constants [START_REF] Feng | New III-V semiconductor InGaAsBi alloy grown by molecular beam epitaxy[END_REF]. The epitaxial splitting is about 45 meV at the x = 0.05 bismuth concentration [reikia cia ideti references kas yra publikuota].

Since the epitaxial strain suppresses a polarization dependence of normal-PL, one can expect that the ordering-induced splitting in the investigated bismide samples is of the order of 40 meV.

Spectroscopic ellipsometry

We assume that the observed photoluminescence spectra are predominantly determined by the E 0,1 optical transitions to the upper valence band branch. To trace the E 0,2 transitions, we have carried out the optical transmission measurements. Figures 6(a)-6(c) present the normal-incidence transmittance, T, and optical density, A = ln(T -1 ), spectra of samples grown on GaAs substrates, the energy gap of which is above that of bismide layers and, therefore, allow for measurements in transmission geometry. Arrows in Figs. 6(a)-6(c) indicate the spectral positions of the E 0,1 and E 0,2 critical points, which were determined from positions of peaks in the optical density differential spectra, dA/d. At the polarization, the optical density experiences a sharp increase at the photon energies of about 1.17, 1.06, and 1.07 eV for samples S1, S3, and S4, respectively, which coincide with the spectral positions of PL peaks (see E g values in Table 1) with an accuracy of a few percent. At the polarization, an absorption at E 0,1 is weaker, in accordance with results of polarized-PL measurements. A sharp increase of optical density begins at the higher photon energies of about 1.23, 1.12, and 1.13 eV (for samples S1, S3, and S4, respectively), which correspond to the expected onsets of the E 0,2 optical transitions from the lower valence band branch to the conduction band.

The carried out phenomenological determination of the E 0,1 and E 0,2 energies from spectral positions of peaks in optical density derivative spectra allows for an estimation of the valence band splitting  = E 0,2  E 0,1 , which in the investigated S1-S5 samples is of about 60 meV.

This result confirms our PL-based assumption that the atomic ordering induced splitting  c is of the order of 40 meV. This estimate is close to experimental values of the CuPt-type ordering-induced valence band splitting  c in the conventional III-V semiconductor alloys.

Namely, in Ga 1x In x P alloy, at x  0.5, the  c value is close to 40 meV when the ordering parameter  approaches unity [START_REF] Zhang | Dependence of the band structure on the order parameter for partially ordered Ga x In 1x P alloys[END_REF]. However, the GaInP ordering-induced splitting is only ~ 2 meV when  = 0.1, i.e., at the maximal values of the ordering parameter in dilute GaAs 1x Bi x bismides (x ~ 0.05). That means that atomic ordering induced changes in the electronic structure of bismides are at least one order of magnitude higher than those in the conventional semiconductors (when they are compared at the same values of the ordering parameter).

To check the influence of atomic ordering on the optoelectronic properties of bismides, we have carried out the Mueller-matrix spectroscopic ellipsometry measurements of the birefringence and linear dichroism . The experimental results obtained are presented in Figs. 6(d)-6(f). As it is seen, the and spectra manifest the Lorentzian-type dispersion centered at the photon energies of E 0,1 and E 0,2 optical transitions. The experimental spectra indicate as well a feature above E 0,2 , which could be related to the ordering-induced changes to the spin-orbit split off valence band. However, we

were not able to access the optical feature in detail due to absorption in GaAs substrate.

The transmission-mode spectroscopic ellipsometry measurements of transmittance and birefringence-dichroism spectra are comparatively easy to perform and can be used as a test for the presence of CuPt-type ordering. We have carried out the test on a dozen of bismide samples with the result that all of them manifested the ordering-induced polarization anisotropy (not shown here). Actually, we have not succeeded in finding a bismide sample which possessed crystalline lattice without the presence CuPt-type ordering anisotropy effects. This preliminary observation suggests that the spontaneous atomic ordering in bismides is induced at their standard growth conditions and is ubiquitous in GaAsBi alloys.

Conclusions

In

Methods

MBE.

Samples S1-S5 were synthesized using Veeco GENxplor MBE chamber, employing metallic Ga and Bi sources, and a two-zone valved-cracker using As 2 . The growth was performed on semi-insulating GaAs substrates supplied by the AXT, Inc. The oxide removal from GaAs substrates was done by heating at 590°C in the MBE chamber. Temperature measurements were done using a thermocouple. Samples S1-S5 contained a standard-grown 100 nm GaAs buffer layer. The substrate temperature was reduced to 330°C and GaAsBi layer was deposited using As 2 /Ga beam-equivalent pressure (BEP) ~1.1 and Bi flux 6*10 -8

Torr. The (2x1) surface reconstructions were clearly seen during the growth.

Sample S6 was synthesized using SVT-A MBE chamber equipped with the same sources and

As-cracker as in the Veeco chamber. A p-type doped exact (001) Ge substrate was employed, native oxides from the substrate were removed at 640°C. The first part of GaAs buffer layer was grown at 400°C using migration-enhanced epitaxy (MEE), which was approximately 5 nm thick. Then the substrate temperature was increased to 500°C and second GaAs buffer part consisting 50 nm of GaAs was grown using standard MBE growth conditions. Finally, the GaAsBi film was grown at 350°C. The BEP ratio of As 2 to Ga was ~1.1, and Bi/Ga ratio was 0.35-0.37. Before the growth was initiated, the shutter was opened to give a bismuth cover GaAs buffer surface. AFM. Surface morphology of the samples were characterized by atomic force microscopy measurements, using a Dimension 3100 SPM system with a Nanoscope IVa controller (Veeco Instruments Inc., USA).

Janas --

SEM.

The morphology and elemental composition of the samples were carried out in FE-SEM-FIB HELIOS Nanolab 650 (FEI) arranged with an X-ray energy dispersive spectrometer (EDX) from Oxford Instruments.

PL. Photoluminescence measurements were carried out using the DPSS laser (532 nm) as an excitation source, with the intensity kept at ~ 0.1 kW/cm 2 , the 0.4 m monochromator SPM-2, and the thermoelectrically-cooled InGaAs photodetector. To minimize polarization effects of monochromator grating, the sample axes [110] and were rotated relative to the monochromator slit by 45. The GlanTaylor prism (with extinction ratio < 5•10 -6 ) was used as polarization analyser and was placed in front of the monochromator slit. The residual prism polarization was measured to be of ~ 2 %, and was taken into account by correcting recorded PL intensities. To correct for an observed degradation of GaAsBi PL intensity, the attenuation curves were acquired for each sample and the subsequent polarized PL measurements were corrected accordingly.

Spectroscopic ellipsometry. The polarized transmittance spectra and spectra of normalized Mueller matrix components were recorded with a dual rotating-compensator ellipsometer RC2 (J.A. Woollam, Co., Inc.) at the normal light incidence. After calibration, a straightthrough measurement of a Mueller matrix in the air was performed. All recorded Mueller matrix components were found to be within 0.005 of the ideal value. The linear in-plane birefringence and linear dichroism spectra were determined from the experimental Mueller matrix spectra by the analytic inversion method [START_REF] Arteaga | Analytic inversion of the MuellerJones polarization matrices for homogeneous media[END_REF].

Figure 1 .

 1 Figure 1. (a) X-ray diffraction rocking-curve (004) reflections showing high-quality S1S5 samples with distinct thickness fringes. Sample S6 is more defective. (b)(d) Skew-symmetric XRD measurements of four ½{111}-type superlattice reflections. Samples S1 and S4 show CuPt B ordering on two and one B-type planes, respectively. Sample S6 manifests superlattice reflections for all four measurement geometries.

Figures 1 (

 1 Figures 1(b)1(d) provide representative XRD-RC data of samples S1, S4, and S6. The measurements on samples S1S3, which were grown on exact (001) GaAs substrates, showed the presence of both ½ and ½ reflections, indicating B + and B  subvariants of CuPt-ordering in the samples [see the representative Fig. 1(b)]. The peak intensities were slightly larger in the thicker sample S1. The samples grown on low-angle offcut GaAs substrates, S4 and S5, demonstrated [Fig. 1(c)] only the ½ reflection, suggesting that a single subvariant CuPt B ordering is stabilized. The sample S6, deposited on exact (001) Ge substrate over a thin GaAs buffer layer, manifested the ½{111}-type superlattice reflections for all four measurement geometries [Fig. 1(d)]. This, however, does not imply that both CuPt B and CuPt A coexist in this sample. Instead, sample S6 is composed of anti-phase domains, in which the CuPt B ordering axes and are rotated by right angles, as will be shown below using STEM data.

  Fig. 2 and Fig. 3) were aligned such that the growth axis was vertical and directed towards the top of the page. The CuPt B -type ordering was detected in STEM images only along [110] zone-axis since the ordered and planes are then viewed edge-on. Bismuth containing columns can be well distinguished in HAADF images due to their higher contrast, as Bi atoms scatter probe electrons much stronger than Ga and As atoms. A representative HAADF image of the GaAsBi samples S4-S5, deposited on offcut GaAs substrates, is shown in Fig. 2(a) and indicates a single variant CuPt B+ ordering. The intensity modulations along direction are clearly visible. Fourier transform (FT) of the image shows a single pair of

Figure 2 .

 2 Figure 2. STEM cross-sectional sample analysis. (a)-(b) Atomic-resolution HAADF images of GaAsBi sample S4 grown on offcut substrate and sample S1 grown on exact GaAs substrate. (c) Lowmagnification LAADF image of sample S6 showing numerous anti-phase boundaries indicated by yellow dashed lines. (d) Atomic-resolution HAADF image of S6 near a kinked {110}-type anti-phase boundary. Note the reversal of dumbbell polarity across the boundary.

Figure 3 .

 3 Figure 3. (a) Image derived from S6-sample HAADF micrograph near anti-phase boundary (APB) using pairs of ½{111} reflections, which are outlined in the Fourier transform inset. (b) Atomic model of the APB region. A region containing two GaAsBi APDs near the GaAs buffer layer is shown in Fig. 3(a). The figure was produced from a HAADF image of S6-sample by applying two pairs of circular masks on the superlattice ½ -type and ½ -type reflections, shown in the FT inset, and forming an image from them. CuPt B -type ordering is thereby highlighted in the left-hand side GaAsBi grain, while no ordering can be seen on the right-hand side. This is because the [001] polarity reversal results in a rotation of B + and B  domains by a right angle around the vertical growth axis. The CuPt B ordered planes in such rotated domains are not viewed edgeon and thus are not visible in the HAADF image. This has been confirmed by imaging the sample along the two orthogonal -type directions. To clarify the reversal, an atomic

Figure 4 .

 4 Figure 4. (a) Atomic structure of CuPt B+ -type ordered GaAsBi. (b) Scheme of optical transitions near the fundamental bandgap.

  photoluminescence (PL) will be determined by the E 0,1 optical transitions to the upper, hh, valence band branch.

Figure 5 .

 5 Figure 5. Polarized-PL spectra of selected samples at room and liquid nitrogen temperatures.

Figure 6 .

 6 Figure 6. Transmission-mode spectroscopic ellipsometry data of selected samples S1, S3, and S4. (a)-(c) Transmittance (full curves) and optical density (dashed curves) spectra at the and polarizations. (d)-(f) Birefringence and linear dichroism spectra.

  summary, we have carried out the structural and optical study of the dilute GaAs 1x Bi x (x = 0.0340.052) bismides. The atomic structure measurements by cross-sectional STEM and XRD techniques directly revealed the spontaneous CuPt B -type atomic ordering. The epitaxial bismide layers, which were MBE-grown on exact GaAs (001) substrates, manifest both B + and B  variants of CuPt B ordering. When bismides are grown on 2 and 6 offcut GaAs substrates, a single variant CuPt B+ ordering is stabilized. The bismide layer grown on an exact (001) Ge substrate is polycrystalline in nature and is comprised of CuPt B ordered domains rotated by 90 with respect to each other around the growth axis.The spontaneous CuPt-type ordering lifts the degeneracy of the valence band and induces the polarization anisotropy of bismide layers at the normal light propagation. The anisotropy was revealed in polarized photoluminescence and transmittance spectra as well as by birefringence and linear dichroism measurements carried out by Mueller-matrix spectroscopic ellipsometry technique. The ordering-induced valence band splitting in bismides is at least one order of magnitude larger than that in the conventional III-V semiconductor alloys. The birefringence/dichroism-based optical test carried out on a dozen of various bismide layers allows for a preliminary assumption that the spontaneous atomic ordering in bismides grown at their standard conditions is wide spread. These findings are important for the future development of GaAsBi-based optoelectronics and also suggest new means to analyze the structurally complex bismide alloys.

  Samples S2-S5.XRD. X-ray diffraction data was measured by a Rigaku Smartlab diffractometer equipped with kW rotating Cu anode X-ray generator (Kα1 line), scintillation detector SC-70, and a linear Dtex/Ultra detector. Reciprocal-space mapping of (004) and (115) reflecitons was used to evaluate layer relaxations, lattice parameters, and Bi concentrations, using lattice constants: GaAs=5.653 Å, GaBi = 6.324 Å. Superlattice ½{111}-type reflections were acquired in skew-symmetric configuration.STEM.Cross-sectional lamella samples nominally 20-25 nm thick were prepared using focused-ion beam beam, employing FEI Helios Nanolab 650. Aberration-corrected JEOL-JEM ARM200CF scanning transmission electron microscope operated at 200 kV was used for sample imaging. Probe convergence semi-angle was set to 22 mrad, high-angle annular dark-field (HAADF) detector inner and outer semi-angles were set to 90 and 370 mrad, respectively. For low-angle annular dark-field, the inner and outer detector semi-angles were 40 and 160 mrad, respectively. All data were processed with Gatan Digital Micrograph GMS3.

Table 1 .

 1 Parameters of investigated samples: the type of substrate, GaAsBi layer thickness d, bismuth concentration, room-temperature bandgap E g , and sample surface roughness R q . GaAs substrate offcut angles are indicated for samples S4 and S5.

	sample	substrate	d (nm)	Bi %	E g (eV)	R q (nm)
	S1	GaAs (001)	485	3.4	1.15	0.3
	S2	GaAs (001)	192	4.8	1.08	0.4
	S3	GaAs (001)	190	5.2	1.04	1.5
	S4	GaAs (001) 2	186	5.0	1.06	0.3
	S5	GaAs (001) 6	180	4.5	1.09	0.5
	S6	Ge (001) / GaAs buffer	500	4.6	1.05	> 5

Funding

This work was supported by European Regional Development Fund (project No 01.2.2-LMT-K-718-02-0020) under grant agreement with the Research Council of Lithuania (LMTLT). This project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement No 823717 -ESTEEM3. This study was also supported by the Operational Program Research and Innovation for the project ITMS2014+ code: 313010T598 co-financed by the European Regional Development Fund. Mathieu

Competing interests

The authors declare that they have no competing interests.


Jean-Théo Jiolat

Christophe D'alessandro

Jean-Loic Le Carrou

Jose Antunes

Towards a physical model of the clavichord

String excitation by the tangent in the clavichord is a unique mechanism. The tangent, keeping in contact with the string after the initial strike, controls continuously the string tension. Four main flexible subsystems are considered in the clavichord: the tangent/key subsystem, the strings subsystem, the bridge-soundboard subsystem, the string damper subsystem. A modal description of the dynamics of these subsystems is proposed. Parameters of the subsystems are estimated on a copy of a historical instrument by Hubert (1784). The different subsystems and their couplings are modeled using a modal Udwadia-Kalaba formulation. The string-tangent interaction is modeled via the intermittent contact dynamics, using the Kirchoff-Carrier string model. Realistic string, soundboard and tangent motions are obtained using a time-domain synthesis scheme that computes the dynamics of the uncoupled subsystems and the constraints resulting from coupling between them. Simulated motions of the model in response to a driving force on the key are analysed. The results are consistent with experimental measurements and published data on the dynamics of the clavichord. The model is able to reproduce the main acoustic features of the instrument: force on the key for intonation control, key velocity for dynamic nuances control, constant spectral slope for varying dynamic nuances.

I. INTRODUCTION

The clavichord is highly prized as a practice instrument among keyboard players, because of its superior ability to encourage a polished technique [START_REF] Crowell | Every players first grammatica[END_REF] . Its sound is weak and sometime a little disappointing at the first contact with the instrument. One explanation for the special appreciation of the clavichord as a wonderful coach of finger technique lies in the refined string control allowed by its simple and direct action. In the clavichord, the string and tangent stay in contact, i.e. are mechanically coupled as long as the key remains depressed (see [START_REF] Brauchli | The clavichord[END_REF] for a thorough presentation of the instrument). This feature has important consequences on the sound and the dynamic of the instrument, allowing for expressive pitch control, a unique feature among stringed keyboard instruments. The aim of the present research is to develop a physical model of the clavichord that is able to account for, and synthesize the specific feature of the clavichord's action, from finger motion to soundboard vibration.

Relatively few acoustic studies on the clavichord have been published so far. The first ones [START_REF] Trendelenburg | Zur Klangwirkung von Klavichord, Cembalo und Flügel[END_REF][START_REF] Hands | A scientific approach to the clavichord[END_REF] mainly report descriptions of sound features (level, spectrum) due to the tangent action compared to the piano and harpsichord actions. A study on the clavichord touch and action is developed in [START_REF] Bavington | Keylever, tangent and string-a preliminary analysis of clavichord touch and action[END_REF] , showing that hardness of touch and pitch stability are related to string tension and key balance parameters. Some aspects of the physics of the instrument are investigated in [START_REF] Thwaites | Some notes on the clavichord[END_REF] : soundboard and cavity coupling, tangent velocity profile and string displacement, sound decay rate, string pair coupling effects. A linear string model is used for qualitative explanation of the tangent velocity profile (modeled as an exponential decaying function), string motion and sound decay rate. Based on these results a comparison with the piano and harpsichord is derived in [START_REF] Donald | String excitation: piano, harpsichord and clavichord[END_REF] . Focusing on the string excitation mechanism, the dynamics of the clavichord is revisited in [START_REF]On the dynamics of the clavichord: From tangent motion to sound a[END_REF] . A mass-spring-damper model of the key/tangent and string system oscillation is developed. Using a quasi-static approximation of the string motion and a delay-line model of the string, the excitation dynamics is studied and compared to experimental data. A linear relationship between tangent velocity and sound pressure level is found. String displacement (and then variation of string tension) has a significant effect on fundamental frequency, but the sound spectral slope does not vary much with tangent velocity or displacement. Measurements on 4 clavichords are in good agreement with these findings [START_REF] Christophe D'alessandro | Tonal portrait of the clavichord[END_REF] . Musical consequences of the needed joint tangent displacement/velocity controls for the clavichord playing technique are studied in terms of the so-called "clavichord's paradox" [START_REF]Le paradoxe du clavicorde et la technique de bach au clavier[END_REF] .

As for sound synthesis, two approaches for physical modeling of the clavichord have been published so far. The signal processing approach by [START_REF] Välimäki | Commuted waveguide synthesis of the clavichord[END_REF] is based on commuted wave-guide synthesis: string/tangent interaction is partly based on sampling of real sounds (for a realistic knock sound) and partly on additional filters accounting for the variable string tension during a tone. Physical models of the piano gave birth to a successful software pi- ano synthesizer able to produce realistic sound synthesis for modern and historical pianos, for harpsichords as well as for hybrid instruments and for the clavichord 12 . For commercial reasons, the physical modeling techniques used and specific features for the clavichord model are not published, but it seems that a hard metallic hammer is simulated for the initial tangent strike, together with an additional after-touch effect accounting for the enforced tangent/string contact. Note that other physical piano models based on finite element techniques have been developed and published [START_REF] Chabassier | Modélisation et simulation numérique d'un piano par modèles physiques[END_REF][START_REF] Chaigne | Modélisation du piano et couplage cordeschevalet[END_REF] .

The aim of this paper is to develop a physical model of the clavichord with special attention to the collision and contact between the tangent and the string, which is of prime importance for explaining the instrument's dynamics and sound features. Intermittent contacts are often encountered in musical instruments and have been studied in numerous works in the past decade. Collisions (see [START_REF] Banerjee | Historical origin and recent development on normal directional impact models for rigid body contact simulation: A critical review[END_REF] for a review on the subject ) are modeled by a penalty approach [START_REF] Chatziioannou | Energy conserving schemes for the simulation of musical instrument contact dynamics[END_REF][START_REF] Bilbao | Numerical modeling of collisions in musical instruments[END_REF] , by nonsmooth contact dynamics, which is particularly adapted for hard contacts [START_REF] Issanchou | Nonsmooth contact dynamics for the numerical simulation of collisions in musical string instruments[END_REF] or dry-friction and stick-slip transitions [START_REF] Thorin | Non-smooth dynamics for an efficient simulation of the grand piano action[END_REF] or by a Lagrange multiplier approach, particularly fit for multibody systems [START_REF] Antunes | A modal udwadia-kalaba formulation for vibro-impact modelling of continuous flexible systems with intermittent contacts[END_REF] . This former method coupled to the modal Udwadia-Kalaba (U-K) formulation seemed appropriate for the clavichord, following recent results for stringed instruments [START_REF] Antunes | Dynamical computation of constrained flexible systems using a modal udwadia-kalaba formulation: Application to musical instruments[END_REF] , including the guitar and Portuguese guitar [START_REF] Debut | Physicsbased modeling techniques of a twelve-string portuguese guitar: A non-linear time-domain computational approach for the multiple-strings/bridge/soundboard coupled dynamics[END_REF][START_REF] Debut | Physical synthesis of six-string guitar plucks using the udwadia-kalaba modal formulation[END_REF] . In Part II, a functional description of the clavichord in terms of vibrating subsystems results in a simplified one string model. Parameters are identified using experimental measurements on a copy of a historical instrument. A model of tangent/string interaction based on the Kirchhoff-Carrier string represen-tation and the U-K formulation for coupled dynamical systems is developed for the one-string clavichord model in Part III. The modal equations of the U-K model can be solved by means of a simple finite difference time discretization scheme. Synthesis results are compared to the measured dynamic behavior of the real clavichord using experiments and published data in Part IV.

II. VIBRATORY AND ACOUSTIC SUBSYSTEMS IN THE CLAVICHORD

A. Principle of the clavichord and tangent action A clavichord and its parts are described in Figure 1-(a) (instrument built in 2007 and inspired by a Hubert (1784) historical model). The main parts of the instrument are indicated on the picture. The strings, organized in pairs, are stretched between the hitch-pins and tuning pins and attached to the radiating soundboard through the bridge and bridge pin. Strings are functionally divided into three sections. The "damped section", between the hitch-pin and the tangent. This section is partially covered by strips of cloth. its vibration is rapidly damped after excitation. The "played section" vibrates between the tangent and bridge, as long as the tangent stays in contact with the string after the initial tangent strike. When the key is released, the tangent contact is lost and the string vibration is damped by the cloth strips. The "resting section" between the bridge and the tuning pin, is not directly excited by the tangent, but as it is not damped in the clavichord (contrary to e.g. the square piano), partial transmission of the played section vibration results in sympathetic vibration [START_REF] Alessandro | Tonal quality of the clavichord: The effect of sympathetic strings[END_REF][START_REF] Jiolat | L'effet acoustique des cordes mortes du clavicorde[END_REF] . The strings are pressed vertically on the bridge and pressed horizontally on the bridge pins laid out along the bridge. This contact leads to the soundboard/string coupling. The soundboard vibrates under the action of the string, and because of its large surface sound is radiated in the air.

A complete study of the instrument is out of the scope of the present work: only the simplified model displayed in Figure 1-(b) is studied in depth. This simplified model is made of only one string (the chosen string is the first string of the G3# choir) and the corresponding key, tangent and damper. The string is stretched between the hitch-pin and tuning pin and in contact with the bridge at the bridge pin. The soundboard and bridge are those of the whole instrument, because it is important to consider this part in its integrity even for modeling a unique string. Other aspects of the instrument, like the case and lid, are not considered in the present study, because their vibratory and acoustic functions are only of second order.

B. Vibratory and acoustic subsystems

For modeling purposes, the model can be considered as an assembling of four vibrating subsystems as displayed in Figure 2. The initial force for playing the instrument is provided by the player's finger. The player moves the tangent/key subsystem, the first subsystem, providing the string excitation. The strings, second subsystem, act as the vibratory engine of the instrument. Attached to the left side of the strings are the cloth dampers, third subsystem. Attached to the right side of the strings, is the bridge/soundboard subsystem, fourth subsystem.

The string is the central subsystem of the instrument. Let Y S (x, t) (resp. Y D (x, t), Y k (x, t), Y B (x B , t)) be the string displacement (resp. damper displacement, keytangent subsystem displacement, bridge displacement) at position x (resp. x, x, x B ). The vibratory subsystems are coupled to the string at points x Di , x Ltg , x B (i-th damper, tangent and bridge). The first subsystem, the tangent/key subsystem, can be considered as a rigid rod which tilts with respect to a pivot. When the tangent strikes the string, the elastic string reacts and the whole system oscillates. The tangent has a mass M T g = 5 g, the key has a mass M k = 30 g. The length of the key is L T = 28.9 cm. The pivot of the key (balance point) is situated at a distance L p = 17.2 cm of the back of the key, the finger presses the key at a distance L f = 27.9 cm of the back of the key, and the tangent is located at a distance L tg = 3.5 cm of the back of the key, see Figure 1-(b). Associated modal parameters for the tangent/key subsystem (k k , m k , c k ) are used. The string, second subsystem, is characterized by its mass, elasticity and damping factors. Between the hitch pin and excitation point, the third subsystem is the cloth strips damper, represented by N parallel dash-pots characterized by their mass and viscous damping coefficients (m D , c D ). The number of dampers is set empirically to N=65, as it seemed enough to ensure effective damping.

In the remaining of this part, parameters of the subsystems in Figure 2 are estimated using the instrument displayed in Figure 1-(a) and the G# 3 string. Three functional parts of the string can be identified: between x = 0 m and x = 0.2 m, is the damped part of the string, with a cloth damper coiled up between x = 3.4 cm and

C. The key-tangent subsystem

The key-tangent subsystem mass (determined in section III B is m k = 1.17 × 10 -2 kg. It is represesented by a rigid body mode, then its stiffness is k k = 0 N.m -1 . Its damping is set to c k = 2.5 kg.s -1 , in order to reproduce the low frequency oscillation damping of the tangentstring contact point observed in section IV B.

D. The string subsystem

The string modal damping characteristics are obtained by measurements of the vibrating string stretched on a string sand bench, allowing for measurements with almost full decoupling of the string from other vibratory structures. The Valette and Cuesta string model (see section III C 3) is used. The string is excited by a copper wire that breaks at a given tension when lifted vertically. The vibratory displacement of the string is measured at the other extremity by means of optical forks [START_REF] Carrou | A low-cost high-precision measurement method of string motion[END_REF] . For the right (sympathetic) part of the G# 3 string (L = 31.7 cm, d s = 0.33 mm, f 0 = 396.9 Hz), damping characteristics for 23 partials between 396.9 Hz to 9354 Hz are analyzed using the high-resolution algorithm ESPRIT [START_REF] Roy | Esprit-estimation of signal parameters via rotational invariance techniques[END_REF][START_REF] Carrou | Sympathetic string modes in the concert harp[END_REF][START_REF] Paté | Predicting the decay time of solid body electric guitar tones[END_REF] .

The measured damping coefficients are matched with the Valette and Cuesta model and displayed in figure 3), with the parameters reported in Table I, where E and ρ corresponds to the Young modulus and density of a brass string respectively, η air and ρ air are taken from the previous data [START_REF] Paté | Lutherie de la guitare électrique solid body: aspects mécaniques et perceptifs[END_REF] , and other parameters are defined in section III C 3. Experimental Q exp and theoretical Q th quality factors are reported together with 95 % error bars for ten measures and relative error

Q = |Qexp-Q th | Qexp .

E. Bridge and soundboard subsystems

To simulate the vibratory motion of the bridge, the modal parameters (mass matrix, stiffness matrix, damping matrix, mode shapes) of this subsystem need to be known. As no analytical solutions can be given for such a system, modal parameters of the bridge are estimated by means of experimental modal analysis [START_REF] Antunes | Dynamical computation of constrained flexible systems using a modal udwadia-kalaba formulation: Application to musical instruments[END_REF] (FRF). The FRF is obtained by measurement of the vibratory response at the coupling point between the bridge and the G# 3 string (using a PCB M352C65 accelerometer and an acquisition system with a sample rate of 51.2 kHz and a 24 bit depth). The system is excited by an automatic impact hammer (force sensor PCB 086E80). The accelerometer is placed above the bridge pin and measures its vertical acceleration. All the other strings are damped by strips of cloth woven on both sides of strings (above the soundboard on the right side and above the keyboard on the left side). Modal analysis between 100 Hz and 3500 Hz is conducted in two steps. A first step is the estimation of physical poles containing the modal frequencies and damping coefficients of the analyzed structure, using the Least square rational function (LSRF) estimation method (Matlab signal processing toolbox [START_REF] Ozdemir | Transfer function estimation in system identification toolbox via vector fitting[END_REF] ). The second step is the estimation of residues which encapsulate the mode shapes and modal masses of the system. Normalizing modal masses to m n = 1 kg for modes n = 1,2,...,N B , the corresponding mode shapes are estimated from the residues. The estimated and measured FRF at the G# 3 string/bridge coupling point are plotted on figure 4. 31 bridge modes are identified between 100-600 Hz, giving a satisfying reconstruction of the FRF. 108 bridge modes are identified between 600-3500 Hz, giving a representation of the bridge mobility for higher frequencies. High frequency modeling is necessary for a realistic simulation, even if it does not comply exactly to physical modes of the system, as previously discussed in [START_REF] Woodhouse | Plucked guitar transients: Comparison of measurements and synthesis[END_REF] . Measured and reconstructed impulse responses are given in sound example MM2.wav and MM3.wav respectively.

F. The damper subsystem, coupling between subsystems and activation

The felt damper is modeled by a series of 65 dashpots. The parameters c D , m D are chosen so that the measured damping effect exerted on the string once the key is released is well reproduced. The values chosen are m D = 1.0 × 10 -2 kg, and cD = 8.0 × 10 2 kg.s -1 . The four subsystems are coupled through the string: the string is coupled to the damper at the damper location; it is coupled to the bridge at the bridge pin location; it is coupled to the tangent/key subsystem during its contact at the tangent location. All these coupling conditions assume a continuity of displacement between the string at the coupling position and the other subsystem involved. Finally, the whole vibratory system is activated by an external force that represents the action of the finger on the key. The typical finger force profile is a step with given attack and release times.

III. PHYSICAL MODELING USING THE U-K FORMULA-TION A. U-K formulation for the clavichord

In this section, the clavichord model in terms of four coupled vibratory subsystems is modeled using a modal U-K formulation [START_REF] Antunes | A modal udwadia-kalaba formulation for vibro-impact modelling of continuous flexible systems with intermittent contacts[END_REF][START_REF] Antunes | Dynamical computation of constrained flexible systems using a modal udwadia-kalaba formulation: Application to musical instruments[END_REF] . Let us consider a mechanical system with mass matrix M which is subjected to a force vector F e (x, ẋ, t). including all constraint-independent internal and external forces. This system is also subjected to constraining forces F c (t). Denoting the dynamical solution y u (t) of the unconstrained system and y(t) of the constrained system, the motion equations of the constrained system derived by Udwadia and Kalaba [START_REF] Antunes | Dynamical computation of constrained flexible systems using a modal udwadia-kalaba formulation: Application to musical instruments[END_REF][START_REF] Arabyan | An improved formulation for constrained mechanical systems[END_REF] are:

ÿ = ÿu + M -1/2 B + (b -Aÿ u ), ÿu = M -1 F e (t) (1)
where A is the constraint matrix and b the constrained vector obtained from P h holonomic and P nh nonholonomic constraints, φ p and ψ p respectively, defined as:

φ p (y, t) = 0, p = 1, 2, . . . , P h (2)
ψ p (y, ẏ, t) = 0, p = P h + 1, P h + 2, . . . , P h + P nh (3)
Time differentiation of equation 2 and 4 gives the matrix-vector constraint equation in terms of accelerations:

A(y, ẏ, t)ÿ = b(y, ẏ, t) (4)
Note that generalized Moore-Penrose inverse matrix B + of B = AM 1/2 can be rendered numerically robust, even for a singular constraint matrix. For a particular external (finger) excitation F e (x, ẋ, t) these equations are solved using a suitable time-step integration scheme. A modal version of the U-K formulation suitable for continuous flexible systems like musical instruments is derived [START_REF] Antunes | Dynamical computation of constrained flexible systems using a modal udwadia-kalaba formulation: Application to musical instruments[END_REF] . Assuming a set of S vibrating subsystems defined in terms of their unconstrained modal basis and coupled through P kinematic constraints, one obtains:

q = W M-1 (-C q -Kq + F ext ) (5)
where q represents the vector of modal displacements, M, K, C are respectively the modal mass matrix, modal stiffness matrix, and modal damping matrix, while W = 1 -M-1/2 B + A is a global transformation matrix (which can conveniently be computed in advance of the time loop ), where A is the modal constraint matrix, and F ext are the external modal forces applied on the system.

B. Key-tangent modeling

A mode shape is associated to the key-tangent subsystem to model the tilting motion of the key. The modal representation of this system is given by :

Y k (x, t) = φ k (x)q k (t) (6)
where Y k is the displacement of the key-tangent subsystem, φ k is its mode shape and q k is its modal amplitude, and,

m k qk (t) + c k qk (t) + k k q k (t) = F ext (t) (7)
where m k , c k and k k are respectively the modal mass, modal damping and modal stiffness of the key-tangent subsystem, and F ext is the modal excitation force that the musician exerts on the key. The key is modeled as a rigid body mode in rocking motion. The mode shape φ k is given in the following way :

φ k (x) = L T -x L T -L p -1, x ∈ [0, L T ] (8)
Let the the linear density of the key be ρ k = M k L T . For any continuous linear system with density ρ, mode shape φ n and length L, the modal mass m n is :

m n = L 0 ρφ n (x) 2 dx (9)
Then, the modal mass of the key-tangent subsystem m k is the sum of the tangent modal mass and that of the key :

m k = M T g φ k (L tg ) 2 + L T 0 ρ k φ k (x) 2 dx ( 10
)
After some calculations, it gives :

m k = M T g L T -L tg L T -L p -1 2 + M k 3L 2 p -3L p L T + L 2 T 3 (L T -L p ) 2 (11)
C. String modeling

Modal description of the string

A modal representation of the string complying with the modal U-K formulation is given in this section. A modal expansion of the string displacement Y S is :

Y S (x, t) = Ns n=1 φ S n (x)q S n (t) (12)
where φ S n are the mode shapes of the string, q S n are its modal amplitudes, and N s is the number of string modes. Considering that the boundary conditions of the clavichord string are pinned-pinned, one obtains the following string's mode shapes, for a string length L :

φ S n (x) = sin nπx L n = 1, 2, . . . , N (13)
For a large enough static displacement, the geometrical non-linear force F S nl related to the string's variation of tension needs to be considered. This yields the following string's modal equations :

M S qS + C S qS + K S q S + F S nl (q S , qS ) = 0 (14)
where M S , C S and K S are the modal mass matrix, the modal damping matrix and the modal stiffness matrix of the string respectively, and q S is the modal amplitude vector of the string. In the next section geometrical nonlinear forces are expressed by means of the Kirchhoff-Carrier model.

Non-linear string dynamics in string tangent interaction

The tangent lifts the string after the initial contact and increases the string tension. This displacement can be quite significant, up to 3-5 mm, inducing a substantial rise in pitch. The string uplift is a geometrical deformation, resulting in non-linear forces that must be considered in the dynamics of the instrument. For dynamic modeling of the non-linear forces, the Kirchhoff-Carrier non-linear string model is used [START_REF] Kirchhoff | Vorlesungen Uber Mechanik[END_REF][START_REF] Carrier | On the non-linear vibration problem of the elastic string[END_REF] , following previous work on the twelve-string Portuguese guitar [START_REF] Debut | Physicsbased modeling techniques of a twelve-string portuguese guitar: A non-linear time-domain computational approach for the multiple-strings/bridge/soundboard coupled dynamics[END_REF][START_REF] Debut | Physical synthesis of six-string guitar plucks using the udwadia-kalaba modal formulation[END_REF] . Note that in this one dimensional string motion model, all couplings between transverse and longitudinal motions are neglected. This seems to be a convenient approximation for 1D simulation, even if for small string vibratory amplitude, these couplings exist and are important for detailed potential energy considerations [START_REF] David R Rowland | The potential energy density in transverse string waves depends critically on longitudinal motion[END_REF] .

According to the mode shapes in equation 13, the Kirchhoff-Carrier model leads to geometric non-linear terms for computing the dynamic tension T dyn :

T dyn (t) = ES 2L L 0 ∂Y S (x, t) ∂x 2 dx ( 15
)
which gives rise to the non-linear differential equation of motion (T 0 being the string tension at rest) :

ρS ∂ 2 Y S (x, t) ∂t 2 -(T 0 + T dyn (t)) ∂ 2 Y S (x, t) ∂x 2 = 0 (16)
The force F nl due to geometric non-linear terms is :

F nl (x, t) = T dyn (t) ∂ 2 Y S (x, t) ∂x 2 (17)
Thereby, it yields the nonlinear modal force terms :

F nl n (t) = L 0 F nl (x, t)φ n (x)dx (18)
Using equations 12 and 13 and calculating the integrals in equation 15, gives the dynamic tension that depends quadratically on the modal response amplitudes:

T dyn (t) = ESπ 2 4L 2 N n=1 n 2 (q n (t)) 2 (19)
Then, calculating the integral in equation 18, the cubic modal force terms are deduced [START_REF] Debut | Physicsbased modeling techniques of a twelve-string portuguese guitar: A non-linear time-domain computational approach for the multiple-strings/bridge/soundboard coupled dynamics[END_REF] :

F nl n = ESπ 4 8L 3 n 2 q n (t) N m=1 m 2 q m (t) 2 (20)
Equation 20 represents the modal non-linear forces for the string due to the vertical displacement resulting from the tangent lift. In contrast to the quasi-static situation, the force in equation 20 can be computed in dynamic modeling of this interaction. The increase in tension due to tangent height Y e is an important parameter for the player, as it is related to the hardness of touch 5 , i.e. the key force feedback felt by the player.

Model of string's modal dampings

The string model by Valette and Cuesta [START_REF] Cuesta | Mécanique de la corde vibrante[END_REF] is chosen to bestow a proper damping coefficient to each string mode. The air friction, the visco-elastic and thermoelastic friction as well as the structural friction are taken into account and represented by the quality factors Q n,air , Q n,ve-te and Q struc respectively.

Q -1 n = Q -1 n,air + Q -1 n,ve-te + Q -1 struc = R 2πρ L (nf 0 ) -1 + 4π 2 ρ L EIδ ve T 2 (nf 0 ) 2 + Q -1 struc (21)
where R stands for mechanical resistance

R = 2πη + 2πd s πη air ρ air f (22)
and where ρ air and η air correspond to the dynamic viscosity and the density of the air respectively, and d s represents the string's diameter. Then, Q -1 n represents the damping coefficient associated to the n th mode of the string. E is the Young modulus of the string, I is the second moment of inertia of the string, T is the string's tension, ρ L is the linear density of the string, and δ ve-te is the imaginary part of the string Young modulus. The term δ ve-te = δ ve + δ te encapsulates visco-elastic effects δ ve and thermo-elastic ones δ te , taking the same approach as [START_REF] Paté | Predicting the decay time of solid body electric guitar tones[END_REF] . Q struc is a constant value. Using

ζ n = Q -1 n D. Bridge modeling
The motion of the bridge is modeled using modal equations. Modal expansion of the bridge displacement is :

Y B (s, t) = N B n=1 φ B n (s)q B n (t) (23)
where s is the curvilinear coordinate position on the bridge, N B is the number of bridge modes, φ B n are the mode shapes of the bridge and q B n are the modal amplitudes of the bridge. Modal equations governing the bridge's vibratory motion are :

M B qB + C B qB + K B q B = 0 ( 24
)
where M B , C B and K B are the modal mass matrix, the modal damping matrix and the modal stiffness matrix of the bridge respectively, and q B is the modal amplitude vector of the bridge. Because of the complexity of the structure, as opposed to the string, analytical expression cannot be derived for bridge modal equations. Numerical values to the bridge modal parameters are obtained by experimental modal analysis (section II E).

E. Damper modeling

The damper is modeled by coupling a portion of string with a number of mass-spring-dampers assuming a continuity of displacement between the dampers and the string at their contact points. All these mass-springdamper are considered independent of one another. The modal equations governing the dampers' vibratory motion are :

M D qD + C D qD + K D q D = 0 (25)
where q D is the amplitude vector of the damper responses. The length of q D is N D , the number of dampers. Matrices M D , C D and K D are square diagonal with identical coefficients m D , c D and k D , respectively. All the mass-spring-dampers associated to the cloth damping device have the same mass, stiffness and damping coefficients. So all the mass-spring-dampers have the same frequency and the same damping. These mass-springdampers representing the cloth damper are coupled with a certain length of the string, as described in section III F 3.

F. Couplings between subsystems

The individual subsystems are described with the help of a modal representation. The modal constraint matrix A and the vector b of the constrained system is given by A Q = b with:

A =    A B A k A D    , b =    b B b k b D    , Q =      qS qB qk qD      (26)
where A B is the matrix coupling the string with the bridge with b B its associated vector, A k is the matrix coupling the string with the key-tangent subsystem with b k its associated vector, and A D is the matrix coupling the string with the damper with b D its associated vector. Given the continuity conditions for coupling, these matrices and vectors are derived in sections III F 2, III F 1 and III F 3.

String and key-tangent subsystem coupling

When the tangent touches the string, coupling between the two subsystems occurs. At the contact location, assuming a continuity of displacement between the two subsystems, the coupling conditions are :

Y S (x ext , t) -Y k (L tg , t) = 0 (27) [Φ S (x ext )] T q S (t) -Φ k (L tg )q k (t) = 0 (28)
where Y k is the displacement of the tangent, x ext is the position where the string is excited. The tangent is initially located below the string with respect to axis y. The whole string is initially at rest at altitude y = 0. At the moment when the tangent reaches altitude y = 0, A k must be modified to couple the two subsystems. The coupling conditions are :

                       If Y k (L tg , t) < Y S (x ext , t) ⇒ b k = 0 and A k = 0 If Y k (L tg , t) = Y S (x ext , t) ⇒ b k = 0 and A k = [Φ S (x ext )] T 0 . . . 0 -Φ k (L tg ) T 0 . . . 0
The U-K formulation apply constraints on the system acceleration. Respecting the constraints on acceleration does not imply respecting the constraints on the system displacement and velocity. Stabilization techniques are needed to avoid numerical drifts during the simulation, because of the displacement and velocity constraint violation. The technique used here is based on a geometric projection approach, applied after each time step t k , when small computational errors lead to violations of the holonomic ϕ(q(t k )) = 0 and non-holonomic constraints Ψ(q(t k ), q(t k )) = 0, expressed in terms of their modal coordinates and velocities. Following Yoon et al [START_REF] Yoon | Geometric elimination of constraint violations in numerical simulation of lagrangian equations[END_REF] , the displacement constraint violations are cancelled by perturbing the solution q(t k ), so that the corrected solution q c (t k ) = q(t k ) + δq(t k ) perfectly copes with the constraints ϕ(q c (t k )) = 0. This yields [START_REF] Marques | On the constraints violation forward dynamics of multibody systems[END_REF] :

ϕ(q c (t k )) = 0 ⇒ ϕ(q(t k ) + δq(t k )) = 0 ⇒ ϕ(q(t k )) + δϕ = 0 (29)
with the constraint perturbation given by:

δϕ = N n=1 ∂ϕ ∂q n δq n = Aδq ( 30
)
Where A denotes the constraint gradient matrix, in terms of the modal coordinates, which is typically nonsquare. Then, from equations 29 and 30, stems the following constraint-enforcing correction:

ϕ(q(t k )) + Aδq = 0 ⇒ δq = -A + ϕ(q(t k )) ⇒ q c (t k ) = q(t k ) -A + ϕ(q(t k )) (31)
where A + is the Moore-Penrose pseudo-inverse of matrix A. When non-holonomic constraints are applied, a procedure similar to 31 may be used for the corrected velocities qc (t k ) = q(t k ) + δ q(t k ) which perfectly cope with the velocity constraints Ψ(q c (t k ), qc (t k )) = 0. Then:

Ψ(q c (t k ), qc (t k )) = 0 ⇒ Ψ(q c (t k ), q(t k ) + δ q(t k )) = 0 ⇒ Ψ(q c (t k ), q(t k )) + δΨ = 0 ( 32
)
with the constraint perturbation given :by

δΨ = Aδ q (33)
and, from 32 and 33, stems the following constraintenforcing correction:

Ψ(q c (t k ), q(t k )) + Aδ q = 0 ⇒ δ q = -A + Ψ(q c (t k ), q(t k )) ⇒ qc (t k ) = q(t k ) -A + Ψ(q c (t k ), q(t k )) (34)

String-Bridge coupling

The string and bridge are coupled, assuming a continuity of the string and bridge displacements at the stringbridge contact point: the string displacement Y S (x B , t) must be the same as that of the bridge Y B (x B , t), where x B is the location of the coupling point on the string and x B is the location of the coupling point on the bridge :

Y S (x B , t) -Y B (x B , t) = 0 ( 35
)
with modal coordinates, it leads to :

[Φ S (x B )] T q S (t) -[Φ B (x B )] T q B (t) = 0 (36)
with the mode shape vectors :

Φ S (x B ) = [φ S 1 (x B )φ S 2 (x B ) . . . φ S N S (x B )] T , Φ B (s B ) = [φ B 1 (x B )φ B 2 (x B ) . . . φ B N B (x B )] T (37)
where N B is the number of bridge modes, N S is the number of string modes. Equation 38shows the stringbridge coupling matrix A B and the associated vector b B :

A B = [Φ S (x B )] T -[Φ B (x B )] T 0 . . . 0 (38)

String-Damper coupling

For string-Damper coupling, continuity of the string's displacement Y S (x D , t) with that of the damper Y D (r s , t) is assumed, x D being the location of the damper on the string:

Y S (x D , t) -Y D (x D , t) = 0 (39) [Φ S (x D )] T q S (t) -q D (t) = 0 (40)
Equation 40 leads to the following matrix A D and vector b D = 0, where Φ S (x Dj ) is the mode shape of the string coupled with the n th damper at the x Dj location :

A D =        [Φ S (x D1 )] T 0 . . . 0 -1 0 . . . 0 [Φ S (x D2 )] T 0 . . . 0 0 -1 . . . 0 . . . . . . . . . . . . . . . . . . . . . . . . [Φ S (x D N D )] T 0 . . . 0 0 0 . . . -1       

IV. SIMULATION AND EVALUATION OF THE CLAVI-CHORD MODEL

In this Section, assessment of the model is performed along three lines: string motion, tangent-string interaction, bridge vibration, by comparing simulated and experimental data. A finite difference approach is chosen for simulation of the one-string model developed in Section III, using the experimental data of Section II. The time step for simulation is chosen according to an energy analysis of the simulation (using N s = 150 string modes). This energy analysis studies the work done by the key-tangent subsystem compared to total energy of the coupled system for different time steps. The energy provided to the system is the work

E e = N n=1
E e,n done by the tangent, with the modal work E e,n :

T 0 F ext,n (t) qk n (t)dt (41)
The sum of all modes energy of the string, of the system tangent-key, of the damper and of the bridge gives the total energy of the modeled system. Results reported in figure 5 shows that the total energy converges towards the key-tangent subsystem work as the time step decreases.

Given the energy analysis of the simulated model, a time step of ∆t = 2 × 10-6 s has been chosen to ensures a realistic simulation. A smaller time step gives a better fit between the work done by the key-tangent subsystem and the total energy of the coupled system but it does not improve much the auditory or visual quality of the simulation, when it increases much the computational load. Assessment of the model is performed in 3 steps: 1/ visualisation of string motion; 2/ tangent-string interaction and dynamics of the clavichord; 3/ bridge vibrations and comparison with experimental data;

A. String motion

The string motion of the clavichord is shaped by the specific excitation mechanism of the instrument. Simulation of the G# 3 string motion is displayed in Figure 6, in response to a 4.2 N excitation force applied on the key. The same data are better visualized in the associated video MM1.mov, displaying evolution of the string and tangent motions in time. The top panel represents the initial 4 ms, i.e. the beginning of the motion. The tangent (represented by circles at x = 0.6 m, sampled with a period of 0.05 ms) comes in contact with the string and lifts the string to a maximum. When the tangent strikes the string, an angular point is created and propagates to the bridge. At the same time, the string is uplifted by the tangent. After the arrival of the angular point at the bridge, it is reflected back and then reflected again by the tangent. As the mechanical impedance of the bridge and that of the tangent are high compared to the string mechanical impedance, most of the wave energy is reflected. Vibratory amplitude (then the sound amplitude level) depends on the angle of the angular point, and then on the ratio of wave velocity in the string and tangent velocity, as discussed in [START_REF]On the dynamics of the clavichord: From tangent motion to sound a[END_REF] , and then on the steepness of the tangent motion slope. In the middle panel of Figure 6 the string motion history is displayed between 7 ms and 57 ms (sampled with a period of 0.5 ms). The low frequency (81 Hz) oscillation of the key-tangent subsystem because of the elasticity of the string is observed. Bottom panel of Figure 6 shows vibration of the sympathetic part of the string between the bridge and tuning pin, between 9 ms and 17 ms (sampled with a period of 0.05 ms). Note that this vibration is two orders of magnitude lower than the played part of the string, between 10 -5 -10 -6 m, and that the string motion looks rather disorganized compared to the more regular motion between the tangent and bridge pin. This Figure and associated Video1 represent the first visualization of simulated clavichord string motion to the best of our knowledge. Comparison with high-speed videos of the string motion in the vicinity of the tangent 40 shows good agreement with the simulation.

B. Tangent-string dynamics

Tangent-string dynamics is an essential feature of the clavichord dynamics, tangent velocity and displacement being the main control parameters for the performer. Figure 7 shows the velocity of the contact point between the tangent and the key. As expected 8 , the key-tangent system oscillates, because of the string elastic reaction to the tangent uplift. The simulated oscillation frequency is 81 Hz (peak velocity 1.279 m.s -1 , top panel in Figure 7), compared to the 40 Hz measured oscillation frequency (peak velocity 1.135 m.s -1 , middle panel of Figure 7). The vertical tangent velocity is measured using a B&K 4374 miniature high-sensitivity accelerometer attached to the key, close to the tangent, and a conditioning amplifier B&K 2635. With an accelerometer mass of 0.75 g and a key mass of 20 g, the mass weighting effect of the accelerometer is neglected. This frequency difference could be explained by the lack of finger weight in the model. Finger weigh would provide additional mass to the key-tangent subsystem, decreasing its oscillation frequency. In this "light key" simulation, the key-tangent subsystem modal mass is 1.17 × 10 -2 kg, with a damping coefficient of 2.5 kg.s -1 . A "heavy key" simulation (bottom panel of Figure 7), where the key-tangent subsystem modal mass becomes 2.87 × 10 -2 kg and its damping coefficient is 3.5 kg.s -1 , is consistent with measurements. Because the key is heavier in this case, the tangent impact velocity decreased (peak velocity 0.928 m.s -1 ).

Figure 8 shows the effect of the tangent uplift on string tension. Simulated tangent uplift changes the measured 41 fundamental frequency. This is consistent with measured data and quasi-static modeling results published earlier [START_REF]On the dynamics of the clavichord: From tangent motion to sound a[END_REF] . The tangent displacement controls pitch.

To study the effect of tangent velocity on acceleration at the bridge, the force F ext applied on the key-tangent subsystem is varied. The impact velocity and the average acceleration are computed by the model. Similarly to the sound pressure level (SPL) the acceleration level is computed as the logarithm of acceleration integrated over 250 ms. Figure 9 shows a linear relationship between the logarithm of the impact velocity and the acceleration level in dB. This is in good agreement with experimental results obtained for impact velocity and sound pressure level and with predictions by a linear model of string-tangent dynamics [START_REF]On the dynamics of the clavichord: From tangent motion to sound a[END_REF] . Playing fasted results in playing louder.

The influence of the impact velocity on the timbre of the bridge acceleration is studied. Spectral slopes of the acceleration spectrum for different impact velocities of the key-tangent subsystem (0.4433 m.s -1 , 0.535 m.s -1 , 0.622 m.s -1 , 0.698 m.s -1 , 0.769 m.s -1 , 0.834 m.s -1 , 0.894 m.s -1 and 0.951 m.s -1 ) are presented in Figure 10. Two effects are noticeable. First, increasing the excitation force leads to increasing the static displacement of the string, hence the fundamental frequency. This ac- counts for the frequency shift of the partials in figure 10. Second, the spectral slopes for the different spectra remain on average identical. This spectral slope of the simulated bridge acceleration is consistent with the SPL spectral slope variation with respect to the tangent impact velocity reported earlier 8 . Playing louder does not change much the clavichord's timbre.

C. Bridge vibration

Bridge vibration assessment is essential, because bridge vibration results in soundboard vibration, and then sound radiation. For assessment of the model, simulated and measured bridge motions are compared. A robotic finger [START_REF] Carrou | Dropic: A tool for the study of string instruments in playing conditions[END_REF] is used for measurement. The robotic finger presses the key following a programmed vertical motion for the string under study. All other strings are muffled using felt strips.

Measurement of acceleration at the bridge pin follows the procedure described in section II E. The first simulation takes into account all the bridge modes identified in section II E up to 3500 Hz. example MM4.wav (measured signal), sound example MM5.wav (modal truncation at 3500 Hz) and sound example MM6.wav (modal truncation at 600 Hz). Figure 12 gives a closer look to details of the three waveforms in Figure 11 at the 0.01 s, 0.1s and 0.31s respectively. The visual and sound results show good general agreement between measurement and simulation. The effect of modal truncation on the high frequency content of the simulated waveforms is clearly observed. The main difference between measurement and simulation is observed in the attack transient. This could be explained by the structural noise noticeable in a real clavichord but not taken into account in the model. Structural noise results from the shock of the tangent on the string, that excites all the case of the instrument, and coupling of the whole the string band through the damper (often coined as "drum noise"), Finally, the forces applied to the key and the response at the bridge are analyzed with the help of Figure 13. A step force of 4.2 N is applied on the key for 1 s. The constraint force at the contact point between the tangent and strings is computed. Note that the tangent force is lower than the force applied to the key, because of the leverage ratio of the pivoting key (since f -x p is smaller than x p -x tg ). Two conditions are studied in the simulation with and without the damper subsystem. The middle panel shows the simulated force at the bridge. As expected, the force is lowered during the tone, because the tangent lifts the string, and then releases the string pressure on the bridge. When the dampers subsystem is removed, the string appears less constrained, and the force lowering is higher. String vibration is apparent in the force signal. The top panel shows the vibration displacement. As expected, the string is raised in response to the tangent lift, and raised higher when the dampers are withdrawn. The vibratory magnitude is surprisingly low (a maximum of about 0.015 mm). For assessment, displacement measurements are performed on the G# 3 string using a Keyence (LJ-V7060) profilometer. The same order of magnitude is observed: a bridge lift of 0.010-0.020 mm and a maximal vibratory amplitude of about 0.010-0.015 mm, a result that is in good agreement with the simulation. After the key release, the tangent loses contact with the string. The remaining vibration after the key release corresponds to the sympathetic vibration between the bridge and tuning pin, and in the non-damped situation to the vibration of all the length of the string. In this latter situation, the magnitude is larger.

V. CONCLUSION AND PERSPECTIVES

Time-domain physical simulation of a simplified clavichord model is developed in this paper, together with measurements on a functional instrument. The clavichord is considered as a constraint system made of four vibratory coupled subsystems: key-tangent, string, damper and bridge-soundboard. In a first part, modal, dynamic and motion parameters for the four subsystems are measured on the clavichord under study. Various measurement techniques and devices are used: experimental modal analysis for the bridge, soundboard ans key (accelerometers and impact hammer), string damping characteristics (isolated string bench and optical forks), tangent and string motions (accelerometer and laser vibrometer), applied force on the key (robotic finger). These experimental results are useful for simulation assessment in this paper, as well as for reference for future studies of other instruments. The second part is devoted to simulation. The Udwadia-Kalaba (U-K) formulation is chosen because of its compatibility with a modal representation for coupled vibratory subsystems. Modal representations for the string, bridge, key-tangent and damper subsystems are developed. Another key point in this work is dynamic modeling of geometrical nonlinear forces resulting from the string uplift by the tangent, with the help of the Kirchoff-Carrier model. Conditions of continuity in displacement between each subsystem are considered to derive the coupling matrix of the constraint system. String excitation is obtained by coupling the key-tangent subsystem and the string at the moment of contact. The constraint force uplifts the string and sets string into vibration. Coupling between the string and the key-tangent subsystem requires a stabilization techniques to correct constraints violation for the displacement and the velocity, because two subsystems are not exactly in the same position at the contact moment. In the third part, simulation results are confronted to experimental results. A suitable time step for time-domain simulation is chosen according to an energy convergence analysis. The key-tangent and string dynamic of the model is compared to measurements. Low-frequency oscillation and dynamic tension variation of this subsystem are successfully rendered by the model. This is on the side of the performer's fingers. Simulated string motion is realistic, and the associated video showing velocity and displacement of the simplified model is in our opinion the peak of the present paper. Bridge motion can be considered as the source the radiated sound. Realistic force, displacement and acceleration waveforms (with associated sound examples) are obtained in response to the force applied on the key. This is on the side of the performer's ears. Results on the dynamics of the clavichord reported earlier in the literature are found in the simulation : linear relationship between the bridge log acceleration and log peak tangent velocity, pitch variation as a function of key uplift height, bridge acceleration with almost constant spectral slope for different tangent impact velocities.

A simplified 1-string instrument is simulated here. Important vibro-acoustic characteristics of the clavichord are therefore missing and must be worked out. For most instruments the tangent strikes a choir of unison strings. This influences considerably the clavichord sound, as the two unison strings are not struck exactly at the same time by the tangent. The sound of the clavichord is influenced by the vibration of the sympathetic strings, creating a reverberation effect and sometimes "wolf" notes. Sound radiation by the soundboard, reflection by the lid, structural noise must be studied. Needless to say, a full model of the instrument should include all the strings. We believe the approach presented here forms a solid bafor such a project, aiming at a fully parametric physical clavichord model that would be desirable for historical instrument simulation and analysis, performance studies and new music instrument design.

FIG. 1 .

 1 FIG. 1. (color online) (a) Description of a clavichord and (b) Sketch of the clavichord's mechanism in three positions of the touch: at rest, the tangent flushing the string and after excitation

FIG. 2 .

 2 FIG. 2. (color online) Schema of the modeled G#3 string being excited by the tangent.

5 × 10 - 4

 5104 Qstruc ηair (kg.m -1 .s -1 ) ρair (kg.m -3 ) String's parameters used to simulate the string damping coefficients x = 13.7 cm. The played part of the string is between x = 0.2 m and x = 0.53 m. The sympathetic part rest of the string is between the bridge pin x = 0.53 m and tuning pin x = 0.84 m.

FIG. 3 .

 3 FIG. 3. (color online) String quality factor: experimental (blue cross) and theoretical (blue line, 95 % error bars (blue bar) and relative error (red star).

FIG. 4 .

 4 FIG. 4. (color online) Bridge and soundboard modal analysis. Top panel: comparison of the spectral magnitude (Ref 1 dB : 1 m.s -2 .N -1 ) of the measured (blue ) and reconstructed (red ) FRF. Second panel, corresponding unwrapped phase),third panel corresponding mode shapes. Bottom panel: corresponding measured (blue) and reconstructed (red) impulse responses.

  FIG. 5. (color online) Convergence of the total mechanical energy of the system compared to the work done by the tangent.

FIG. 6 .

 6 FIG. 6. (color online) History of the transverse motion for the G#3 string. x axis: time (s); y axis : string length; z axis: transverse string motion. Top panel: full string between t = 10.0 ms and t = 14.0 ms. Middle panel: full string between t = 7.0 ms and t = 57.0 ms. Bottom panel: sympathetic part alone, between t = 9.0 ms and t = 17.0 ms.

FIG. 7 .FIG. 8 .

 78 FIG. 7. Comparison of simulated and measured velocity at the tangent-string contact point. Top panel: "light key" simulation; middle panel: measurement; bottom panel: "heavy key" simulation.

FIG. 9 .

 9 FIG. 9. (color online) Logarithmic bridge acceleration (ref 1 dB : 1 m.s -2 ) as a function of the logarithmic tangent impact velocity.

FIG. 10 .

 10 FIG. 10. (color online) Logarithmic spectral magnitude of bridge acceleration (ref 1 dB : 1 m.s -2 ), for different impact velocity of the key-tangent subsystem ), with the average spectral slope.

FIG. 11 .FIG. 12 .

 1112 FIG. 11. (color online) Bridge acceleration waveforms (1s). Top: measurement. Middle: simulation (3500 Hz modal truncation). Bottom: simulation (600 Hz modal truncation).

FIG. 13

 13 FIG.13. (color online) Bottom: force applied on the key and resulting force on the tangent. Middle: force on the bridge with (red) and without (blue) dampers. Top: bridge displacement with (red) and without (blue) dampers.

  The second simulation takes into account the first 31 identified bridge modes, leading to a modal truncation of the bridge at 600 Hz. Measured and simulated signal waveforms are displayed inFigures 11, together with corresponding sound

	acceleration [m.s-2]	-50 0 50
			0	0.2	0.4	0.6	0.8	1
	acceleration [m.s-2]	-50 0 50
			0	0.2	0.4	0.6	0.8	1
	acceleration [m.s-2]	-50 0 50	0	0.2	0.4	0.6	0.8	1
					time [s]

Preprint J. Acoust. Soc. Am. / 7 October 2021Towards a physical model of the clavichord

Preprint J. Acoust. Soc. Am. / 7 October 2021 Towards a physical model of the clavichord

, one can obtain the damping ζ n coefficients of the string.

Preprint J. Acoust. Soc. Am. / 7 October 2021Towards a physical model of the clavichord

ACKNOWLEDGEMENTS

The authors wish to thank the two anonymous reviewers and Thomas Steiner (Basel) for their criticisms and suggestions that helped much in clarifying and improving both the content and form of this work.