\section{Introduction}

Laser based diagnostics provide detailed information on plasma parameters through measurements of the phase or spectrum of light passing through the plasma.
Using diagnostics which work on fundamentally different physical principles provides a more complete picture of a given plasma, and allows results and analysis to be cross checked, validated and constrained.
In this paper we present and compare results from laser interferometry and Thomson scattering used to study a relatively simple plasma from a pulsed source.

In laser interferometry,  the phase difference between a probing beam, which passes through a plasma, and a reference beam, which passes around a plasma, can be used to infer the refractive index of the plasma.
In cold plasmas with low ionisation fractions, this refractive index can include significant contributions from both the neutral refractive index (which increases with wavelength) and the electron refractive index (which decreases with wavelength).
Two-colour interferometry is a technique which exploits these different dependencies of refractive index on probing wavelength by combining phase maps from two interferograms produced at the same time at two distinct wavelengths.\cite{Muraoka2000}
This technique is therefore used to separate the electron and neutral contribution and give separate estimates for the electron and neutral densities, and hence can provide additional constraints when analysing other diagnostics such as Thomson scattering.

Thomson scattering is a powerful tool for studying the properties of plasmas over a vast range of densities and temperatures, from cold, low temperature industrial plasmas\cite{Carbone2015} to magnetic confinement devices\cite{Peacock1969} and high energy density plasmas.\cite{Swadling2014a}
Unlike many line integrated diagnostics, such as self-emission imaging, spectroscopy, interferometry and Faraday rotation imaging, Thomson scattering provides a truly local measurement of the plasma properties inside a small volume defined by the focal spot and the collection optics.

Thomson scattering operates in two principle regimes --- collective scattering, in which the wavelength of the scattering light fluctuations is larger than the Debye length ($\lambda_{De}=\sqrt{\epsilon_0 k_B T_e/n_e e^2}$) and the light scatters from collective modes within the plasma, and non-collective scattering in which the wavelength is smaller than $\lambda_{De}$ and hence the laser scatters from individual electrons, whose motion is uncorrelated.

In the collective regime, the scattered light at a given angle can be analysed to determine the energy and momentum of these resonant collective modes, which provides information on the plasma conditions, such as the electron and ion temperatures, the electron density and the bulk velocity in the direction of the resultant scattering vector.
In the non-collective regime, the scattering spectrum is a direct measurement of the distribution function of the electrons or ions, though the noise in a realistic system limits measurements to velocities close to the thermal velocity.

Both the non-collective and collective scattering regimes can sometimes be obtained in the same experiment when the scattered light is collected from two different directions.
If the light is collected from the same plasma volumes, the parameters inferred from both regimes should be consistent.
However, in any real experimental set-up there are additional effects on the spectrum which must be carefully accounted for in order to reach agreement.
These effects may be significant in experiments which collect the scattered light from only one direction, and therefore have no inherent checks on the validity of the plasma parameters inferred from Thomson scattering.

In addition to cross checking results from Thomson scattering spectra taken at different angles, it is important to place additional constraints during fitting the spectra, as the models used are highly non-linear.\cite{Kasim2018a}
These constraints can take the form of reasonable heuristics about possible densities and temperatures, or information taken from other diagnostics.
In this paper we will demonstrate how we used two-colour interferometry to produce maps of the electron density, which can be compared with Thomson scattering fits obtained simultaneously at two different scattering angles.
Although we observe good agreement between the electron densities inferred from the interferometry and Thomson scattering, the temperatures inferred from the two Thomson scattering spectra disagree considerably, and we consider some additional broadening mechanisms which partially account for this discrepancy.

This paper is organised as follows: in \autoref{sec:experimental_setup} we describe the set up of the plasma gun and the diagnostics: laser interferometry and Thomson scattering.
In \autoref{sec:results} we present results from two-colour interferometry, which gives unphysical negative neutral atom densities in some regions of the plasma jet, but fortunately not in the region in which we perform Thomson scattering.
We then present Thomson scattering spectra, and show that we can resolve electron plasma waves in the collective and borderline non-collective regime.
In \autoref{sec:discussion} we compare the electron densities obtained from interferometry and Thomson scattering, and show there is good agreement.
We also present the temperatures measured by Thomson scattering at two different collection directions, which give wildly different results.
We propose some possible broadening mechanisms, including finite aperture effects, density gradients within the collection volume, and collisional broadening, which together account for a significant proportion of the inferred electron temperature, though not the full discrepancy measured.
Finally we look at the blue/red asymmetry in the intensity of the electron plasma waves, and discuss the possible role of heat transport and non-kinetic effects in this plasma.

\section{Experimental Setup}\label{sec:experimental_setup}

\begin{figure}[t]
	\Fig{setup}
	\centering
	\caption{Experimental set-up. a) Side view of the plasma gun, showing a cross-section of the co-axial cable and the current path. Ablation of the dielectric and copper produces a plasma jet. The Thomson scattering probe passed along a horizontal chord through the centre of the plasma jet, at \z{5}. The origin of the coordinate system is marked with a black dot. b) Top view of the plasma gun, with the jet coming out of the page. The Thomson scattered light is observed using two fibre optics arrays using a two-lens telescope. The positions of the fibre arrays defines two scattering vectors, $\textbf{k}_{s1}$ and $\textbf{k}_{s2}$. c) The resultant scattering vectors $\textbf{k}=\textbf{k}_{s}-\textbf{k}_{in}$ for fibre arrays at \deg{90} and \deg{24} to the probe beam.}
\end{figure}

\subsection{The Plasma Gun}

The plasma gun used in these experiments was originally designed for use in plasma opening switches, providing a simple, low jitter plasma source.
The gun is made from copper semi-rigid coaxial cable, with an outer diameter of 6.25 mm and a Teflon dielectric insulator (\sref{fig:setup}{a}).
The end of the cable is machined to produce an inverted cone with an angle of \deg{60},\cite{Goyer1993} and is driven using a 0.7 $\mu$F capacitor charged to 21 kV.
As the capacitor discharges, the dielectric insulator breaks down, allowing current to flow across the dielectric surface.
The current pulse rings, forming a decaying sinusoid with a peak current of 20 kA and a period of around 4 $\mu$s  --- the data presented in this paper was taken at 940 ns after current start, near the maximum current.
The surface plasma is formed by the ablation of material from the Teflon insulator and the copper conductors by the current pulse, and is therefore mostly carbon and hydrogen with some copper and fluorine impurities.

This plasma gun has previously been studied using multi-time, single colour interferometry and time-resolved single chord two-colour triature interferometry.\cite{Macdonald2015}
These diagnostics show that the plasma accelerates from the surface, reaching velocities of around 1 km/s as the plume expands into vacuum.
The electron density rapidly drops away from the nozzle of the plasma gun, initially at $n_e=$ \xcmcubed{1}{19} but reaching \xcmcubed{1}{17} after a few mm.
Spectroscopy suggests that the plasma is known to be relatively cold ($T_e=$ 5 to 10 eV) with a significant neutral gas density.

\subsection{Diagnostic Set-up}

In this paper, we extend this previous work using two diagnostics: a two-colour spatially and temporally resolved Mach-Zehnder interferometer, and spatially and temporally resolved Thomson scattering of the collective and non-collective electron feature.

We used Mach-Zehnder laser imaging interferometry to measure the line integrated electron and neutral density in the side-on ($z-x$) plane.
Interferometry was carried out using the 2nd (532 nm) and 3rd (355 nm) harmonics of a Nd-YAG laser (EKSPLA SL321P, 500 ps, 500 mJ)\cite{Swadling2014a}, with the beams imaged onto Canon 500D DSLR cameras.
The 3rd harmonic was delayed by 20 ns with respect to the 2nd harmonic, but this time-scale is short compared to the experimental time-scales and so it is reasonable to treat the two-colour probing as simultaneous.

The shutters of the DSLR cameras were left open for 1.3 seconds, allowing them to be triggered before the application of the current pulse.
The laser energy was set such that the laser intensity significantly exceeded the intensity of the self-emission from the plasma --- long exposure self-emission images (not shown in this paper) suggest that the time-integrated self-emission was three orders of magnitude smaller than the laser light.
The probing paths for the two harmonics were identical, with the beams combined using a dichroic beamsplitter, travelling through the same optics until they were split with a second dichroic beamsplitter just before the DSLR cameras.

The fringe shift in the interferograms is due to contributions to the phase from the free electrons, the neutrals and the ions.
These contributions are line integrated, and due to the azimuthal symmetry of the plasma gun we are able to use Abel inversion to estimate the phase shift as a function of radius.
The phase shift due to free electron density is linearly proportional to the laser probing wavelength, and the phase shift due to neutrals is inversely proportional to the wavelength.
By solving simultaneous equations for the phase shifts at the two different wavelengths, we can obtains maps of the electron and neutral densities.
As we will show later in this paper, this simple treatment does not always give physical results, though the reasons for these discrepancies are unclear.

Thomson scattering was used to measure the electron temperature and electron density in spatially resolved locations across the plasma jet.
A laser beam (532 nm, 2 J, 8 ns) was focused using a 750 mm achromat through the plasma, and the scattered light was imaged onto two fibre optics arrays at different scattering angles (\sref{fig:setup}{b}) by two two-lens telescopes which sit inside the vacuum chamber.
Each telescope consists of two achromat lenses (50 mm diameter, 150 mm focal length), with the first at 150 mm from the plasma, and the second lens a further 150 mm away and connected by lens tube to avoid any stray light entering the optical set up.
The second lens focuses light onto a linear fibre optic array located outside of the vacuum chamber.
The two arrays consist of seven evenly spaced 200 $\mu$m diameter fibres, which are separated by an identical ``dead'' fibre, leading to a 400$\mu$m centre-to-centre spacing which reduces ``cross-talk'' between fibres when imaged onto a CCD.
These arrays therefore collect light from seven spatially distinct volumes within the plasma.

The light from the fibres is combined onto a single 14 fibre linear array, which is imaged by a two-lens telescope onto the entrance slit of an Andor Shamrock 500i Czerny-Turner spectrometer using a grating with 1200 lines/mm.
The entrance slit width was chosen as a compromise between the signal to noise ratio and spectral resolution, with a slit width of 120 $\mu$m for the results presented in this paper, which corresponds to a Voigt-profile response function with FWHM of 2\AA{}.
The light is collected by an intensified CCD (Andor iStar 334) with a 4 ns gate, triggered by a signal from a diode which sees a fraction of the Thomson beam earlier in the amplifier chain.

Although the ion feature dominates the Thomson scattering spectra, for this plasma gun the bulk and thermal motions provide only small Doppler shifts and Doppler broadening, which cannot be resolved with this spectrometer.
Instead, we chose to saturate the ion feature on our ICCD and study the electron feature of the Thomson scattering spectra.
The electron modes are high frequency, which corresponds to a larger Doppler shift, and so can easily be distinguished from stray light at the initial laser wavelength.

Scattering can occur in the collective or non-collective regimes, determined by $\alpha=1/|\textbf{k}| \lambda_{De}$ where $\textbf{k}=\textbf{k}_s-\textbf{k}_{in}$ is the wave--vector of the mode which couples $\textbf{k}_{in}$ to $\textbf{k}_s$ and $\lambda_{De}$ is the electron Debye length.
As Thomson scattering is an almost elastic process, $|\textbf{k}_s|\approx|\textbf{k}_{in}|= k_{in}$, and so $|\textbf{k}|=k\approx2 k_{in} \sin{(\theta/2)}$, where $\theta$ is the angle between the incoming laser beam and the scattered light.
For a plasma at a specific temperature and electron density, probed by a set wavelength, there is only one free parameter left to control the value of $\alpha$ and hence determine whether the scattering is collective or non-collective, and this is the angle from which the scattered light is collected.

In the collective regime, the wavelength of the scattering fluctuations is larger than the Debye length, and hence the free electrons react collectively.
The high frequency mode is the Electron Plasma Wave (EPW) or warm longitudinal plasma wave, which follows the Bohm-Gross dispersion relationship:
\begin{eqnarray}
\omega^2&=&\omega^2_{pe}+3 k^2 v^2_{th,e}\\
&=&\omega^2_{pe}\left(1+3/\alpha^2\right) \label{eqn:eqw_dispersion}
\end{eqnarray}
where $\omega^2_{pe}=e^2 n_e/\epsilon_0 m_e$ is the electron plasma frequency.
Hence the position of the resonances, at a frequency shift of $\pm\omega$, depends on both the electron density and the electron temperature.
For a given spectrum in the collective regime, it is difficult to uniquely specify both the electron density and the electron temperature within experimental uncertainty --- this means is it is useful to have measurements from at least two angles, which allow the measurements to be constrained.
In our experiment, the spectra at \deg{24} are at large $\alpha\approx4$, where the position of the EPW resonance most strongly depends on the electron density alone, and the spectra at \deg{90} have an intermediate $\alpha\approx1$, where the EPW resonance is more sensitive to the electron temperature.
In this paper we will refer to these two fibre arrays as $k_{s1}$ and $k_{s2}$, which have scattering fluctuation vectors $k_1$ and $k_2$ respectively.

\section{Results}\label{sec:results}

All the results presented in this paper are from a single shot.
The interferometry and Thomson probing lasers simultaneously collected data at 940 ns after the start of the current pulse.
We present results from this shot because the plasma conditions provided a measurable fringe shift on both interferometers, and the Thomson scattering spectra were in a regime which allowed us to infer both the electron density and electron temperature.
During this shot series, we also gathered interferometry and Thomson scattering data from 350 ns to 1710 ns after current start, and the analysis of these results will be the subject of future work.

\subsection{Interferometry: Electron and neutral densities}

\begin{figure*}[!t]
	\Fig{interferometry_results}
	\centering
	\caption{Results from the two-colour interferometry diagnostic. The nozzle of the plasma gun is at \z{0} and between \x{\pm3}.
		a) Vacuum interferogram at 532 nm taken in the absence of any plasma.
		b) Interferogram from a 532 nm laser at 940 ns after current start. The fringes near the origin (bottom, centre) are clearly shifted compared to the fringes at the top or sides.
		c) The phase map computed by interpolation of the isophase contours in b) and subtracting the background phase. The locations of the Thomson scattering volumes are shown in cyan ($k_1$) and red ($k_2$), to scale.	d) Vacuum interferogram from a 355 nm laser in the absence of any plasma.
		e) Interferogram from a 355 nm laser at 960 ns after current start. The fringes are also shifted near \x{0}. f)  Phase map computed from e), showing significant artefacts due to large, curved fringes.}
\end{figure*}

Results from the two-colour Mach-Zehnder imaging interferometry are presented in \sref{fig:interferometry_results}.
The raw interferograms (\sref{fig:interferometry_results}{a, b and d, e}) show the distortion of the interference fringes caused by the increased electron and neutral density caused by the plasma gun.
Away from the nozzle of the gun [at $(x,z)=(0,0)$], the fringe shift decreases, and near the corners of the image, the fringes are the same as in the interferograms taken in the absence of the plasma (the ``reference'' or ``vacuum'' interferograms seen in \sref{fig:interferometry_results}{a and d}.)
To determine the line-integrated phase along the probing laser direction, it is necessary to determine the phase shift between a vacuum interferogram and a shot interferogram, and we discuss this process in more detail in Appendix \ref{sec:magic}; further examples demonstrating the fidelity of this procedure can be found in Refs. \onlinecite{Swadling2014a} and \onlinecite{Swadling2013}.

The interpolated phase maps are shown in  \sref{fig:interferometry_results}{c and f}), produced by manually tracing the constructive and destructive interference fringes using Adobe Photoshop.
The phase map for the 532 nm interferometry (\sref{fig:interferometry_results}{c})  is quite smooth, and appears to be left-right symmetric, making it suitable for Abel inversion.
The data is cut off in the bottom left corner due to the reduced fringe contrast making it difficult to reliably reconstruct the phase map.
Thomson Scattering collection volumes for arrays $k_{s1}$ (cyan) and $k_{s2}$ (red) are overlaid on this phase map, along the chord \z{5}.
The different fibre spacings are due to different magnifications for the two arrays, caused by the projection of the fibre array onto the probing laser beam from different angles.
The collection volumes for $k_1$ are elongated along the beam and collect light from an elliptical volume.

The interpolated phase map for the 355 nm interferometry (\sref{fig:interferometry_results}{f}) contains significant artefacts due to the large fringe spacing, which makes it difficult to track the exact maxima or minima of the fringes.
This leads to different interpolated phase maps between the reference and shot interferograms, and the artefacts occur during the subtraction step.
These artefacts make this phase map unsuitable for numerical Abel inversion, as the artefacts will strongly distort the final result.

In order to determine the electron and neutral density independently, we need to use the phase maps from both the 532nm and 355 nm interferometry.
This means that we need to find a way of modelling the data from the 355 nm interferometry that gives a smooth, left/right symmetric function.
We initially fit each row of pixels with a Gaussian profile multiplied by polynomials with only even powers\cite{Deutsch1983}, but found that only the lowest order polynomial (that is, a simple Gaussian profile fit) was significant, and so discarded the higher order terms.
The profiles are by construction smooth and symmetric, tending to 0 as $x\rightarrow\pm\infty$, as required for Abel inversion.

\begin{figure*}[t]
	\Fig{electron_density}
	\centering
	\caption{Electron and neutral densities obtained via Abel inversion. a) Line-out of phase along \z{5} from \sref{fig:interferometry_results}{c} (solid) and fit with a Gaussian distribution (dashed). b) Electron density obtained by analytical Abel inversion of Gaussian fits followed by simultaneously solving for electron and neutral density using phase maps from 532 nm and 355 nm interferometry. c) Neutral density - note significant region of unphysical negative neutral density around \z{2}.}
\end{figure*}

We perform the Abel inversion using the analytical result for the transform of a Gaussian profile.
The suitability of the Gaussian profile fit is demonstrated in \sref{fig:electron_density}{a}, where the dashed green fit is a good approximation to the solid green line-out along \z{5} in \sref{fig:interferometry_results}{c}.
The smoothed phase maps are first analytically Abel inverted row-by-row using
\begin{eqnarray}
A\cdot G(x,\sigma)+ C\rightarrow (A/\sigma) \cdot G(r,\sigma)
\end{eqnarray}
where $G$ is the standard Gaussian or normal distribution, A is a scaling factor and C is an additive constant which accounts for any non-zero phase shift due to vibrations in the interferometry set-up, and which vanishes during the Abel inversion.
Note the change in the horizontal coordinate from $x\rightarrow r$, which is reflected in the relabelled axes in (\sref{fig:electron_density}{b and c}).
The resulting Abel inversions are not shown, as visually they appear very similar to \sref{fig:interferometry_results}{c and f}, but the phase maps (now in units of rad/cm) are used in the standard two-colour interferometry formulas to determine the electron and neutral density:\cite{Muraoka2000}
\begin{eqnarray}
n_e&=&\frac{\phi_1\lambda_1-\phi_2\lambda_2}{r_e(\lambda_2^2-\lambda_1^2)}\\
n_a&=&\frac{\phi_1/\lambda_1-\phi_2/\lambda_2}{\alpha_a(\lambda_1^{-2}-\lambda_2^{-2})}
\end{eqnarray}

Here the subscript 1 refers to the 532 nm interferometry and 2 refers to the 355 nm interferometry, $r_e=2.83\times10^{-13}$ cm is the classical electron radius and $\alpha_a=6.4\times10^{23}$ cm$^{-3}$ is the atomic polarisability for the dominant neutral species, carbon.\cite{Schwerdtfeger2014}
The electron and neutral densities calculated through these formula are shown in \sref{fig:electron_density}{b and c}.

The electron density of the plasma jet takes the form of a broad plume (\sref{fig:electron_density}{b}).
The colour scale has been chosen to emphasise the features further away from the gun nozzle, near where the Thomson scattering was performed, at \z{5}.
As expected, the electron density drops off rapidly both in the $\pm x$ direction and in the $z$ direction.

The neutral density (\sref{fig:electron_density}{c}) has a region of unphysical negative neutral density around $(x,z)=(0,3)$ mm.
The cause of this is unclear, but it is important to note that it is not related to the value of $\alpha_a$, the static polarisability, which only determines the neutral density up to a constant.
One possible explanation is the varying atomic mix in the plasma plume --- the coaxial cable is made from a copper conductor with a plastic insulator which contains traces of fluorine.
If the ratio of these elements changes throughout the plume due to different mobilities or ablation rates, then the neutral polarisability could also change across the plume.
In particular, if there are atomic absorption resonances close to either laser probing wavelength, these would strongly alter the polarisability.

This negative neutral density has also been inferred using time-resolved continuous wave triature interferometry on this plasma gun, and has been reproduced in many previous experiments [see, for example, Ref. \onlinecite{Macdonald2015}].
These results underline the care needed when using two-colour interferometry.

Importantly, the neutral density along \z{5} is positive, and so we believe that the two-colour interferometry correctly infers the electron density in this region, and this electron density also agrees with that inferred from Thomson scattering, which will be discussed in the next section.

\subsection{Thomson Scattering: Electron densities and temperatures}

\begin{figure*}[t]
	\Fig{thomson_results}
	\centering
	\caption{Results from Thomson scattering at 940 ns after current a) Raw spectrogram showing variation of the scattering spectra in space. The top seven spectra correspond to the $k_{s1}$ scattering at \deg{24}, and the bottom correspond to $k_{s2}$ scattering at \deg{90}. Fibre 4 was aligned close to \x{0} for both fibre arrays. The change in electron density with position is clearly visible in both fibre arrays as the location of the electron plasma wave satellites are visible by eye. b) Fibre 4 from the $k_{s1}$ scattering direction, with the response of the spectrometer show in green, the ion feature (cut off) near 0 \AA{} shift and the electron plasma wave features at around $\pm50$\AA, showing scattering in the collective regime. A fit to the data is shown in red. c) Fibre 4 from the $k_{s2}$ scattering direction, response shown in green. Here the scattering is intermediate between the collective and non-collective regimes and the peaks are very broad.}
\end{figure*}

The Thomson scattering diagnostic collected scattered light from two angles in the ($x,y$) plane (\sref{fig:setup}{b}), simultaneously with the laser interferometry.
The scattered light was collected along \z{5} --- the locations of the fibres are shown in \sref{fig:electron_density}{b}.
The raw spectrogram is shown in \sref{fig:thomson_results}{a}, with the fourteen fibres visible as bright horizontal bands, and the unshifted laser light as series of bright regions down the centre of the image.
The top seven bands correspond to the seven fibres aligned in the $k_{s1}$ direction, at \deg{24} from the probing laser, and the bottom seven are from the fibres in the $k_{s2}$ direction at \deg{90}.
The spatial extents covered by the two fibre arrays are different due to the different angles --- the $k_{s1}$ array spans an extent $1/\sin{24^\circ}=2.5$ larger than the $k_{s2}$ array, and the collection volumes are elongated by the same factor in the $x$ direction.
This means that the locations of the scattering volumes do not match up for the two arrays, except for fibre 4, which was aligned close to $(x,z)=(0,5)$ mm for both fibre arrays.
In a separate shot we also collected the self-emission spectrum for the plasma with the same time after current start, exposure, camera gain etc.
The atomic lines and continuum emission were smaller by a factor of over five compared to the Thomson scattering signal shown in \sref{fig:thomson_results}{a}.

A great deal of information can be determined by simply looking at \sref{fig:thomson_results}{a}.
The most visually striking feature are the bright regions close to 532 nm, the initial laser frequency.
This includes contributions from stray light (light scattered from metal surfaces), Rayleigh scattering from neutrals (much lower cross section than electrons, but can be significant for partially ionised plasmas) and the ion feature (light scattered from electrons which Debye shield the ions, and hence determined by the ion properties).
Scattering from the electron feature occurs at larger shifts from the initial wavelength, due to the higher mobility of free electrons.
It is not possible to simultaneously resolve the electron and ion feature on our spectrometer due to the inverse relationship between resolution and bandwidth, and so we focus on the electron feature.
This means we do not need to consider contributions from stray light and Rayleigh scattering, which only occur close to the initial wavelength.

The electron feature is very clear in the $k_1$ fibres, forming a pair of evenly split satellites in each fibre.
The shift of each electron plasma wave satellite is a function of both the electron density and electron temperature, as shown in eqn.  \ref{eqn:eqw_dispersion}.
The satellites are at the largest shift in fibres 3 and 4, close to the centre of the plasma jet, and become closer to the initial wavelength as the collection volumes move further from the centre of the jet.
This corresponds to the decrease in electron density seen in the interferograms in \sref{fig:electron_density}{b}.

The electron feature can also be seen in the $k_2$ fibres, which are at a larger angle, giving a smaller value of $\alpha=1/k\lambda_{De}\approx 1$, and hence are closer to the non-collective regime $\alpha<1$.
For true non-collective scattering, we would observe a Gaussian profile, as the scattering comes from individual electrons within the Maxwellian distribution without collective effects.
Instead, we observe satellites which are far broader, but occur at larger shifts than the $k_1$ satellites due to the more significant contribution of the plasma temperature to the dispersion relationship for the $k_2$ vector (as in eqn.  \ref{eqn:eqw_dispersion}.)
The trend in shift with fibre is less clear in these fibres for two reasons: firstly, the smaller spatial range covered by this fibre array and secondly, the significant contribution of the temperature (which barely changes over the range of volumes, see next subsection) to the overall Doppler shift of the satellites.

The overall spectrum for each fibre is found by integrating in the vertical direction over the 17 pixels closest to the centre of the emission.
The spacing between fibres is 55 pixels, so this procedure avoids cross talk from other fibres contaminating the spectrum.
This unshifted light is obtained for a given slit width from a separate shot in which the spectrometer gain is set to maximum and the probe laser is used in the absence of the plasma --- in this case the spectrum comes from stray light scattered from the focal point halo by nearby metal surfaces.
The response function of the spectrometer is found from fitting the unshifted light with a Voigt profile (the convolution of a Gaussian profile with width $\sigma$ and a Lorentzian profile with width $\gamma$), resulting in $\sigma=0.6$ \AA{} and $\gamma=0.5$ \AA{}, which gives a FWHM of 2 \AA{}.
This response function is used to convolve the results from a multi-parameter fit in which specific parameters can be fixed or varied by the fitting routine.\cite{Hare2017a}
The data and model are multiplied by a ``notch filter'' to avoid fitting to the ion feature, which is not resolved in these experiments.

\sref{fig:thomson_results}{b} shows the spectrum for fibre 4 in the $k_{s1}$ array, shown with translucent blue circles.
The response function of the spectrometer is shown with a solid green line.
The red line shows the best fit, allowing the electron density, electron temperature and amplitude to vary --- the notch filter is clear here, and fitting does not occur within this region.
Note that the electron plasma wave satellites are much broader than the response function of the spectrometer, and so we are capable of resolving broadening due to electron thermal motion (temperature) or other effects.

The spectrum for the same scattering volume (fibre 4) from the $k_{s2}$ array is shown in \sref{fig:thomson_results}{c}.
The signal for this scattering angle has a lower signal to noise ratio due to the increased broadening.
For fitting purposes, we fixed the electron density for this fit to the value given from interferometry --- as the dispersion relationship is sensitive to both $n_e$ and $T_e$ for this scattering geometry, these two parameters are strongly inversely correlated in the fit and so a constraint must be added using prior knowledge from another diagnostic.
Hence for this fit, the electron temperature and the amplitude of the signal were the free parameters.

\begin{figure}[h]
	\Fig{thomson_profiles}
	\centering
	\caption{Profiles of density and electron temperature from interferometry and Thomson scattering. a) Line-out of electron density along \z{5}, averaged over $\pm$ 0.2 mm perpendicular to the line-out (solid line, orange) and density inferred from collective Thomson scattering through shift of electron plasma waves in $k_1$ array (markers, cyan). The x error bars represent the horizontal extent of the collection volume. b) Electron temperatures inferred from broadening of electron plasma wave feature for $k_1$ (cyan) and $k_2$ (red). The inferred temperatures are significantly different. Note that the volumes from which light is collected by $k_1$ and $k_2$ are not the same due to the angles at which the light is collected.}
\end{figure}

The results of these fits for each of the seven fibres in the two fibre arrays are shown in \sref{fig:thomson_profiles}.
In \sref{fig:thomson_profiles}{a}, a line-out of
 the electron density taken along \z{5\pm0.5} is shown in orange.
Overlaid is the electron density obtained from fitting the fibres from the $k_{s1}$ array, allowing $n_e$ and $T_e$ to be varied.
The horizontal error bars here indicate the spatial extent of the collection volume.
The vertical error bars are provided by the fitting code, but provide an underestimate of the true error due to the inverse correlation between $n_e$ and $T_e$, which was usually around -0.8.
The electron density inferred from interferometry is closely matched by that inferred from the collective Thomson scattering of the electron feature.

\sref{fig:thomson_profiles}{b} shows the electron temperatures inferred from the two scattering directions.
The results from the fibres in the $k_{s2}$ direction are shown in red.
To fit these fibres, the electron density was set using the mean electron density within the collection volume, as determine by interferometry (\sref{fig:thomson_profiles}{a}).
The horizontal error bars again indicate the spatial extent of the collection volume, and the vertical error bars come from the fitting routine (with no correlation to $n_e$, which is fixed).
The electron temperature profile is flat, in the range 9--11 eV across the 3 mm spanned by this fibre array.

The results from the $k_{s1}$ array are in marked contrast to those from $k_{s2}$, and are shown in cyan.
The temperature inferred is significantly higher, reaching above 30 eV in the same spatial range in which the $k_{s2}$ array infers $\sim$10 eV.
This significant discrepancy is the focus of the remainder of this paper.

\section{Discussion}\label{sec:discussion}

The agreement between the electron densities from Thomson scattering and interferometry is unexpected but reassuring, especially in light of the issues of negative neutral density seen in some regions.
More surprising, however, is the difference between the temperatures inferred from two different angles, and by a factor of up to three.
Clearly there is some additional peak broadening mechanism not accounted for by our simple analysis.
In this section we will consider three distinct mechanisms: the finite aperture effect, density gradients, and collisional broadening.
By combining these effects we can obtain significant peak broadening for the electron plasma wave satellites.

\subsection{The Finite Aperture Effect}
\begin{figure}[h]
	\Fig{thomson_broadening}
	\centering
	\caption{
		a) The scattering geometry for finite aperture effects. Note that $(x',y',z')$ are not the same as the $(x,y,z)$ in \sref{fig:setup}{}. $k_{in}$ is the probing vector at a polar angle $\theta_0'$ in the $x'-z'$ plane, $k_{su}$ is the unshifted $k$-vector directed to the centre of the aperture, and $k_s$ is the scattering vector which enters the aperture at a polar angle $\theta '$ and an azimuthal angle $\phi '$.
		b) The aperture geometry, with $\theta '$ and $\phi '$ corresponding to the scattering vector $k_s$ in a). The scattered spectrum is evaluated for $\theta '$=[0, $\theta_{max}'$] and $\phi '$=[0, $2\pi$].
		c) Simulated Thomson scattering spectra for a reference pinhole aperture with $T_e=9$ eV and $n_e=$ \xcmcubed{25}{16} (\textit{green dotted line}), the broadening due to a
		a $\mathbf{f/3}$ aperture ($\theta_{max} ' =$\deg{10}) (\textit{blue solid line}), and a fit to this broadened spectrum which infers $n_e$=\xcmcubed{22}{16} and $T_e$=17 eV (\textit{red dotted line}).
	}
\end{figure}

The equations of Thomson scattering are derived for a single value of $k$, the resultant $k$-vector.
This implies that light is collected from only a single angle, which specifies an infinitesimally small aperture or a pinhole.
In any realistic Thomson scattering system, the apertures are finite --- this is necessary to collect more scattered light in order to improve the signal to noise ratio.
In our experiment, we placed a 50 mm diameter achromat at a distance of 150 mm from the plasma, so that this first optic subtends a cone with a range of $k$-vectors from \deg{14} to \deg{34} in the plane defined by the probing beam and the centre of the lens, and $\pm$\deg{10} out of this plane.
As the central scattering vector is only at \deg{24} and the scattering spectrum changes strongly as a function of $k$, this means the lens collects scattered light with very different spectra.

In order to calculate the size of this effect, we adopt the scattering geometry shown in \sref{fig:thomson_broadening}{a and b}.
This geometry is adapted from a more complex model intended to study the effects of finite apertures on the Thomson scattering spectrum for magnetised electrons.\cite{Hare2017a}
In \sref{fig:thomson_broadening}{a}, the scattering vector which points towards the centre of the lens is labelled $k_{su}$, and lies along the $z'$ axis (which is unrelated to the $z$ axis in our experimental configuration).
The probing laser $k_{in}$ lies in the $(z',x')$ plane at an angle $\theta_0'$ to $k_{su}$ --- this is the angle of the fibre array we have been quoting in the rest of this paper.
The range of scattering vectors which can enter the lens are labelled $k_s$, and they are directed along a polar angle $\theta'$ with respect to the $z'$ axis and at an azimuthal angle $\phi'$ with respect to the $x'$ axis.

The aperture geometry is shown in \sref{fig:thomson_broadening}{b}.
The angles $\theta'$ and $\phi'$ are the same as those in \sref{fig:thomson_broadening}{a}, and $\theta_{max}'$ is the largest angle ray which can enter the lens --- \deg{10} for our system.
In order to calculate the total spectrum seen by the lens, the lens is broken up into small areas $dA=\sin{\theta'}d\theta' d\phi'$.
The spectrum is evaluated at each point on the aperture, and summed:
\begin{eqnarray}
S_{Tot}(k,\omega)=\int_0^{2\pi}\int_0^{\theta'} S(k(\theta', \phi'),\omega) \sin{\theta'}d\theta' d\phi'\label{eqn:finite_aperture}
\end{eqnarray}
For speed and clarity, the various $k$-vectors are treated as 3-vectors, and vector operations such as the dot product are used, rather than trigonometric identities.
For a typical calculation, we used eleven steps in $d\theta'$ and eleven in $d\phi'$.
The small area centred on $\theta'=0$ does not need to be calculated as $\sin \theta'=0$ here, which avoids calculating functions with singularities when magnetic fields are included.

After $S_{Tot}(k,\omega)$ is calculated, it is convolved with the response function of the spectrometer to take into account the intrinsic broadening from the instrument.
To demonstrate the effect of this broadening, we consider the parameters found for fibre 3 in the $k_{s1}$ array, in which the temperature inferred by $k_{s1}$ is 32 eV, much larger than the 9 eV inferred from the $k_{s2}$ array.
To begin with, we generate a synthetic spectrum for $n_e=$ \xcmcubed{25}{16} and $T_e$=9 eV, which is shown as a dashed green line in \sref{fig:thomson_broadening}{c}.
These parameters correspond to the average density with the collection volume of fibre 3 in the $k_{s1}$ array, and the average temperature inferred from the $k_{s2}$ array for this volume.
We note that for this collection angle the response function of the spectrometer has the same width as this convolved synthetic spectrum, and so we cannot resolve temperatures of less than 9 eV using this spectrometer and fibre geometry.

We calculate a synthetic spectrum for a finite aperture with a half-angle of \deg{10} ($f/3$) and show the spectrum of the blue shifted electron plasma wave satellite, at around -45 \AA{} to the initial wavelength (blue line, \sref{fig:thomson_broadening}{c}).
The satellite is significantly broader than the reference result (which uses an infinitesimal or pinhole aperture), and it is shifted further from the initial wavelength.
This is due to the different intensities of the electron feature as a function of $k$-vector --- for these parameters, the intensity of the scattering at larger $k$-vectors (larger angles) is stronger, and so when summed together, the spectrum has more of this character.

To test how large an effect this can have on our experimental results, we then fit the resulting broadened spectrum with our standard Thomson scattering model, which gives $n_e$=\xcmcubed{22}{16} and $T_e$=17 eV --- the inferred density has actually reduced, but this has been compensated for by a significant increase in the inferred temperature, which overall results in the peak shifting further from the initial wavelength (see eqn. \ref{eqn:eqw_dispersion}).
The significant increase in inferred temperature for the $k_{s1}$ array shows the importance of considering finite aperture effects when the collection angle is large.

\subsection{Density gradients}

\begin{figure}[h]
	\Fig{thomson_gradient_combined}
	\centering
	\caption{
		a) Simulated Thomson scattering spectra for a reference pinhole aperture with $n_e=$ \xcmcubed{25}{16} and $T_e=9$ eV (\textit{green dotted line}),
		a pinhole aperture with a range of $n_e$ taken from \sref{fig:thomson_profiles}{a} (\textit{blue solid line}), and a fit to this broadened spectrum (\textit{red dotted line}) which infers $n_e$=\xcmcubed{21}{16} and $T_e$=16 eV.
		b) The same reference spectrum (\textit{green dotted line}), a synthetic spectrum combining the finite aperture and density gradient effects(\textit{blue solid line}), and fit to the combined broadened spectrum (\textit{red dotted line}) which infers $n_e= $ \xcmcubed{21}{16} and $T_e=20$ eV.
	}
\end{figure}

A second important effect is the variation of plasma parameters throughout the collection volume, in the form of gradients in temperature and density.
As we have information on the nature of the density gradients from our interferometry diagnostic (\sref{fig:thomson_profiles}{a}), we will focus on this effect here, though temperature gradients may also play a role.
Again we consider fibre 3 in the $k_{s1}$ array, for which the interferometry shows densities between 21 and \xcmcubed{28}{16}.
Calculation of the broadening here is more straightforward --- we take the density at each point within the collection volume from the interferometry diagnostic, and calculate the spectral density function $S(k,\omega)$.
Each spectra is calculated assuming $T_e=9$ eV, using the value from the $k_{s2}$ array.
The contributions are summed and then convolved with the response function of the spectrometer, as above.

The results of this calculation are shown in \sref{fig:thomson_gradient_combined}{a}.
The model (blue solid line) is very asymmetric, clearly showing the different intensities of the scattered light from different points within the collection volume and the best fit (red dotted line) gives $n_e$=\xcmcubed{21}{16} and $T_e$=16 eV.
Again, the reduction in density from the mean density within the volume is compensated by a significant increase in temperature.
The spectrum is significantly distorted from the usual form of the Electron Plasma Wave satellite (green dotted lined), and we do not observe this distortion in our results.
However, the finite aperture effect and the density gradients work in concert, and so we must consider both together.

\subsection{Combined broadening from finite apertures and density gradients}

To combine these two effects, we calculate the finite aperture spectrum through eqn. \ref{eqn:finite_aperture} for each of the density points within the collection volume.
As before, the reference $T_e=9$ eV was used, and the densities were taken from \sref{fig:thomson_profiles}{a}.
This calculation required 7,500 spectra to be computed, each over 10,000 points in wavelength space.

The results are shown in \sref{fig:thomson_gradient_combined}{b}.
The model result (blue solid line) is significantly broader than the reference result (green dotted line).
The peak is also shifted further from the initial wavelength, and fitting (red dotted line) gives $n_e=$ \xcmcubed{21}{16} and $T_e=$ 20 eV, compared with 17 eV for the finite aperture and 16 eV for the density gradients alone.
Note that the spectrum is now quite symmetric around the electron plasma wavelength, and it is well fit by the model for a pinhole aperture at a single density.
It is difficult to distinguish experimentally between this broadened spectrum and the real, unbroadened spectrum given by a pinhole aperture in the absence of density gradients.

The inferred temperature has increased by more than a factor of two, demonstrating the importance of these two effects.
However, we still infer a temperature less than the 30 eV found by simply fitting the experimental data from the $k_{s1}$ array (\sref{fig:thomson_profiles}{b}).
This additional broadening could be due to small scale density gradients, not resolved by our interferometry, such as those caused by turbulence.
The effects of a finite aperture and density gradients were also calculated for the $k_{s2}$ array and found to be small --- these effects are far less important for this array as the both collection volume (density gradient effect) and fractional range of angles (finite aperture effect) are much smaller.

\subsection{Effect of collisions}

Collisions can also lead to broadening of the electron plasma wave features in Thomson scattering.
Using the model of Rozmus et. al (Ref. \onlinecite{Rozmus2017}), we calculate the broadening of the signal from a plasma at $T_e=9$ eV, $n_e=$ \xcmcubed{25}{16} and at an angle of \deg{24} --- the same conditions as the dashed green line in \sref{fig:thomson_broadening}{c}.
There is significant additional broadening of the unconvolved spectrum, but once the response function of our spectrometer is accounted for through convolution, the collisional broadening gives an effective $T_e=11$ eV.
This is less than the broadening from the finite aperture and density gradient effects, but does not depend on the collection optics, and so may be an important effect in other experiments.\cite{Davies2019}

\subsection{Asymmetry of the electron plasma waves}

\begin{figure}[t]
	\Fig{thomson_asymmetry}
	\centering
	\caption{Asymmetry of the blue and red shifted electron plasma waves. Data taken from fibre 4 of the $k_2$ array, with the red shifted spectrum reflected around the central wavelength so that it overlies the blue spectrum. Separately fitting symmetrised spectra gives a red shifted signal strength which is 27\% more intense than the blue shifted signal.}
\end{figure}

We note that there is a significant asymmetry in the intensities of the electron plasma waves in some collection volumes, as shown in \sref{fig:thomson_asymmetry}{} (and also visible in \sref{fig:thomson_results}{c}, the same spectrum.)
Here we have mirrored the blue shifted and the red shifted electron plasma waves around the initial wavelength and then fit them as two separate spectra to show the difference in peak intensity.
It is clear by eye that the red shifted wave is more intense than the blue shifted wave, despite the low signal to noise ratio, and the inferred intensity ratio is $I_{red}/I_{blue}\approx 1.3$.

Recent work on the Omega laser facility by Henchen et al. (Ref. \onlinecite{Henchen2018}) showed an asymmetry in electron plasma waves which was due to heat transport.
They observed measurable asymmetries when $\lambda_{ei}/L_T\sim10^{-2}$, where $\lambda_{ei}$ is the electron-ion mean free path and $L_T$ is the length scale of the temperature gradients.
In our experiments, we find $\lambda_{ei}/L_T\approx2\times10^{-3}$, which makes this plasma around an order of magnitude more collisional than that reported in Henchen et al.
Further work, including kinetic simulations, is therefore necessary to determine whether these kinetic heat transport effects could be important in the plasma gun outflow.

\section{Conclusions}

We have presented and compared temporally and spatially resolved measurements of the electron density and temperature in the outflow from a plasma gun.
Using two colour interferometry we determined the electron density in a region where the neutral density has a physical (positive) value, in contrast to other regions where the standard two-colour interferometry techniques give an unphysical (negative) neutral density.
These physical results suggest that additional sources of refractive index, such as resonances or ion polarisability, may be important in some regions of the plasma gun outflow.

We compared the results from the interferometry to those from Thomson scattering.
In the collective regime, the Thomson scattering spectrum contains resonances related to the electron plasma waves, with frequency shifts which are due to both the electron density and the electron temperature.
The electron densities measured from interferometry agreed with the measurements from Thomson scattering ($n_e\approx$ \xcmcubed{20}{17}), but the electron temperatures from Thomson scattering ($T_e\approx$ 20 to 30 eV) were much higher than previous measurements and simulations for plasma guns ($T_e\approx$ 5 to 10 eV).

We also obtained Thomson scattering spectra at a different angle, giving scattering in the non-collective regime.
The temperatures inferred from these spectra ($T_e\approx$ 9 eV) were in agreement with previous results, and in contrast to the other Thomson scattering measurements at a different scattering angle.
We attempted to reconcile these two measurements by taking a variety of broadening mechanisms into account, such as finite aperture effects, density gradients within the collection volume and collisional broadening.
Together, these mechanisms can give a inferred temperature of $T_e\approx 22$ eV, which is a significant increase of over 100\%, but is significantly less than the measured discrepancy of 200\%.
Clearly these mechanisms are important and may be significant in other Thomson scattering experiments, but the remaining broadening mechanisms are not yet understood.

Finally, we showed there was a significant asymmetry in the intensities of the electron plasma wave resonances which are red or blue shifted from the initial wavelength.
This effect has been previously reported as due to kinetic effects in electron heat transport, though the plasma in this current paper is more collisional than in the previous report.
However, the low cost and high repetition rate of a plasma gun makes it an ideal object for further investigation, which could include resolving the ion feature of the Thomson scattering or magnetising the plasma flow.

\section*{Acknowledgements}

This work was supported in part by AWE Aldermaston PLC, the Engineering and Physical Sciences Research Council (EPSRC) Grant No. EP/N013379/1, by the U.S. Department of Energy (DOE) Awards No. DE-F03-02NA00057 and No. DE-SC-0001063 and by a grant from DOE Fusion Energy Sciences user FWP100182.


\section{Introduction}
\label{sec:introduction}

\dropcap{P}rotein transport within cells relies heavily on membrane-enveloped vesicles that ferry packets of enclosed cargo  \cite{Ivanov:2008gw,Vassilieva:2008hk,Jahn:2012kq,Sudhof:2009cw}.
The content of the  vesicles is released \emph{via} their fusion with target membranes.
This transition is impeded by repulsive forces acting when the distance between the membranes is in the range of \SI{\sim 1}{\nano\meter}.
The encountered energy barrier is of the order of \SI{30}{\kbt}, implying that spontaneous fusion would take minutes, which is not fast enough in most biological situations \cite{Rand:1989iy,Leckband:2001bq,Ryham:2016gu,FrancoisMartin:2017db}.
For this reason, the process is assisted by SNARE proteins (soluble N-ethylmaleimide-sensitive factor attachment protein receptors, SNAREpins) whose conformational change (zippering) exerts forces pulling the vesicle membrane towards the target membranes.

While the total free energy change associated with the zippering process is of the order of \(  \SI{\sim70}{\kbt} \) \cite{Zhang:2017fm}, most of this energy is consumed as the SNAREpins bring the membranes into close apposition.
Biologically, the initial assembly prior to fusion provides compartmental specificity (pairing the correct SNAREs together) and
allows for temporal regulation (clamping).
Terminal zippering is then the  process which utilizes the remaining  energy  for bilayer fusion at the small (\(\SIrange{\sim1}{2}{\nano\meter} \)) separations where the repulsive forces become relevant.
Recent studies suggest that each SNAREpin can deliver only about \SI{5}{\kbt} of mechanical work at this stage \cite{Gao:2012gu,Zhang:2008bp}, which  explains why it takes about 1 second for a single SNAREpin to fuse two bilayers \cite{Xu:2016jda,Domanska:2009cw}.

It is known however that the release of neurotransmitters from synaptic vesicle occurring at nerve endings happens considerably faster, in a 10th of a \si{\milli\second} as is necessary to keep pace with action potentials, and ensure synchronous release \cite{Sudhof:2004gh,Sudhof:2009cw,Camacho:2017gd,Liu:2016cf,Li:2016isa,Krishnakumar:2015de}.
A widely accepted explanation for this remarkable difference in time scales is  that multiple SNAREpins would need to cooperate to accelerate fusion after being   synchronously released from a clamped state.
There have been indirect indications that the number of SNAREpins necessary to achieve a sub-millisecond fusion may be relatively small, ranging from 2 to 6 \cite{Hua:2001bb,Sinha:2011gxa,Mohrmann:2010eg}.
Very recently, cryo-electron microscope tomography of synaptic vesicles \emph{in situ} revealed an underlying 6-fold symmetry suggesting that exactly 6 SNAREpins are involved in such processes  \cite{Li:2019}.

How so few co-operating SNAREpins manage to accelerate fusion ten thousand times  (from \( \SI{\sim1}{\second} \) to \(\SI{\sim0.1}{\milli\second} \))   has been a mystery.
Previous modeling attempts   have suggested that more than 16  SNAREpins would be required  \cite{Mostafavi:2017gd,McDargh:2018jd}.
Here we show  that the key to understand how  only a  few  SNAREpins can achieve such rapid fusion  is the simple fact that they are mechanically coupled through effectively rigid common membranes.
The account of such mechanical coupling leads to a striking prediction  that the  number of SNAREpins must be highly constrained to ensure  sub-millisecond release  of neurotransmitters. Quite remarkably, the predicted  optimal range,     \numrange{3}{6}, is in excellent agreement with most recent experimental results \cite{Li:2019}.

We draw a fundamental analogy between the collective zippering of the SNAREpins and the  \emph{power-stroke} in a bundle of elastically coupled muscle myosin II proteins which is known to also take place at \SI{1}{\milli\second} timescale.
Building upon the seminal theory of the myosin power-stroke proposed by Huxley and Simmons \cite{Huxley_1971,Caruel:2013jw,Caruel:2016kw,Caruel:2017eba}, we model the fusion machinery as a mechanical system where the SNAREpins are represented as snap-springs interacting  through supporting membranes.
The implied  bi-stability is supported by recent experiments showing the presence of a metastable half-zipped state \cite{Gao:2012gu,Zhang:2008bp}.

The theoretical approach developed in this paper highlights  the essential role of mechanical coupling among proteins undergoing conformational changes in ensuring swift, highly  synchronized mechanical response.  This is likely a general biological principle  \cite{Caruel:2016kw,Caruel:2018jg}.

\begin{figure}[tbp]
	\centering
	\includegraphics[scale=0.22]{Model.pdf}
	\caption{
	The fusion machinery. (A) schematic of the two membranes with two attached SNAREpins. (B) Mechanical model with \( N=4 \) SNAREpins in parallel bridging the two membranes separated by the distance \( y \); two SNAREpins are in state \( c \) and two are in state \( n \) so \( N_{c}=2 \). (C) Model of a single SNAREpin. (D) Fusion energy landscape. See Table~\ref{tab:parameters} for the complete list of parameter values.
	}
	\label{fig:model}
\end{figure}

 \begin{figure*}[!htbp]
 	\centering
 	\includegraphics[scale=0.34]{Results.pdf}
 	\caption{
 	Main results. [(A) and (B)] typical stochastic trajectories of the inter-membranes distance \( y \) (A) and the number \( N_{c}(t) \) of SNAREpins in state \( c \) (B) obtained from the numerical simulation. (C) Average of the waiting times \( \tau_{1} \) (black), \( \tau_{2} \) (blue) and \( \tau_{\mbox{\tiny fusion}}=\tau_{1}+\tau_{2} \) (red) obtained from the numerical simulations (symbols) and from our effective chemical model (lines). (D) Effective free energy landscape \( \Phi \) showing the three stages of fusion and the associated transition rates. Parameters are listed in Table~\ref{tab:parameters}.
 	}
 	\label{fig:results}
 \end{figure*}

\section{Fusion Machinery}
\label{sec:fusion_machinerey}

The goal of the model is to describe the dynamic coupling between the individual SNAREpins zippering  and to study the associated evolution of the distance between the vesicle and a the target membrane.
The assembled SNARE machinery  is represented as a bundle of $N$  parallel SNAREpins bridging the two membranes, separated by the distance \( y \),  see Fig.~\ref{fig:model}[(A) and (B)].
We assume that irreversible fusion occurs when this distance reaches a critical value \( y_{f} \).
  The characteristic length associated with the deformation of the membranes generated by a zippering SNAREpin is large compared to the typical size of the SNARE bundle, see Materials and Methods \ref{sub:rigid_membrane_assumption}. Hence, the membranes can be viewed as two rigid backbones cross-linked by \( N \) identically stretched SNAREpins.

\paragraph{Single SNAREpin as a bi-stable snap-spring}The experimental work conducted in Refs.~\cite{Sutton:1998cw,Gao:2012gu,Zhang:2008bp}, suggests that a single SNAREpin can switch randomly between two metastable conformations: $n$ (half-zippered) when only the N-terminal domain of the SNAREs is zippered and $c$ (fully zippered),  when both the C-terminal and the linker domain are zippered.
To describe this process, we assume that the half-to fully zippered transition in a SNARE complex  is similar to the pre- to post-power-stroke conformational change in a myosin motor, \cite{Huxley_1971,Caruel:2016kw}.

Suppose that each SNAREpin is equipped with an internal spin type degree of freedom characterizing  the state of the protein, \( n \) or \( c \).
We denote by \( a \) the amount of shortening resulting from the \( n\to c \) transition in the absence of external load,
and by \( e_0 \) the energy difference between the two states.
This parameter can be interpreted as the typical amount of mechanical work necessary to force the \( c\to n \) transition (partial unzipping), see Fig.~\ref{fig:model}(C).

When the SNAREs are bound to the membranes, we assume that the rates $k_{+}(y)$---associated with the \( n\to c \) transition---and \( k_{-} \)---associated with the \( c\to  n\) transition---depend on the mechanical load induced by the variations of the inter-membrane distance.
To specify this dependence, both states are assumed to be  ``elastic'' in the sense  that they  exist, as phases, over an extended range of separations $y$, due to elongations of the zippered and unzippered SNARE residues, internal bonds rearrangement etc..., see Fig.~\ref{fig:model}(C).
For simplicity we assume that the deformations remains in the elastic regime so that states \( n,c \) can be associated with quadratic energies \( e_{n,c}(y) \), with minima located \( y=\{0,a\} \), and with the lumped stiffnesses \( \kappa_{n,c} \).
The transitions rates are defined so that, for a given separation, they favor the state with the lowest energy and verify detailed balance. For the detailed expressions of \( e_{n,c} \) and \( k_{\pm} \), see Materials and Methods \ref{sub:model_of_a_single_snarepin}.

\paragraph{Dynamics of the fusion machinery}
The parallel arrangement of the SNAREs implies that the conformational state of the bundle is fully characterized by $N_c$, the number of SNAREpins in state $c$.
This variable evolves according to the stochastic equation $N_c(t+dt) =N_c(t)+\{1,-1,0\}$,   with the outcomes \( \{1,-1,0\} \), characterized  by the probabilities $W_{+1}(y,N_c)=(N - N_c)\,k_+(y)dt $,  $W_{-1}(y,N_c)=N_c\,k_-(y)dt $ and $W_{0}=1-W_{+1}-W_{-1}$.
While  the SNAREpins can switch independently, the transition rates \( k_\pm \) are functions of  the  collective variable  $y$, whose dynamics in turn depends on \( N_{c} \).

 \begin{table}[tbp]
 \centering
 \caption{Physical parameters adopted in the model and references}
 \begin{tabular}{lp{0.5cm}ccp{1cm}}
 \toprule
 Parameter & Symbol & Value & Units & Reference \\
 \midrule
 Zipping distance & $a$ & \num{7} & \si{\nano\metre} & \cite{Gao:2012gu}\\
 Energy bias & $e_0$ & 28 & \si{\kbt} & \cite{Gao:2012gu}\\
 Fully zipped stiffness & $\kappa_{c}$ & \num{12} & \si{\pico\newton\per\nano\metre} & SI** \\
 Half-zipped stiffness & $\kappa_{n}$ & \num{2.5} & \si{\pico\newton\per\nano\metre} & SI  \\
 Maximum zippering rate & $k$ & \num{1} & \si{\mega\hertz} & \cite{Gao:2012gu}\\
 \midrule
 Drag coefficient & $\eta$ & \num{3.8e-07} & \si{\newton\second\per\meter} &  \\
 \midrule
 FB* position & $y_f$ & \num{2} & \si{\nano\metre} & \cite{Evans:1991fz, Rand:1989iy}\\
 FB width & $\sigma_f$ & \num{0.3} & \si{\nano\metre} & \cite{Leckband:2001bq,Donaldson:2011cn} \\
 FB height & $e_f$ & \num{26} & \si{\kbt} & \cite{Ryham:2016gu,FrancoisMartin:2017db}\\
 \bottomrule
 \label{tab:parameters}
 \vspace{-0.5cm}
 \end{tabular}
 \flushleft \vspace{-0.3cm}
 \addtabletext{*FB: Fusion Barrier; **SI: SI Appendix; \( \SI{1}{\kbt}\approx\SI{4}{\zepto\joule} \)
 }
 \end{table}

To specify the coupling between the two degrees of freedom \( N_{c} \) and \( y \), and thereby formulate the complete model of the fusion process, we first recall that the  motion of the vesicle in the overdamped regime results from the balance between the force applied by the \( N \) SNAREpins, the membrane repulsion and the viscous drag.
Taking into account the thermal fluctuations, this force balance translates into the stochastic equation,
\begin{equation}
	\label{eq:dynamics_y}
		\eta \dot y  =    -\frac{\partial}{\partial y} (E_{\mbox{\tiny snare}}+E_{\mbox{\tiny fusion}})  + \sqrt{2\eta\,k_{\mbox{\tiny B}} T}~\xi(t),
\end{equation}
where $\xi(t)$ is a standard  white noise, and  $\eta$ is a drag  coefficient  representing the friction opposing the motion of the vesicle.
At a given \( y \), the force applied by the bundle derives from the sum of individual SNAREpin energies $E_{\mbox{\tiny snare}}(y,N_c) = N_c\, e_c(y) + (N-N_c)\, e_n(y)$.
Finally, the inter-membrane repulsion, due to short-range   forces  between the two membranes, is schematically modeled by a Gaussian energy barrier \cite{Rand:1989iy} $E_{\mbox{\tiny fusion}}(y)= e_f \exp[-(y-y_f)^2/(2\sigma_f^2)]$, where $y_f$ is the critical separation,  while $e_f$  and $\sigma_f$ are the height and the width of the barrier, respectively, see Fig.~\ref{fig:model}(D).
  The ensuing dynamics of the  system unfolds in the space of two variables:  the continuous  one, $y(t)$, and the integer-valued one, $N_c(t)$.
  The associated energy landscape has a multiwell structure that accounts for the configurational states of \( N \) individuals.
  The response is governed by the two stochastic equations, for which initials conditions still need to be specified.
 We consider the initial state \( N_{c}(t=0)=0 \), and  \( y(t=0) = a \), which corresponds to the configuration where the SNAREpins are at the bottom of the energy well describing  state $n$.
 This configuration characterizes the system immediately after the calcium induced collapse of Synaptotagmin which triggers the full zippering of the SNAREs \cite{Sudhof:2013fw}. This point is discussed in more details in Section~\ref{sec:discussion}.

\paragraph{Model parameters}
The model is calibrated as follows, see Table~\ref{tab:parameters} and SI Appendix~(A) for further details. The mechanical parameters characterizing a single SNAREpin, \( a \), \( e_0 \), \( \kappa_{n} \) and \( \kappa_{c} \), are determined by using our model to reproduce the experimental results obtained from stretching tests with optical tweezers, see Refs.~\cite{Gao:2012gu,Zhang:2017fm}.
The energy bias \( a \) and \( e_0 \) are chosen to be compatible with the results obtained from these studies.
The procedure used to estimate the stiffnesses \( \kappa_{n,c} \) is more complex and explained in details in SI Appendix.
The value of the rate \( k \), is fixed in accordance with estimates from Refs.~\cite{Gao:2012gu,Yang:2003ib}.
The drag coefficient is computed using the Stokes formula $\eta = 6\pi \mu R$, where $R=\SI{20}{\nano\meter}$ is the vesicle radius, and $\mu=\SI{e-03}{\pascal\second}$ the fluid viscosity.
 The corresponding characteristic timescale is $\tau_{\eta} = \eta a^2 / (k_{\mbox{\tiny B}} T) \approx \SI{4.5}{\micro\second}$.

The values of the parameters \( y_{f} \), \( \sigma_{f} \) and \( e_f \) are chosen to be compatible with the current literature \cite{Ryham:2016gu,Lentz:2009jv, Zhang:2008bp,Cohen:2004ds,Mostafavi:2017gd,McDargh:2018jd,FrancoisMartin:2017db,markvoort2011,finkelstein1986,grafmuller2007,grafmuller2009,lee2008,Xu:2016jda}.
In particular, values of \( e_f \) between 26 and \SI{34}{\kbt} have been reported for various types of lipids. We chose \SI{26}{\kbt} [POPC\footnote{1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine lipid}, see Ref.~\cite{FrancoisMartin:2017db}] which leads to a single SNAREpin average fusion time of \SI{1}{\second}.

\section{Results}
\label{sec:results}

 \paragraph{Numerical simulations} Typical stochastic trajectories, \( y(t) \) and  \( N_{c}(t) \) obtained from numerical simulations, are shown in  Fig.~\ref{fig:results}[(A) and (B)]. They indicate that the fusion process can be decomposed into two stages, characterized by the times \( \tau_{1} \) and \( \tau_{2} \).
 During the first stage, the system remains in its initial configuration \( (y\simeq a, N_{c}=0) \) with only isolated  \( n\to c\to n \) transitions. After a   time \(  \tau_{1} \), the   inter-membrane distance drops abruptly to \( y\simeq \SI{2.5}{\nano\metre}\), while all the SNAREpins collectively switch from state \( n \) to state \( c \).
 The inserts in Fig.~\ref{fig:results}[(A) and (B)] show that this transition occurs within \SI{10}{\nano\second} after the inter-membrane distance has reached the value \( y=y_* \); the irreversible collective zippering  itself (\( N_{c}=0\to N_{c}=4 \)) lasting about \SI{1}{\nano\second}.
After the synchronized \( n\to c \) transition the inter-membrane distance remains above the threshold \(y= y_f \) for a time \( \tau_2 \) before fusion. The duration of the whole process is therefore \( \tau_{\mbox{\tiny fusion}}=\tau_{1}+\tau_{2} \).

The mean timescales \( \bar{\tau}_{1} \) and \( \bar{\tau}_{2} \) (obtained by averaging \num{e03} stochastic trajectories) are represented as functions of the number of SNAREpins in Fig.~\ref{fig:results}(C) on a semi-logarithmic scale. Observe  that  \( \bar{\tau}_{1} \) increases exponentially with \( N \) and   \( \bar{\tau}_{2} \) decreases exponentially with \( N \).
These antagonistic  \( N \)-dependencies result in  the average fusion time \( \bar \tau_{\mbox{\tiny fusion}}=\bar\tau_{1}+\bar\tau_{2} \) exhibiting a remarkably sharp minimum, see Fig.~\ref{fig:results}[(C), red].
With the set of parameters values reported in Table~\ref{tab:parameters}, this minimum is attained at \( N_{*}=4 \) and is associated with a fusion timescale of \SI{\sim 100}{\micro\second}. In addition we obtain a fusion time of the order of \SI{1}{\second} for a single SNAREpin.
Both values are consistent with \emph{in vitro}  \cite{Xu:2016jda} and \emph{in vivo} \cite{Meinrenken:2002dc,Sudhof:2004gh} experimental measurements.

\paragraph{Fusion as a two-stage  reaction} To elucidate the mechanism of fusion in two stages, we present here a `toy' model where the whole process is recast as two successive reactions
\begin{equation}
\label{eq:chemical_reaction}
\ce{unfused <=>>[$k_{1}$][$k_{-1}$] IS ->[$k_{2}$] fused},
\end{equation}
where IS stands for an intermediate state whose characteristics depend on the mechanical properties of the zippered SNAREpins. In this representation,  the fusion is viewed as the outcome  of two distinct sub-steps: the collective zippering and the topological membrane merger.

To justify such model reduction, we assume that the timescale of the  \( n\rightleftharpoons c  \) transition is negligible compared to the timescale describing the relaxation of the vesicle position.
In the corresponding limit ($kt_{\eta}\gg 1 $) \eqref{eq:dynamics_y} can be averaged with respect to the equilibrium distribution of the variable \( N_{c}(t) \) [see Materials and Methods~\ref{sub:adiabatic_elimination_of_the_variable_boldsymbol_n__c}], so the original system reduces to the one-dimensional stochastic equation
\begin{equation}
	\label{eq:effective_langevin}
	\eta\dot{y} = -\,\frac{\rm d}{\rm{d} y} \big[N f_{\mbox{\tiny snare}}(y)  + E_{\mbox{\tiny fusion}}(y) \big]+ \sqrt{2\eta\,k_{\mbox{\tiny B}} T}\,\xi(t),
\end{equation}
where, the energy \( E_{\mbox{\tiny snare}}(N_{c},y) \) appearing in \eqref{eq:dynamics_y}---which depends on \( N_{c} \) and \( y \)---is replaced by the equilibrium free energy $f_{\mbox{\tiny snare}}(y)=-k_{\mbox{\tiny B}} T \log\,\{
	\exp\,[-e_c(y)/(k_{\mbox{\tiny B}} T)\,] + \exp\,[-e_n(y)/(k_{\mbox{\tiny B}} T)\,] \,\} $---which depends only on  \( y \). This free energy is illustrated in Fig.~\ref{fig:results}[(D) dashed line].
The overall potential $\Phi(y)= N f_{\mbox{\tiny snare}}(y)+E_{\mbox{\tiny fusion}}(y)$, driving the effective dynamics (\ref{eq:effective_langevin}) is shown by the solid red line in Fig.~\ref{fig:results}(D). It exhibits two local minima representing two metastable states.
The first metastable state [point $A$ in Fig.~\ref{fig:results}(D)] is located at \( y\simeq a \), where on average all the SNAREpins are in state $n$.
The second one  [point $B$ in Fig.~\ref{fig:results}(D)] is located at \( y_{f}<y_{2}<y_{*} \) and represents the  \emph{intermediate} state where on average all the SNAREpins are in state $c$, still confronting a reduced fusion barrier.

The system evolving   in this energy landscape from the initial---\( y\simeq a \)---to the final---\( y=y_f \)--- state faces two successive  energy barriers $\Delta\Phi_1$ and $\Delta\Phi_2$.  With each barrier $\Delta\Phi_{1,2}$ one can associate a waiting time $\bar\tau_{1,2}$, that can be approximated by the Kramers formula \cite{Kramers_1940, Risken_1988,Schuss_2010}
\begin{equation}
	\label{eq:kramers}
	\bar\tau_{1,2} = \tau_{\eta}\alpha_{1,2} \exp\left[\Delta\Phi_{1,2}/(k_{\mbox{\tiny B}} T)\right].
\end{equation}
The values of the numerical prefactors \( \alpha_{1,2} \) are determined by the local curvatures of the potential \( \Phi \) at its critical points and depend weakly on $N$, see Materials and Methods~\ref{sub:adiabatic_elimination_of_the_variable_boldsymbol_n__c}.

The approximated timescales $\bar\tau_{1,2}$  are compared with the  numerically computed values $\tau_{1,2}$  in Fig.~\ref{fig:results}[(C) solid lines].
The excellent agreement between the two sets of results suggests that the whole fusion process can effectively be described by two successive `chemo-mechanical' reactions, and that the rates in \eqref{eq:chemical_reaction} can be computed from the formulas \( k_{1,2} = \bar\tau_{1,2}^{-1} \)   while the remaining rate \( k_{-1} \) is prescribed by the condition of detailed balance.

Finally, note that with the parameters reported in Table~\ref{tab:parameters}, \( k\tau_{\eta}=4.5 \), which shows that our effective model is accurate even if the condition \( k\tau_{\eta}\gg 1 \) is not fully satisfied.

The peculiar dependencies of the waiting times $\tau_{1,2}$ on the number of SNAREpins $N$ can be now  understood  by  referring to  the  $N$-dependence of the energy barriers \( \Delta\Phi_{1,2} \).

\paragraph{ Timescale \( \boldsymbol{\tau_{2}} \): The cooperative action of the SNAREs reduces the time for crossing the fusion barrier.}

In the intermediate state [point $B$ in  Fig.~\ref{fig:results}(D)], the SNAREpins are all in  state $c$ and the pulling force they apply on the membranes is exactly balanced by the short range repulsive forces.
The system remains trapped in this state until a thermal fluctuation delivers the energy \( \Delta\Phi_{2} \) allowing the system to reach the distance \( y=y_{f} \), where the fusion occurs.

In the absence of SNAREs, this energy difference is simply the bare fusion barrier $e_f$, see Fig.~\ref{fig:model}(D).
When the SNAREpins are present, the total force they apply brings the two membranes in close contact, which reduces the energy barrier.
This effect is amplified by an increase in  the number of SNAREpins: the larger  the number of SNAREpins, the larger the overall force, so the closer the membrane can be brought together, see Fig.~\ref{fig:results}(D).

Since the inter-membrane potential  \(E_{\mbox{\tiny fusion}}(y)\) decays rapidly as \( y \) increases, we can approximate the second energy barrier by \(\Delta\Phi_{2}\simeq e_f-Nw \), where \( w \) represents the amount of mechanical work that a single SNAREpin can deliver, see Materials and Methods~\ref{sub:estimation_of_the_mechanical_work_w} for the derivation of this result and the mathematical expression of \(w\).
According to \eqref{eq:kramers}, we then have
\begin{equation}
	\label{eq:tau_two}
	\bar\tau_{2} (N) \propto \exp\left[-Nw/(k_{\mbox{\tiny B}} T)\right],
\end{equation}
hence the exponential decay of the time \( \tau_{2} \) with the number of SNAREpins.

With the parameters of Table~\ref{tab:parameters}, each SNAREpin provides a mechanical work \( w\simeq \SI{4.5}{\kbt}\) when it encounters the fusion barrier, which reduces the average time for fusion \( \bar\tau_2\) by a factor of \(\sim 100\), see Fig.~\ref{fig:results}(D).
This multiplicative effect allows fast fusion at the sub-millisecond timescale with as few as three SNAREpins.
For large enough number of SNAREpins (here \(N>7\)), the overall applied force  surpasses the membrane repulsion and the remaining fusion barrier disappears.
The obtained exponential decay of the timescale \( \tau_{2} \) with the number of SNAREpins suggest that the fusion could in principle proceed much faster than \( \SI{\sim 100}{\micro\second} \), being only limited by viscous forces.
Considering that each vesicle can accommodate up to \( \num{\sim100} \) SNAREpins, one cannot rule out the possibility of neurotransmitter release occurring much faster than \SI{100}{\micro\second}.
 Next we argue that such scenario is unlikely by showing that the fusion process gets slowed down if the number of SNAREpins becomes too large.

\paragraph{Timescale \( \boldsymbol{\tau_{1}} \): Increasing the number of SNAREpins slows down the synchronous zippering.}

The average time \(\bar\tau_1\), taken for all the SNAREpins to switch from the \( n\) to the \(c\) conformation and then pull the membranes toward the bottom of the fusion barrier, exponentially increases with the number of SNAREpins \(N\), see Fig.~\ref{fig:results}C.
This dependence can be explained as follows.

As long as \( y_{*}<y<a\), the individual transition rates are such that \( k_+<k_- \) which implies that, on average, all the SNAREpins are in state \( n \) and therefore under compression, see Fig.~\ref{fig:model}(C). This idea is in agreement with the experimental results from Ref.~\cite{Zhang:2017fm}, that revealed the presence of the half-zipped metastable state.
Consequently in the interval \( y_{*}<y<a\) the average force \(- \rm{d}f_{\mbox{\tiny snare}}/\rm{d}y\), that the SNAREpins collectively exert on the membranes, is repulsive.
Beyond the point \( y\simeq y_{*} \), the state \( c \) is stabilized (\( k_{+}>k_- \)) and the average force becomes attractive.
Since this force is proportional to the number of SNAREpins, the waiting time before a fluctuation can provide enough energy to surpass the repulsion---and overcome the barrier \( \Delta\Phi_{1} \)---increases with \( N \).
This constraint results from the mechanical feedback
induced by the membranes. The latter  play the role of a rigid backbone that forces  the SNAREpins to bridge approximately the same intermembrane distance, see Refs.~\cite{Caruel:2015im,Caruel:2017eba,Caruel:2018jg}.

To specify the \( N \)-dependence of \( \tau_{1} \), we use the fact that for \( y>y_{*} \) we can consider that \( E_{\mbox{\tiny fusion}}=0 \), so \( \Delta\Phi_{1} \) can be approximated by \(\Delta\Phi_{1}\simeq N[f_{\mbox{\tiny snare}}(y_{*})-f_{\mbox{\tiny snare}}(a)] \simeq  \Delta e-k_{\mbox{\tiny B}} T\log(2) \), where \( \Delta e=e_{n}(y_{*})-e_0\). According to \eqref{eq:kramers}, we can then write
\begin{equation}\label{eq:tau_one}
\bar\tau_{1}(N)\propto \exp\left[N\Delta e/(k_{\mbox{\tiny B}} T)\right],
\end{equation}
which shows that the timescale \( \bar{\tau_{1}} \) increases exponentially with the number of SNAREpins.

 \begin{figure}[tbp]
 	\centering
 	\includegraphics[scale=0.4]{EffectOfDeltaV.pdf}
 	\caption{
	Effect of the intrinsic energy barrier \( \Delta e \) on the optimal number of SNAREpins (A) and on the associated fusion time (B).The parameters values are taken from Table~\ref{tab:parameters} with \( \kappa_{n} = \SIrange[range-phrase = -,range-units = single]{0.11}{24}{\pico\newton\per\nano\meter}\).
 	}
 	\label{fig:EffectOfDeltaV}
 \end{figure}

From \eqref{eq:tau_one} we obtained that the first energy barrier is fully controlled by a single parameter  \( \Delta e \) which therefore has a strong influence on both the existence and value of the optimal number of SNAREpins.
To study the  effect  of $\Delta e$ on the fusion time, we varied the parameter \( \kappa_{n} \)  describing  the curvature of the energy \( e_{n} \). The results of our parametric study are summarized in Fig.~\ref{fig:EffectOfDeltaV}.
These data were obtained by using \eqref{eq:kramers} to compute the intersection of the curves \( \bar\tau_{1,2}(N) \) for each value of \( \Delta e \).
We checked that the  results are in good agreement with  direct numerical simulations. Despite the broadness of the interval of parameter values tested, the optimal number of SNAREpins remains below \num{10}.
If we consider only the cases corresponding to sub-millisecond fusion times, we obtain \( N\geq 3 \) with \( \Delta e<\SI{4}{\kbt} \).
The latter value is compatible with recent estimate of \( \Delta e\approx \SI{5}{\kbt} \) for the \( n\to c \) transition energy barrier, see Refs.~\cite{Li:2016isa,Zhang:2017fm}.
Note also that the predicted optimal  number of SNAREpins   is robust  because it corresponds to a  \emph{plateau}  on the  \( N_{*}(\Delta e )\) curve, see Fig.~\ref{fig:EffectOfDeltaV}(A).

\paragraph{Robustness of the predictions. }

\begin{figure}
	\centering
	\includegraphics[scale=0.4]{Robustness.pdf}
	\caption{
	Robustness of the prediction. Influence of the parameters \( e_0 \) (A), \( \kappa_{c} \) (B), \( e_f \) (C) and \( \sigma_f \) (D), on the timescales \( \bar{\tau}_{1,2} \) and \( \bar \tau_{\mbox{\tiny fusion}} \).
	The results were obtained using \eqref{eq:kramers}.
	}
	\label{fig:Robustness}
\end{figure}

 The results presented above, were obtained for the  parameters values listed in Table~\ref{tab:parameters}. For some of these parameters only a rough estimate  is available  at this stage, see SI Appendix.
To test the robustness of our theoretical predictions, we computed the average waiting times $\bar\tau_{1,2}(N)$ from \eqref{eq:kramers} for different values  of four key parameters of the model, $e_0$, $\kappa_{c}$, $e_f$, $\sigma_f$, see Fig.~\ref{fig:Robustness}.
For each of these parameters, the lower and the upper bounds delimit broad intervals covering the values obtained from different experimental studies.

A comparison between Fig.~\ref{fig:Robustness} and Fig.~\ref{fig:results}(B) shows that our  results are  only marginally affected by changes in the parameter values. In particular, the existence of a sharp minimum of the fusion time, associated with an optimal number of SNAREpins, is a robust prediction.
In addition, the value of the optimal number of SNAREpins, is weakly sensitive to the parameters: it always remains in between $3$ and $6$.
Remarkably, despite the large difference between the upper and lower bounds for each of the parameters, the average fusion time remains in the sub-millisecond scale.

The energy landscape associated with the zippering of the SNARE complexes is the object of intense current research   \cite{Li:2016isa,Walter:2010dt,Zhang:2017fm}. In our model this landscape is fully characterized by only four parameters: the distance \( a \), the energy  \( e_0 \), and the stiffnesses \( \kappa_{n,c} \).
While the distance \( a \) has been measured with precision in recent works \cite{Gao:2012gu,Zhang:2017fm}, the values of the other three parameters are still not known with certainty.  Several estimates of the energy bias \( e_0 \) lying between \num{20} and \SI{40}{\kbt} can been found in the literature \cite{Zhang:2017fm,Li:2007ka}. We show in Fig.~\ref{fig:Robustness}(A) that variations within this interval  affect  mostly \( \tau_{1} \), change the fusion time by one order of magnitude,  but have  almost no effect on the optimal number of SNAREpins.
Currently, only indirect evaluation of the stiffnesses \( \kappa_{n,c} \) can be obtained from the available data, see SI Appendix.  Within the broad range  of values tested in our numerical simulations we again observed only small variations of the optimal number of SNAREpins, see Fig.~\ref{fig:EffectOfDeltaV}, and Fig.~\ref{fig:Robustness}(B).

One of the most documented physical phenomena involved in the fusion process is the merging of the two membranes. The amplitude of the associated repulsion force  depends, in our model, on the parameters \( e_f \) and \( \sigma_f \), whose influence  on the fusion time is illustrated in Fig.~\ref{fig:Robustness}[(C) and (D)], respectively.
As expected from the analysis presented in Section~\ref{sec:results}, changing the values of these two parameters   affect only the height of energy barrier \( \Delta\Phi_{2} \) and therefore the timescale \( \tau_{2} \).
Increasing $e_f$ raises the height of the maximum of $E_{\mbox{\tiny fusion}}$ [see Fig.~\ref{fig:model}(D]), while decreasing $\sigma_f$ deepens the second energy well [point B in Fig.~\ref{fig:results}(D)], which results in both cases in the increase of $\tau_2$.
This leads \emph{in fine} to the increase of the optimal number of SNAREs.
Notice that $e_f$ depends on the type of lipids, on the membrane curvature and is also strongly sensitive to the membrane tension \cite{markvoort2011,finkelstein1986,grafmuller2007,grafmuller2009,lee2008}. Therefore, its value can be different in different cells or experimental set-ups. In particular, we expect the \emph{in vivo} value to be smaller than the value measured in artificial systems (\SI{35}{\kbt} ), which in general use low tension and low curvature membranes, see Ref.~\cite{FrancoisMartin:2017db}.

In conclusion, while additional experimental studies are needed to refine the calibration of the  model, the above parametric study shows the robustness of the effects of the mechanical crosstalk between the SNAREpins.

\section{Discussion}
\label{sec:discussion}

In this paper, we elucidated the central role played by mechanical coupling in synchronizing the activity of  SNAREpins, which is necessary to enable sub-millisecond release of neurotransmitters.
Our approach to the problem   complements previous studies focussed predominantly on the molecular details of the single SNARE zippering transition \cite{markvoort2011,risselada2012,Walter:2010dt,Acuna:2014km,Hernandez:2012cg,Li:2014dg,Sutton:1998cw,Stein:2009jy, Gao:2012gu,Zhang:2017fm}.

As a starting point we used  a previously unnoticed analogy between the activity of SNARE complexes  and the functioning of myosin II molecular motors.
Viewed broadly, both systems ensure ultra-fast mechanical contraction.
In the case of muscle, destabilization of the  pre-power stroke state is the result of a mechanical bias created by an abrupt shortening of the myofibril \cite{Huxley_1971, Ford_1977}. In the case of SNAREs, similarly abrupt destabilization is a result of the calcium-induced removal of the synaptotagmin-based clamp, most likely when \ce{Ca^2+} triggers disassembly of the synaptotagmin ring \cite{Ramakrishnan:2019}.

To pursue this analogy, we developed a variant of the power-stroke model of  Huxley and Simmons ~\cite{Huxley_1971,Caruel:2016kw}, in which the zipping is viewed as a transition between two discrete states endowed with different elastic properties.
This representation is supported by recent experiments \cite{Gao:2012gu,Zhang:2017fm}, which provided essential data for the calibration of the model.

Our analysis of the collective behavior of $N$ ``switchers'' of this type suggests that the main function of the SNARE machinery is to bring the two membranes to a distance beyond which the fusion process can proceed spontaneously.
The emerging intermediate configuration, where the two membranes are sufficiently closely tethered, can be then viewed as an intermediate state in the reaction process linking the fused and unfused states.
The result is a representation of the SNARE mediated fusion as a two stage reaction.

We linked the first stage of the process with the collective zippering of the SNAREpins and showed that this step gets exponentially more sluggish as the number of SNAREpins increases. This phenomenon was studied previously in the context of muscles, see Refs.~\cite{Caruel:2013jw,Caruel:2018jg}. It originates  (i) from the experimentally suggested presence of a metastable half-zipped state along the zippering free energy landscape \cite{Zhang:2017fm}, and (ii) from the long-range mechanical interactions mediated by the scaffolding membranes, which create a negative feedback preventing a fast collective escape from the metastable half-zippered state.

The second stage of the process is the transition from the intermediate state to the fused state. The associated timescale \( \tau_{2} \) decreases exponentially  with the number of SNAREpins because the larger the number of acting SNAREpins, the closer the membranes can be brought together in the intermediate state and therefore the higher is the energy of this state. This results in an exponential decay of the timescale \( \tau_{2} \) with the number of SNAREpins.
Behind this phenomenon is the presence of a residual force in the configuration where the SNAREpins have reached the intermediate state. This perspective is supported by the results of Ref.~\cite{Gao:2012gu,Zhang:2008bp}.

The antagonistic \( N \)-dependence of the rates characterizing the two stages reveals the existence of an optimal number of SNAREs \( N_{*} \) allowing the system to perform fusion at the physiologically appropriate timescales. Our prediction \( N_{*}=\numrange{4}{6} \) is supported by recent in situ cryoelectron microscope tomography observation, see Ref.~\cite{Li:2019}.

We remark that our result strongly depend on the initial configuration of the system, which we link with the structure of the fusion machinery immediately after Synaptotagmin removal by calcium.
Notice that the position \( y_{*} \) of the barrier separating the half-zippered and the fully zippered states is such that \( y_{*}<a \). Therefore the timescale \( \tau_{1} \) exists only if the initial membrane separation \( y_{0}>y_{*} \).
This assumption seems to be supported by experiments, see Ref.~\cite{Gruget:2018dr}.
It has previously been reported that, upon approach of two membranes devoid of SNAREs, Synaptotagmin exerts repulsive force from \SI{10}{\nano\meter} down to \SI{4}{\nano\meter} where it becomes a repulsive wall \cite{Gruget:2018dr}. According to this result, \( y_{0} \) should range between \num{4} and \SI{10}{\nano\meter}. However it is probably slightly larger  under physiological conditions because of the presence of the SNAREs.
With the parameters adopted in our simulations (see Table~\ref{tab:parameters}), the position of the barrier is \( y_{f}\simeq \SI{4.5}{\nano\meter} \) in accordance with Refs.~\cite{Gao:2012gu,Zhang:2017fm}, see Fig.~\ref{fig:results}. Therefore, in all likelihood, \( y_{0} \) is larger than \( y_{*} \) and our predictions should be valid.

Finally we mention the fact that the timescale \( \tau_{2} \) exponentially decrease with \( N \) seems to be supported by experimental studies reporting sub-millisecond fusion time with \( N=\numrange{3}{6} \) \cite{Hua:2001bb,Sinha:2011gxa,Mohrmann:2010eg}
However, in a recent theoretical study the decay was also found to be exponential but with a much slower decay: the cooperation of at least sixteen SNAREs was predicted to be necessary to reach the physiological fusion time \( \SI{\sim 100}{\micro\second} \) \cite{Mostafavi:2017gd,McDargh:2018jd}. The difference is explained by the fact that the residual work in this study is \(w=\SI{0.48}{\kbt}\) instead of \SI{4.5}{\kbt} in our model, see\eqref{eq:tau_two}.
This difference originates from the assumption made by the authors that the zippering energy of a SNARE complex is entirely dissipated before the membranes encounter the fusion barrier. In other words, the authors have implicitly assumed that after the calcium entry, the zippering of the SNAREpins does not generate any pulling force to assist fusion and concluded that the remaining  residual force is of entropic nature.
Recent direct microscopic observations implying that synaptic fusion involves only six SNAREpins \cite{Li:2019} would appear to invalidate this assumption.

In conclusion,  our model describes membrane fusion by a team of mechanically interacting SNAREpins as  a two stage process. We show that conventional  biochemical and  biophysical measurements cannot be used directly  to predict the
 associated rates and that mechanical modeling is crucial for linking these rates with independently measured parameters.
  Our work emphasizes the importance of identifying  mechanical pathways and  specifying mechanistic feedbacks.  The main conceptual  outcome of our study is the realization   that in the case of synaptic fusion,  SNARE proteins can  perform optimally only if they act collectively.
  The remarkable fact is that when the team is of the optimal size,  such synchronization is not deterred by thermal fluctuations, which guarantees that the collective strike is  simultaneously fast, strong and robust.

Finally, we mention that the synaptic fusion  is  only one  of many biophysical processes involving mechanically-induced collective conformational changes. Other examples include  ion gating   in   hair cells \cite{Martin_2000,Bormuth:2014hh}, collective decohesion of adhesive clusters \cite{Yao:2006de,Erdmann_2007}, folding-unfolding of macromolecular hairpins   \cite{Woodside:2008uz,Bosaeus:2012kp,Liphardt:2001fp} and  folding of ParB-ParS complexes in DNA condensation \cite{Chen:2015cv,Funnell:2016bo}.
In each of these situations one can identify a dominating long-range mechanical interaction making  the  theoretical  framework  developed  in this paper potentially useful.

\matmethods{
\subsection{Rigid membrane assumption}
\label{sub:rigid_membrane_assumption}
We assume for simplicity that the vesicle and the target membranes are rigid, which implies that all the SNAREpins share the same inter-membrane distance \( y \).
This approximation is valid if the characteristic length \( \ell \) associated with the deformation generated by a single SNAREpin is large compare to the size of the SNARE bundle.
We can use the following estimate \( \ell=\sqrt{\kappa/\sigma} \), where \( \kappa \) is the membrane rigidity and \( \sigma \) is the membrane tension.
We have typically \( \kappa\sim\SIrange{20}{50}{\kbt} \) and \( \sigma\sim\SIrange{e-04}{e-06}{\newton\per\meter} \) so \( \ell\sim\SIrange{30}{120}{\nano\meter} \). Since the size of the SNARE bundle is less than \SI{10}{\nano\meter}, our assumption should be valid.

\subsection{Model of a single SNAREpin}
\label{sub:model_of_a_single_snarepin}
We set, for simplicity, that the energies \( e_n \) and \( e_c\) of the SNAREpins in the states $n$ and $c$, respectively, depend on $y$ quadratically, so that
\begin{equation}
	\label{eq:en_ec}
\begin{split}
	e_n(y) &= (\kappa_{n}/2)(y-a)^{2} + e_0,\\
	e_c(y) &= (\kappa_{c}/2)y^{2},
\end{split}
\end{equation}
where \( \kappa_{c,n} \) represent lumped stiffnesses parameters.
We denote \( y* \) the distance where \( e_n(y_{*})=e_c(y_{*}) \), see Fig.~\ref{fig:model}(C).
In the absence of external load (zero force), the stable states are located at \( y=a \) and \( y=0 \). In this situation the entire energy associated with the zippering process is consumed when the SNAREpin reaches state \( c \) at \( y=0 \), which can then be considered as a ground state with zero energy.

The rates \( k_{\pm} \) of the \( n \rightleftharpoons c \) transitions obey the detailed balance relation $k_{+}/k_{-} = \exp\left[(e_c -e_n)/(k_{\mbox{\tiny B}} T)\right]$, with the bias towards the direct  transition  $n\to c$, i.e. \( k_+ > k_-\),     at \( y < y_{*} \) and conversely,  in the direction of the reverse transition $c\to n$ at \( y > y_{*} \), see Fig.~\ref{fig:model}(C)

For simplicity, and following \cite{Huxley_1971}, we consider that the transition from the high to the low energy state occurs at a constant rate $k$, which fixes the characteristic timescale of the conformational change. This assumption could be  easily
replaced by a more adequate one
 at the expense of introducing two additional parameters but with only a minimal impact on the results, see Ref.~\cite{Caruel:2017eba}.
With this assumption and using the detailed balance, we write the transition rates as
\begin{equation*}
\begin{split}
k_-(y)=k,&\; k_+(y)=k\exp\{[e_n(y)-e_c(y)]/(k_{\mbox{\tiny B}} T)\},\;{\rm if}~ y>y_{*}\\
k_+(y)=k,&\; k_-(y)=k\exp\{[e_c(y)-e_n(y)]/(k_{\mbox{\tiny B}} T)\},\;{\rm if}~ y<y_{*}.\\
\end{split}
\end{equation*}

\subsection{Numerical implementation of the model}
\label{sub:numerical_implementation_of_the_model}
The discrete stochastic process associated with the variable \( N_{c} \) was simulated as a two-state Markov chain with a fixed timestep \( \Delta t=10^{-6}\,t_{\eta} \). At each timestep the transition probabilities \( W_{+1,-1,0}\Delta t \) are computed and the next event is chosen based on an acceptation-rejection condition using a random number uniformly distributed between 0 and 1.
The Langevin equation was simulated using a first order explicit Euler scheme.
More details about the computer algorithms can be found in SI Appendix.

\subsection{Adiabatic elimination of the variable \( \boldsymbol{N_{c}} \)}
\label{sub:adiabatic_elimination_of_the_variable_boldsymbol_n__c}
We consider the situation where \( t_{\eta}\gg k^{-1} \): the characteristic time of the conformational changes is negligible compared to the timescale associated with the relaxation of the vesicle's position. In this limit, the conformational state of each SNAREpin can be considered at equilibrium. Therefore for a given position of the vesicle \( y \), the probability of a configuration with \( N_{c} \) SNAREpins in state \( c \) follows the Boltzmann distribution
\begin{equation}
	\label{eq:distrib_nc}
\rho(N_{c};y)=\frac{1}{Z(y)}\binom{N}{N_{c}}\exp
\left\{
-\left[
N_{c}e_c(y)+(N-N_{c})e_n(y)
\right]/\left(k_{\mbox{\tiny B}} T\right)
\right\},
\end{equation}
where \( \binom{N}{N_{c}}=\frac{N!}{N_{c}!(N-N_{c})!} \).
We then integrate \eqref{eq:dynamics_y} with respect to the distribution (\ref{eq:distrib_nc}) and obtain \eqref{eq:effective_langevin}. Since the energy \( E_{\mbox{\tiny snare}} \) is linear in \( N_{c} \), our approximation results in replacing \( N_{c} \) by its average \( \langle n_{c}\rangle(y)=\sum_{N_c}N_{c}\rho(N_{c};y) \) in \eqref{eq:dynamics_y}.
In \eqref{eq:kramers}, the prefactors are given by \cite{Schuss_2010} \( \alpha_{1,2} = \frac{2 \pi k_{\mbox{\tiny B}} T} {a^2 \sqrt{\Phi_{1,2}''(y_{\mbox{\tiny max}}) |\Phi_{1,2}''(y_{\mbox{\tiny min}})|}} \), where \( y_{\mbox{\tiny max}} \) and \( y_{\mbox{\tiny min}} \) denote the positions of the considered barrier and minimum, respectively.

\subsection{Estimation of the mechanical work \( \boldsymbol{w} \)}
\label{sub:estimation_of_the_mechanical_work_w}
In the intermediate state, inter membrane distance \( y_{2} \) is sufficiently lower than the threshold \( y_{*} \), so that the free energy can be well approximated by the energy of the state \( c \). We then write \( \Phi(y)\simeq E_{\mbox{\tiny fusion}}(y)+N\frac{\kappa_{c}}{2}y^{2} \) which leads to the following expression for the energy barrier separating the intermediate state and the fused state,
\[
\Delta\Phi_{2}=e_f\left\{(1-\exp\left[-(y_{2}-y_{f})^{2}/(2\sigma_f^{2})\right]\right\} + N\frac{\kappa_{c}}{2}(y_f^{2}-y_{2}^{2})
.\]
By noting that \( y_{2} \) verifies \( \left.\frac{{\mbox{d}}\Phi(y)}{{\mbox{d}} y}\right|_{y=y_{2}}=0 \), we obtain
\(
\Delta\Phi_{2}=e_f - N w,
\)
with \[ w\simeq\kappa_{c}\big(y_{2}^{2}-y_f^{2}+\frac{\sigma_f^{2} y_{2}}{y_{2}-y_f}\big)\geq0. \]
 Notice that since the energy \( E_{\mbox{\tiny fusion}} \) decays rapidly for \( y>y_f \), the parameter \( y_{2} \) depends weakly on \( N \).
}
\showmatmethods{}

\acknow{We would like to thank Yongli Zhang  for providing us with the data used for the calibration of our model and Ben O’Shaughnessy for stimulating discussions.
This work was supported by a European Research Council (ERC) funded grant under the European UnionQs Horizon 2020 research and innovation programme (grant agreement no. 669612) to J.E.R
}

\showacknow{}