Breaking of Josephson junction oscillations and onset of quantum turbulence in Bose--Einstein condensates

We analyse the formation and the dynamics of quantum turbulence in a two-dimensional Bose-Einstein condensate with a Josephson junction barrier modelled using the Gross-Pitaevskii equation. We show that a sufficiently high initial superfluid density imbalance leads to randomisation of the dynamics and generation of turbulence, namely, the formation of a quasi-1D dispersive shock consisting of a train of grey solitons that eventually breakup into chains of distinct quantised vortices of alternating vorticity followed by random turbulent flow. The Josephson junction barrier allows us to create two turbulent regimes: acoustic turbulence on one side and vortex turbulence on the other. Throughout the dynamics, a key mechanism for mixing these two regimes is the transmission of vortex dipoles through the barrier: we analyse this scattering process in terms of the barrier parameters, sound emission and vortex annihilation. Finally, we discuss how the vortex turbulence evolves for long times, presenting the optimal configurations for the density imbalance and barrier height in order to create the desired turbulent regimes which last as long as possible.


I. INTRODUCTION
The Josephson junction (JJ) is an experimental set-up designed to showcase the Josephson effect [1]. This quantum mechanical effect, which describes particle tunnelling through a barrier and periodic oscillations, is well studied in the context of Bose-Einstein condensates (BECs) in both theory [2][3][4] and experiments [5][6][7][8][9]. In this article, we discuss the dynamics when a BEC JJ is pushed to its limit and "goes bad". Namely, we study the regime where the periodic oscillations in the superfluid density break down and quantum turbulence arises in the system. Current methods to create vortex turbulence include optical spoons [10] and shaking confining traps [11]. Both of these mechanisms have drawbacks, the methods require much energy to create a single vortex in fluids with high density. Moreover, if the density is low the resulting vortices annihilate very quickly; therefore, the creation rate will be similar to the annihilation rate. We propose a method to create vortex dipoles with initial imbalance and sustain them by using a Kibble-Zurek [12,13] like mechanism to prolong the vortex turbulence, that is, create vortices in a region of low density and then increase the density in a controlled manner to maintain the topological defects and decrease the relative strength of the acoustic waves.
The Josephson junction consitsts of a barrier or weak link separating two wells of superfluid or superconductor. We consider two wells separated by a potential barrier with an initial density imbalance between the wells. As the system evolves the fluid moves through or over the barrier, which leads to oscillatory dynamics. We show that when the initial density imbalance is pushed to high values, the regime of regular oscillations break down. The system then enters into chaotic behaviour with interesting nonlinear dynamics consisting of chaotic motion of vortices coupled with turbulent acoustic waves following the break down of a soliton train caused by a dispersive shock [14]. Studies on critical parameters for vortex generation have been undertaken [15]. However, these do not consider high numbers of vortices an turbulence.
In place of the predictions by Josephson, we see interesting turbulence characterised by a separation of acoustic and vortex turbulence. Such chaotic dynamics appear when an initial train of solitons are formed then break down, causing the generation of vortices. This is due to the instability of quasi-1D solitons.
We demonstrate that the Josephson junction is a nonlinear system and that by tuning experimental parameters, we can produce rich and controllable non-linear behaviour, which gives an ideal set up for turbulence. Readily available experimental apparata in BECs [5,11], atomic vapors [16] and photorefractive crystals [17] can implement such a system. We show that a crucial element in the turbulent dynamics in such a system is due to the interaction of vortex dipoles with the barrier. These are responsible for the separation of turbulence within the two wells, with one retaining most the vortices and the other containing weak acoustic wave turbulence [18]. Certain parameters (for instance the incidence angle) control the ability for vortices to cross the barrier [19,20].

II. THE MATHEMATICAL MODEL
To model theoretically the JJ we perform direct numerical simulations (DNS) of the Gross-Pitaevskii (GP) equation. The GP equation describes the dynamics of a BEC made of dilute ultra-cold gas of bosons [21,22]. For simplicity we consider the case of a quasi-two dimensional BEC, that is we simulate the GP equation in two spatial dimensions as well as time. Lengths are expressed in units of the healing length ξ = / √ 2mρ 0 g 2D where ρ 0 is the mean density, g 2D is the effective two-dimensional interaction constant between bosons of mass m, and µ is the chemical potential of the system. Time is rescaled by ξ/c, where c = ρ 0 g 2D /2m is the speed of large scale density/phase fluctuations (sound) in the bulk. The external potential V is given in units of ρ 0 g 2D . The non-dimensional form of the GP equation reads as follows: The complex wave function ψ(x, y, t) is the BEC order parameter. The BEC order parameter can be expressed in fluid like variables via the Madelung transformation ψ = ρ(x, t)e iφ(x,t) , where ρ(x, t) = |ψ(x, t)| 2 and v(x, t) = 2∇φ(x, t) are the density and velocity of the superfluid respectively. In our dimensionless variables we rescale density by the initial density ρ 0 .
Our aim is to model an elongated JJ domain, large with respect to the healing length, in order to observe the formation of several quantised vortices and, eventually, fully developed quantum turbulence. We thus choose a twodimensional spatial domain x = (x, y), with x = [−256ξ, 256ξ] and y = [−128ξ, 128ξ], setting a computational uniform grid of spacing 0.25ξ. We have Dirichlet boundary conditions with ψ = 0 at the boundary, this is effectively confining the fluid in an abrupt rectangular tap with an external potential at the boundary with infinite strength. The JJ barrier is modelled using an external potential V JJ (x, t) that separates the domain into two equally sized boxes, labelled B L and B R corresponding to the position left or right of the potential. The potential is given by a Gaussian function centred at x = 0 and stretched along the entire y-axis. To create the initial superfluid density imbalance needed to trigger the JJ oscillations at t , an extra non-zero external potential (almost) uniform V d is present in the right box. Mathematically, the JJ barrier thus results in where V 0 and σ control the intensity and the width of the JJ barrier, respectively, V d sets the initial density imbalance and H(·) is the Heaviside step function. In all the simulations reported in this work we keep σ = 1.2 so that the width of the barrier is always similar to the healing length in the system. To create the initial condition, we use imaginary time propagation for t < t . This method involves evolving the GP model with the substitution t = −iτ while imposing conservation of the superfluid density to effectively minimise the energy of the system. Once the desired energy stagnation has been reached, we call the obtained state the ground state, set t = t and evolve the system in real-time. Because of the sudden jump in the V JJ at t = t , the initial condition is no longer the ground state; therefore, dynamics follow. The GP equation (1) is integrated using the 4th order central finite difference method in space and a Runge-Kutta 4th order method time-stepping scheme. Due to the nature of the method, only waves with wavelength of the order of half the healing length or higher are well resolved; smaller waves will be numerically dissipated. This natural dissipation of high energy high-frequency waves is ideal for such a study as dissipation occurs in experiments similarly due to the interaction with the thermal cloud and/or bosons losses due to finite-amplitude confining potential [23,24].

A. Measurable quantaties
We define the relative superfluid density imbalance using the total superfluid mass per box B i where i is an index for the box left or right of the separating potential, Note that the total density N = N L + N R is an integral of motion and it is numerically conserved in all simulations up to numerical dissipation. Analogously, the energy per box reads The total energy E = E L + E R is also an integral of motion, but due to the intrinsic high-frequency numerical dissipation, its value decreases in time as much as 28% when Z 0 = 0.49 around 10% when Z 0 = 0.88. For the results in section III C, where there is no acoustic turbulence, the energy is conserved to 0.0015%. The energy is naturally decomposed in (5), with the second term corresponding to energy from the external potential and the third term the internal energy of the fluid. We can further decompose the first term into kinetic and quantum energy by applying the Madelung transformation ψ = √ ρe iφ to the first term where the first term corresponds to the kinetic energy density and the second term is the so called quantum energy density. By performing a Helmholtz decomposition on the kinetic energy density we can further decompose into the compressible and incompressible energies. That, is ε kin = ε c kin + ε i kin , where the incompressible component of the energy corresponds to the vector field satisfying ∇ · ( √ ρv) i = 0. Further details on the calculation can be found in [25]. Thus the energy can now be written The total incompressible energy E in kin = B L +B R ε in kin dxdy is a measure of the energy in large scale incompressible potential flow and vortices, whereas the total compressible energy E c kin = B L +B R ε c kin dxdy is the energy in the acoustic component.
It is also convenient to define the local healing length for each box as and we will call the natural healing length of the system ξ, that is, the healing length if the initial density imbalance is set to zero. The length of each box in units of ξ is given by L.
Finally, it is instructive to measure the total number of vortices in the system versus time, with N V L and N V R the number for the left and right boxes respectively. Each quantised vortex is numerically identified using the pseudo-vorticity defined as follows, with where j is the density flux. Then we find the maxima in simply connected regions ignoring the field below a chosen cut-off value, see [26] for further details. Ghost vortices are phase fluctuations in large regions where ψ is close to zero, these do not show the same dynamics as hydrodynamical vortices. We add an extra filter to our vortex tracker, namely we only consider vortices with substantial density surrounding them, this is to remove the ghost vortices and to track only the hydrodynamic vortices. The numerical scheme calculates the average density around any point identified by the vortex tracking routine and discards the vortex if the average is below a threshold value.

B. Creation of vortices
In the GP model, 1D dark or grey solitons are unstable to transverse perturbations in two spatial dimensions, the instability is known as the snake instability. This instability is a result of the speed of a soliton being proportional to its amplitude. The instability can be understood as a result of smaller solitons having larger speeds. A small transverse perturbation introduces a local difference in speed of the soliton. Such a difference in speed will cause the soliton to bulge in the direction of motion if the perturbation is negative, or the opposite direction if the perturbation is positive. For instance, if we introduce a bulge in the direction of motion, since the soliton will move perpendicular to its tangent there will be a focussing effect on either side of the perturbation. Since the speed of the soliton is reduced when its amplitude increases, the focused parts of the solitons slow down producing an inverted bulge. The process then continues along the length of the soliton with bulges and inverted bulges forming along the soliton, for more details and mathematical analysis see [27]. When a grey-soliton's amplitude becomes as large as the density around it, i.e. points at which ψ = 0 will appear: phase defects in the form of vortices are nucleated. In this section we will discuss how we take advantage of this instability to produce vortices.
The potential is such that a train of solitons is produced within a dispersive shock wave in the right well which has low density N L (0) < N R (0), this can be seen in Fig. 1 which shows an example simulation for the entire domain. Fig. 1(a) shows the production of the train of solitons, seen as the stripes in the low density region. The solitons will then decay by the snake instability into alternating signed vortices, the process begins at the boundary and can be seen in Fig. 1 (b) with later stages in Fig. 1 (c). Due to the large number of vortices and to the fact that the local density is small, the initial vortices have large cores and do not interact like hydrodynamic vortices. These vortices are often referred to as ghost vortices and are not counted by the tracking algorithm. The continuous flow of solitons carrying density into the right well reduces the local healing length, which in turn transitions the ghost vortices towards hydrodynamic vortices. After the ghost vortices are produced in Fig. 1 (b), we see chaotic motion with a proliferation of hydrodynamic vortices in Fig. 1 (c) when the density becomes larger due to the fluid flux from the left box. As the process continues, the healing length tends to ξ in both boxes, the healing length when N L = N R . Some opposite signed vortices annihilate which results in continuous decay of the total number of vortices. As a result, the remaining vortices become even better formed as hydrodynamic because the mean distance between them becomes much greater than ξ. Later stages of the dynamics are shown in Figs In order to show the snake instability better, we present an example of the early stage 2D density and phase fields with potential strength V 0 = 1.5/µ and initial imbalance Z 0 = 0.88. We zoom in on the subregion of B R in Fig. 2 (a), the region depicted by the red rectangle in Fig. 1 (a). The train of quasi 1D solitons within a dispersive shock region is seen in Fig. 2 (a). Fig. 2 (b) shows the snake instability forming on the solitons. As they travel, the solitons begin to snake until they break up into a chain of ghost vortices with circulations of alternating signs, which then interact to form vortex turbulence. The ghost vortices can be seen in the corresponding phase plots. For instance, in Fig. 2 (f) we see discontinuities in the phase where the phase winds from −π to π around them. The ghost vortices correspond to the dark (blue) circles in Fig. 2 (b). They are not counted as in the number of vortices or marked by the black and white circles corresponding to the well-formed hydrodynamic vortices.
In the case of many atoms, where the Laplacian term in the Gross-Pitaevskii equation can be neglected, we can analytically find the stationary profile for the wave function. The latter is known as the Thomas-Fermi (TF) profile. We choose to take a slice along y = 0 as there is no y-dependence on the potentials. When the fluid has enough energy to flow over the barrier, that is when the TF profile |ψ(x) T F | 2 > 0 for all x, the dynamics are different to that of the classical Josephson junction as the flow is not only a consequence of tunnelling but also due to the fluid which has enough energy to pass the barrier. The stationary TF profile is given by:  Fig. 1(a). Panels (a,e): t = 30ξ/c, (b,f): Thus the barrier TF width, W T F , is given by: From this we can see that flow has no tunnelling when V 0 < ρ 0 + V d . Once the vortices are introduced, confining them to the right box only should make them interact more due to the smaller inter-vortex distance. The main process in which vortices spread out is by forming dipoles which move away from the vortex bulk at a nearly constant speed. These vortices can penetrate the barrier in certain cases, namely when V 0 is small and/or the vortex dipole is fast. The barrier width W T F will also affect the ability for vortices to penetrate the barrier. Such dependences of the penetrability on the barrier and vortex properties are not obvious, so we will now present a study to classify and quantify different outcomes of the dipole-barrier interactions.

C. Vortex dipole scattering off the barrier
The interaction between quantised vortices and the JJ barrier plays a vital role in the dynamics discussed in section III B and is an interesting problem in itself. The barrier can trap vortices as well as assist in their annihilation. In certain regimes of 2D vortex turbulence, vortices tend to couple into vortex dipoles [18]. This process is the result of random vortex motion; it includes inter-vortex collisions which can re-couple or scatter vortex pairs. As a result of this motion some vortex pairs will move away from the turbulent bulk. As we have boundaries, the vortices will be incident either on the outer boundaries or the barrier. If the dipole is incident on the outer boundary it, will get split into the two vortices moving along the boundary in opposite directions. If the dipole is incident on the barrier, it can be transmitted, annihilated or trapped, depending on the barrier height and the dipole size.
In this subsection we present a study of a vortex dipole interacting with a JJ barrier where no initial density imbalance between the the left and right boxes is present, that is, Z 0 = 0 and V d = 0. Initially, we position a vortex dipole centred at (−25ξ, 0), and define θ as the angle between the x-axis and the direction of propagation of the vortex dipole. Clearly by symmetry, we expect the dynamics to be mirror-symmetric with respect to the x−axis i.e. with respect to the change θ → −θ. The vortices are initially separated by a distance of d 0 , and the vortex with positive circulation is in the upper-half plane so that the dipole moves towards the positive x-direction. Examples of the resulting vortex trajectories are shown in the different panels of Fig. 3. In Fig. 3 we see that in all examples the vortices lose energy to sound during the interaction with the barrier. We quantify the amount of sound emitted by calculating the change in the incompressible energy We measure this change in energy at a final time t f after the interaction has happened. As a criteria for determining t f we choose the following three criteria: (i) the vortex dipole passed the line x = 25ξ corresponds to Fig. 3 (a,d); (ii) either vortex becomes within 5ξ of the system boundary corresponds to Fig. 3 (b); (iii) if case (i) and (ii) are not fulfilled we allow a maximum time of t f = 750ξ/c. We assume that the vortices have annihilated if they do not fulfil case (i) or (ii) after such a long time, so case (iii) corresponds to vortex annihilations an example is shown in Fig. 3 (c). Fig. 3 shows four different examples of the dipole-barrier scattering for different values of the scattering parameters V 0 , d 0 and θ. The images show the superfluid density plots after the scattering with the dipole with the trajectories overlaid. For the dipole-barrier scattering experiments presented in this subsection we have chosen to shorten the L x side and increase the L y side of the JJ system, compared with the one presented for instance in Fig. 1, in order to give the vortices more space to interact with the barrier. Fig. 3 (a) shows a the path two vortices take when passing a barrier. Notice that the vortices after the interaction are closer together. As the vortices' motion was perpendicular to the barrier there was no deflection. Fig. 3 (b) shows the case of a dipole not being able to pass the barrier. This corresponds to case (ii). The vortices separate from one another and move along the barrier. The vortices effectively see images of themselves in the barrier. This motion is similar to that when a vortex-pair is incident on an external boundary, where the boundary conditions are similar to that of an infinite barrier. Fig. 3 (c) shows an annihilation, this is case (iii). Fig. 3 (d) shows the vortices being deflected during the interaction with the barrier. In this example the vortices were not moving perpendicular to the barrier but had an incidence angle θ = 0.4π. We perform two sets of simulations for two different incidence angles θ = 0, 0.4π, sweeping over different vortex separations d 0 , and vortex barrier heights V 0 . We observe and classify for each set of parameters, what kind of dipole barrier interaction takes place, (i), (ii) or (iii), and also measure sound released in the interaction. The results are presented in Fig. 4. Firstly we calculate whether or not a vortex dipole of separation d 0 will pass the JJ barrier of given strength V 0 ; the cases when the dipoles pass are marked with asterisks.
If the dipole cannot pass, this can be due to two reasons: the dipole is annihilated by the barrier producing sound (case (iii) marked by squares) or the interaction of the dipole and the barrier is not over when the boundary effects become relevant(case (ii) marked by circles). In the latter case, it is not possible to say whether the dipole would have annihilated in an infinite domain or not, hence it is a consequence of the finiteness of the system. We see that when annihilations happen, there is more sound energy (measured by the change in the incompressible energy) released from dipoles with larger separations. Fig. 4 shows that there are clear connected regions in which each of the possible cases happens.
As we mentioned in the introduction our aim is to explore the best parameters for vortex turbulence. Annihilations are key events which determine the vortex decay rate, which will be further discussed in section III B. We see that the barrier can cause annihilations. In Fig. 4 (a) we see that the annihilations are numerous and they also correspond to a high emission of sound. As well as removing vortices from the system, the extra sound is known to increase the vortex decay rate [28]. On the other hand, the highly energetic sound produced can penetrate the barrier and be spread over the adjacent box which is void of vortices. As a result, the sound is distributed over twice the area (both the wells). Thus, the second well acts as a sound absorber i.e. as an effective heat sink. By introducing a finite angle of incidence θ, the amount of vortex annihilations is greatly reduced. At the same time vortex splitting becomes much more frequent. However, the amount of vortices that pass the barrier does not change much. We see that small dipoles (d 0 < 7) are annihilated by the barrier in the region of the barrier strength V 0 ∼ 1. For the region V 0 = 1.3 − 1.6 with large dipoles (d 0 > 7) we see that the dipoles separate. Compared to the splitting region in Fig. 4(a) in this case they emit much more sound. The vortex-barrier interaction discussed here is more simple than the interaction when the densities are unequal, the background condensate is saturated with sound, and there are more than two vortices involved in the interaction. The results are informative as a rough measure of the vortices ability to cross the barrier.

Optimal parameters
We now discuss the optimal choice of parameters to produce vortex turbulence. Our aim is not only to produce the highest number of vortices, but to also minimise the secondary by-products of the method, namely, sound and large density waves. We also introduce the mean number of vortices over timē and for the left and right box,N V L andN V R respectively. We use the mean (opposed to the maximum) as a measure as it also takes into account the sustainability of the vortex turbulence. We not only want to create many vortices we want them to persist for as long as possible. Another measure we use to classify the quality of the turbulence is the amount of interfering compressible waves, these account for the large scale density sloshing and the small scale acoustic component. Obviously, it is desirable to minimise such compressible motions.
The mean number of vortices produced by the proposed method depends both on the barrier strength and the initial imbalance. In Fig. 5 we show the effects of varying the two control parameters, the barrier strength V 0 and the initial imbalance Z 0 . In Fig. 5 (a) we see that for a large imbalance (Z 0 = 0.88) increasing V 0 monotonically reduces the mean number of vortices, and for a large enough V 0 no vortices will be produced; this is due to the barrier disrupting the creation of solitons by reducing the rate at which fluid can cross the barrier. On the other hand, for a lower initial imbalance (Z 0 = 0.49), Fig. 5 (a) shows that there is a clear maximum for the mean number of vortices in B R at around V 0 /µ = 0.6. This is due to the mechanism discussed in section I: a steady influx of density 'freezes' the vortices in the right well and the potential reduces oscillations of the density imbalance which would otherwise strongly interact with the vortices. For a high initial imbalance (Z 0 = 0.88), the mean number of vortices is higher without a barrier as shown in Figs. 5 (a). However, as we will discuss below, the vortices are accompanied by more small-scale (acoustic) and large scale (sloshing) compressible waves, which is an undesirable effect if we want 'clean' vortex turbulence. Even when the initial imbalance is high in Fig. 5 (a) and there is no barrier (V 0 = 0), the vortices are evenly distributed over both boxes. As a consequence the amount of vortices in the right box B R , is not so much larger than the cases with higher values of V 0 . When V 0 is fixed, by increasing the initial imbalance we reach a plateau in the mean number of vortices in Fig. 5  (b). This plateau, along with the increase in compressible energy in Fig. 5 (c), indicates that after a certain imbalance the energy is more swiftly converted into compressible sound waves. Using this insight, we propose that it is preferable to choose a lower imbalance such that the vortex dynamics are cleaner, that is, there is less acoustic turbulence and large-scale density sloshing interacting with the vortices.
A secondary effect also complements the longevity of the vortex turbulence. The highly-energetic small scale sound, produced during the vortex creation and subsequent interactions, can easily pass the barrier. The sound energy from the vortex turbulence of the right well is then distributed over twice the area, reducing the interaction with the vortices and, therefore, slowing down their annihilation rate.
In the absence of a barrier, many vortices are produced; however, we identify two critical issues with the turbulence that follows. The first is that the local healing length oscillates for a long time, see Fig. 6 (a). The fluctuation of the local healing length causes vortices to annihilate. Secondly, the large density waves are seen to interact with vortices, this adds additional complexity to the interactions, and it is not clear what the effect on the vortex interactions this will have. Introducing a barrier addresses both of these problems as seen in Fig. 6. When a barrier is present, as in Fig. 6 (b), the oscillations are dampened, with only small oscillations remaining once the boxes have equilibrated. This is due to the amplitude of the large wave simultaneously being reflected and transmitted on each barrier interaction. This multiplies the number of waves decreases the size of the local wave amplitude. Also, it causes the density to fill the right box more smoothly. The smoother descent is due to barrier reflecting more of the wave when the barrier is stronger. That is, the amplitude of the transmitted wave is reduced on each wave-barrier interaction, therefore, the density flux across the barrier is also reduced; this essentially dampens the overshooting of the oscillations. This can be seen in Fig. 6 (c) where we plot the standard deviation of the oscillations, σ, against barrier strength: this reduces quickly when the barrier strength increases and a transition seems to occur at a value close to V 0 /µ = 1. We choose the time period to measure the standard deviation over to be the same for each simulation. The period is chosen to be after the all the simulations have reached Z(t) = 0 for the first time.

Vortex decay rates
Recent discussions [29][30][31], indicate that the number of vortices in a homogeneus condensate decay as t −1/3 , this corresponds to four-vortex interactions. The arguments in [29][30][31] where α here corresponds to the number of colliding vortices causing an annihilation. Equation (15) is a crude approximation which does not take into account any correlation between the vortices, nor does it take into account spatial inhomogeneity in the system. Following this approach, for two-vortex annihilations we have N v is proportional to t −1 , for three-vortex interactions -to t −1/2 and four-vortex interactions -to t −1/3 .  Fig. 3 (a,b,d) (however, the barrier can be chosen such that the vortices cannot penetrate or so that only vortex dipoles of a certain size can exist in the left well, this is discussed further in section III C); (4) vortices annihilate from interacting with the background sound. Fig. 7 shows examples of the evolution of the number of vortices in a log-log plot. In the figure the red line represents the number of vortices in the right box N V R , the blue the number in the left box N V L and the black the total in the entire domain N V R + N V L . We overlay a vertical line for which the imbalance has become almost zero (Z < 0.05) for the first time in each simulation, this is an indicator of when the initial vortex creation period has ended.
We first focus on the decay rate in the later stages of the dynamics where we expect to see the exponents predicted above. We present in Fig. 7 cases with two values of initial imbalance, Z 0 = 0.49 (a-c) and Z 0 = 0.88 (d-f). In the cases with no barrier (a, e) we see that the vortices move freely between the two boxes. In Fig. 7 (a) in particular, we do not produce enough vortices to see vortex turbulence, thus we do not see a clear decay rate. In Fig. 7 (e) homogeneous vortex turbulence develops quickly with an almost equal amount of vortices in each box. In this case we see a decay closer to t −1/3 which is predicted for a four-vortex process. In Fig. 7 (f), for small value of V 0 , see also a decay closer to t −1/3 , whereas a t −1 decay is seen for higher barrier height in Fig. 7 (g), which corresponds to a two-vortex collision in the logistic equation.
Clearly the amount of vortices in the left box decreases with the barrier height. It also appears that the number of vortices in the right box remains almost the same whilst the number in the left tends to zero. We also note that for cases with a high initial amount of vortices there is a steeper decay rate. We propose that some transition in the dominant type of the vortex collision occurs at a vortex density which is dependent on the mean density of the vortices. Indeed, it is natural to think that the higher-order vortex collision (e.g. four-wave) dominates over the lower-order collision (e.g. three-wave) for higher vortex densities and vice versa. The mean vortex density is dependent on the barrier height in a non-trivial way as shown in section III C.
It is difficult to distinguish the best fit from the decay in Fig. 7 by visual inspection. We instead calculate the the best fit numerically and present the results in Fig. 8. We choose a value t s such that all the simulations are safely in the late stage vortex decay regime and calculate the best fit of the exponent over the same range in all the simulations.
In Fig. 8 (a) we present the exponent plotted against increasing initial imbalance Z 0 . We notice that the decay rate is faster for a higher imbalance. We conjecture that this is due to the higher mean vortex density as well as an increased acoustic component (see Fig. 5 (c)) which interacts with the vortices helping them to annihilate. In Fig. 7 (b) as the barrier height increases the total number of vortices (black line) and the vortices in the right box (red line) converge due to all of the vortices being in the right box. We also see in Fig. 8 (a) and (b) that the decay rate fluctuates. However, for most of the parameters, it seems to be steeper than t −1/2 and shallower than t −1 which may be due to the combination of all of the annihilation mechanisms. For instance, the decay rate may be explained by the interaction with the boundary and barrier, which can aid annihilation. For instance in Fig. 9 we see a vortex dipole scatter of a third vortex (a three-vortex precess) until the dipole is small enough such that the boundary or the barrier will annihilate them (a two-vortex process). Also in Fig. 9 (c,d) the red ellipses show the dynamics during (c) and after (d) the interaction with the barrier. We highlight that the rarefaction pulse caused by the annihilation is directed into the left box, thus the pulse is not likely to interact with other vortices causing more annihilations. We note that by increasing dipole size, we do not necessarily increase the chance of annihilation; see Fig. 4 (a). If we take V 0 = 1.6/µ we see that smaller dipoles are trapped whereas larger dipoles annihilate.

IV. CONCLUSIONS
In this paper, we have explore the use of a Josephson junction set-up for generating BEC vortex turbulence. We have shown that the generation and decay of vortices in a Josephson junction BEC configuration, can be altered by controlling the barrier height and the initial density imbalance parameters, which can be readily applied to existing experimental apparatus. We discussed the critical advantage of this method for generating vortex turbulence is the creation of vortices in low-density regions, with the density then being increased to solidify (shrink) the vortex cores. We have provided ranges of parameters to tune the barrier to produce a certain optimal number of vortices and also shown parameters which allow for vortex penetration. We showed that for a higher imbalance we produce more vortices. However, most of these vortices decay quickly leaving much of acoustic noise and large-scale density sloshing, and for this reason it may be preferential to have a smaller initial imbalance. For instance, if one would like to create vortices confined to a single box with the acoustic component less than 10% of the total energy one would choose an initial imbalance of Z 0 = 0.6, with a barrier with strength (measured in chemical potential µ) in the range of µ to 1.5µ with the width of the order of the healing length. Although our simulations are larger than current experiments, there are experiments that are not so much smaller [32,33]. Another disadvantage of the vortex turbulence without the barrier is the large wave-vortex interaction. A key consideration is the control of the amount of sound energy released as part of the snake instability. Such sound can be shared over the two boxes while confining the vortices to one box, and thereby separate the vortices and sound in physical space better so that the sound does not adversely affect the vortex dynamics. We show that the interaction of vortices with a quasi-1D barrier is non-trivial and that the barrier also can work as another effective mechanism for vortex dissipation. This is particularly interesting in terms of using such a barrier to filter vortex turbulence and produce vortices of a certain size, a more detailed study will be the subject of future work. In this work we focused on the vortex turbulence; however, the acoustic turbulence present is also worth further study.