Josephson Junction-Cavity quantum simulation codes and data.

Ben Lang & Andrew Armour
University of Nottingham
ben.lang@nottingham.ac.uk or andrew.armour@nottingham.ac.uk

This simulation data describes the quadratic quantum solutions for the photons in a Josephson-Junction+Cavity system, as described in the paper: Multipartite entanglement in a Josephson Junction Laser. The relevant simulation and graphing codes are also given.

Data is of three catgories. First, classical fixed point amplitudes. Second, the quadratic solutions themselves. Finally, for V=2, N=3 the full quantum solution from QuTiP for comparison. (Note that the code uses V where the paper uses p, but they are the same thing).


Codes : These codes were not written with the intention that anyone else might use them. Consequently they are not very clear. Sincere apologies. They are not going to "just work" as-is, because they need to be placed in directories relative to one in a particular way to load functions from one another. However, if you have a specific aspect of our methodology that you want to check, understand or reproduce then these codes should be very useful material for that. For copy/pasting the parts you want to incorporate into your own project the code should be fine.

Data Categories:

Fixed Point Data:
	Classical fixed point data. These are the classical (mean filed) amplitudes. Only collected for V=1, and N=1 to 9. A re-scaling allows these fixed points to represent any situation. For example with V=3, and N=9 one takes the V=1, N=3 data for the mode amplitudes, and (after scaling) the three mode amplitudes give the amplitudes of modes 3, 6 and 9. The nessesary scalings are described in Phys. Rev. B 107, 144509, and an implementation can be seen in the code.
	
	The folder "Fixed_Point_data" holds the fixed points.
	N, V as in file names. "Attempt=1" is historical and unimportant. (For the curious, when going above threshold multiple fixed points exist, so multiple initial conditions "attempts" should be tried to find them all. However, the present data concerns the regime where there is only one fixed point where additional attempts are unnessesary.)
		Format: 
		All files in ".npy" format. Can be loaded by numpy.
			fp_sols - Each row is a list of amplitudes for the N resonant harmonics, up to scaling factors (as described above, see code).
			fp_stabs - Each row is True if the fixed point is stil stable, considering only instabilities in the resonant harmonics.
			fp_LambTs - Each row is the Lambda value at which the equivalent rows of fp_sols and sp_stabs were calculated. These re-scale to EJ based on V.

	The code "Run_quadratic_solve_N_geq_2p.py" shows how the fixed point data can be loaded and re-scaled for use.
	The code "Semiclassical_fixed_points.py" shows how new fixed points can be calculated.
	
	Also inlucded here in "V_not_1" are classical fixed point solutions where V was not 1. Calculations of this type are needed to find discontinuous symmetry breaking, where a new stable solution forms "out of the blue" far in parameter space from the existing solutions. (See Phys. Rev. B 107, 144509). For the two data sets given the location of the discontinuous transition was relevant as we wanted to restrict our analysis of the entanglement to the region before the discontinuity.
	Additionally, the folder "gamma_sqr" contains fixed point data, but where instead of the loss rate being assumed equal for all modes, it increases quadratically (so that the second mode has 4 times the loss of the first and the third 9 times). Reality is probably somewhere between the two extremes (all loss same, and quadratically increasing loss). See appendix C1 of Phys. Rev. B 107, 144509 for further details.
	
Quadratic Solution Data:
	Information about the solutions to the quadratic (Gaussian) model around the fixed points.
	For historical reasons this data is split between two folders, "Delta_over_root_n" and "Multi_alpha". The rule is that if there is exactly one resonant harmnic (IE if V < 2N) then it goes in "Delta_over_root_n", otherwise in "Multi_alpha".
	(Specifically, it was assumed initially that having 2N > V would not be possible, and the case where 2N is > V was added as an experimental feature that was eventually found to work well. Hence many code file names contain "N_geq_2p" meaning they are the versions of the code to work in that regime.)
	
	Inside these folders the sub-folders are named according to N and V. The files found inside (all .npy) are:
		Lambs.npy - The drive strengths at which all the other data were collected.
		Stabs.npy - Whether or not the fixed point has lost stabiity at each Lamb.
		covmats.npy - The covariance matrices at each Lamb. If $R = (\hat{x}_1, \hat{p}_1, \hat{x}_2, \hat{p}_2, \hat{x}_3, \hat{p}_3,...)$ the covariance matrix elements are given by: $C_{nm} = \langle R_n R_m + R_m R_n \rangle - 2 \langle R_n \rangle \langle R_m \rangle$. These will be nonsensical if the stability has been lost. Each covmat is a matrix of size 2N by 2N. For a coherent state it is the identity.
		free_squeezes.npy - Half the smallest eigenvalue for each covmat.
		expects_quad.npy - <a_n a^{dag}_n> for each of the N modes for each Lamb.
		gen_group_entangles.npy - For each lamb has entires eequal to the number of K-groups. Usually "0" or "False". Returns a Tuple if there is genuine multipartite entanglement accross ALL the modes in that K-group. Tuple is (Ture, Method, UV values used to find entanglement).
		sub_gen_entangles.npy - Similart to above, however, considers not just each K-group, but every possible fragment of a K-group containing 3 or more modes. Ordering is determined by the "get_pods" function in "tripartite_range_figure.py".
		extended_*.npy - The solver initially guesses how high it will need to go in Lambda to reach the instability, and divides this into 100 steps. So the above files are 100 steps long. However, if it does not reach the instability in its sweep (it under-guessed) an extension is carried out. The extended files have the same form as the equivalent normal files.
		(Inconsistently) - Plots, summarising the key data for that V, N combination.
	
	The code "Run_quadratic_solve_N_geq_2p.py" is informative of how this data is generated. "extend_quadratic_solve_N_geq_2p.py" is similar, showing how a sweep can be extended.
	The relevant functions used are in "quantum_quad_solver.py"
	For an example of how the data can be read and processed see "tripartite_range_figure.py"

	Additionally, "Multi_alpha_quad_loss" contains the same types of data, but for the situation where the losses are not the same for every mode, but increase quadratically. We did not carefully study how the entanglement behaved with these "soft cut-offs", and this data can reasonably be described as as pre-empting referee comments.
	
Non-Gaussian quantum data:
	For N=3, V=2 a single set of calculations were run using the "QuTip" python package (https://qutip.org/). These do not make a quadratic approximation, so give a way of testing it. However, they were very time consuming.
	
	Quitp_solving contains all the codes/data relevant.
		Delta_over_root(n)_QUTIP_Compare contains the data.
		V=2, N=3, [15,6,15] contains the QUTIP data using a Hilbert space size of 15,6,15 (relatively small and fast).
		V=2, N=3 contains QUTip data for a larger state space ([28, 6, 26]), and only for the last 4 values of Lambda.
		
		The QuTip data is:
			Q_Lambs - The Lambda values at which QuTiP data was collected.
			Displaced_SS - The steady state density matrices (displaced from the fixed point to the origin) at those Q_lambs.
			Displacements - The classical fixed point amplitudes at those Q_lambs, that were used to displace the quantum equation of motion.
			aops_disp - The displaced anhilation operators of each of the 3 modes. (Inlcuding tensoring with identiy acting on the other modes)
			Displaced_expects - <a_n a^{dag}_n> for each of the N modes for each Lamb (displaced picture).
	
	"make_qutip_quad_compare_fig.py" is a good example of how this data can be loaded and used.
	"Displaced_RWA_Ham.py" shows how this data was generated. Note that it uses many functions from "RWA_Ham.py". Also note that, running with the [28,6,26] sized state space will be extremely time consuming on a normal pc (days).
	

K_plot shows which subspace looses stability first for each N, p combination. This can be helpful for making sense of the other data.

Figure_making contains code that makes useful plots with the data. If you want to make use of the data, canibalising sections of these codes is probably the best way to start.


For questions or issues please contact either ben.lang@nottingham.ac.uk or andrew.armour@nottingham.ac.uk. 


Have a superb day,

Ben