Skip to content

Conversation

@jviatteau
Copy link

@jviatteau jviatteau commented Aug 11, 2024

Added code to solve equations of motion for frequency conversion based on the 2013 paper by Andreas Christ, "Theory of quantum frequency conversion and type-II parametric down-conversion in the high-gain regime".

Created new folder FrequencyConversion containing new code files.

Added FC_solve_EoM.py which provides functions to calculate the numerical solution to the FC equations of motion.

Added test_solve_EoM.py and example_solve_EoM.ipynb as tests and examples for the functions provided by FC_solve_EoM.py.

Added test_fermionic_bloch_messiah.py which provides tests for the fermionic Bloch-Messiah decomposition (to be added later).

@MHoude2
Copy link
Collaborator

MHoude2 commented Aug 12, 2024

Few comments for the "Fc_solve_EoM.py".

import numpy as jnp
from scipy.linalg import svd
from scipy.linalg import expm

why jnp and not just np?
Numpy has its own SVD, why use the scipy one?

If you are done loading things to the PR, you need to add reviewers. I can add Nico if you would like.

We should also probably not include the Fermionic Bloch-Messiah since it fails and wait until we get it to work properly.

@jviatteau
Copy link
Author

I changed the jnp back to np. That was an artifact from the time when I was asked to use jax in my code, so that you could run the code faster on GPUs. I switched back to ordinary numpy (without switching jnp to np) when I needed to work with double precision with the fermionic bloch-messiah function.

I also replaced the scipy svd by the numpy one.

I removed fermionic_bloch_messiah.py but kept the test file, so you have it somewhere when the time comes to test it.

Also, I may not have the permissions or something, but it seems I can't add reviewers myself. I would then be grateful if you could indeed add Nicolas as reviewer.

Comment on lines 154 to 159
Uc = np.diag(np.exp(1j * dw * kc * L)) @ Uc
Vc = np.diag(np.exp(1j * dw * kc * L)) @ Vc
Va = np.diag(np.exp(-1j * dw * ka * L)) @ Va
Ua = np.diag(np.exp(-1j * dw * ka * L)) @ Ua
return np.block([[Ua, Va], [Vc, Uc]])

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thsi removal of phases in incorrect. You need to multiply phases on both sides.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

As Nico said, this is incorrect. The correct method is outlined in https://journals.aps.org/pra/abstract/10.1103/PhysRevA.102.033519, paragraph between Eq.65 and 66. I also do it in my "phases" function in 'propagator.py'. Note that these are for SPDC where the idler mode is ^\dagger. For F.C. one needs to take the opposite sign for that mode.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I remember thinking about that and looking at it when looking for how to remove the phases. But then, I thought that since I have my waveguide start at z=0 and finish at z=L I only had phases on one side, since my z_0 would be 0. Is that incorrect or is there a profound reason that I should have z_0 = -L/2, z_1 = L/2, or is it just a convention? In any case, I can change it for sure, just so it remains consistent with what you have in propagator.py.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Right ok, yeah we should keep it all consistent.

Kc = np.diag(1j * dw * kc)
# Bottom left and upper right (non-diagonal) blocks

J = jnp.block([[Ka, -1j * jnp.conj(D) * pump], [-1j * D * pump.conj().T, Kc]])
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This way of generating the interaction matrix is incorrect. By combining everything into this 'D' you are omitting a \Delta \omega term which depends on the step size used. This term comes from discretizing the integral into a sum. (see Sec.4 of https://journals.aps.org/pra/abstract/10.1103/PhysRevA.102.033519, more precisely Eq. 59a for this off diagonal term). If you also look at my function "Hprop" in "propagator.py" I have a very convoluted way of expressing that \Delta \omega.

Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You will also need to divide the JCA by \Delta \omega. (See line 161 in propagator.py)

Returns:
array: pulse shape
"""
gaussian = np.exp(-(dw**2) / (2 * sig**2)) / np.sqrt(2 * np.pi * sig**2)
Copy link
Collaborator

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This seems to be L1 normalization, whereas we want to use L2 normalization. Denominator should be (pi*sig**2)^(1/4)

@jviatteau jviatteau requested review from MHoude2 and nquesada October 7, 2024 07:26
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

3 participants