Discrepancy between experimental data and LWE simulations

[timk110]

Hi Nick,

I have been using LWE to model second harmonic generation of a transform limited 1030 nm, 230 fs laser pulse. I have been using BIBO (1 mm thick)and BBO (0.5 mm thick) as my nonlinear crystals. When using LWE to generate plots of how the 1030 nm-515 nm conversion efficiency varies as a function of 1030 nm pulse energy I seem to be getting some very odd results.

I made sure to follow your youtube tutorial on setting the spatial and temporal grid, as well as doing a convergence test in dz (using the 3D radial symmetry propagation mode). I also made sure to download the latest version of LWE.

For the BIBO, when I scan the 1030 nm pulse energy from 1-90 microjoules I get very odd oscillatory behaviour in the conversion efficiency which is shown in the plot below. Are you aware of any papers in the literature explaining this phenomenon ? If not is it possible that the simulation is not functioning correctly?

I have attached the txt file summarising my simulation parameters (the raw data was too large too upload).
2hg bibo 1030nm pulse energy scan 361 micrometer bw 1-90 microjoule PE.txt

I also tested BBO with the same variation in fundamental pulse energy (1-90 microjoule) and instead I get a very different covnersion efficiency curve, which is much closer to what I expect . A monotonic increase until a plateau is reached. However here the maximum conversion efficiency reaches nearly 80%, which is very different from what I can find in the literature and measure in the lab (60% seems to be around the maximum, and this is with a experimental system with M2 <1.2 and near transform limited pulses)). What could explain the difference?

Here are the simulation parameters again
2hg bbo 1030nm pulse energy scan 0.5mm 90 deg phi.txt

Any help would be greatly appreciated!

Cheers,

Tim

[NickKarpowicz]

Hi Tim,
thanks for pointing this out!

I think you might have found a physics bug rather than a code bug (which is more interesting). I would indeed expect a more monotonic increase of efficiency with intensity for SHG (for parametric amplifiers you can indeed see heavy oscillations back and forth between signal/idler and pump, but SHG doesn't exhibit it as strongly).

My first suspicion was something with the nonlinear coefficients being off. The ones in the LWE database were taken from M. Ghotbi and M. Ebrahim-Zadeh, Optics Express 12, 6002-6019 (2004), and look like this:
2.53 2.93 -1.93 1.63 0 0
0 0 0 0 1.67 3.48
0 0 0 0 -1.58 1.67

Now, first thing to look at if you think something is off with nonlinear coefficients is Boyd. On p. 48 (3rd ed), you can see the symmetry maps of various crystal classes:

BiBO is space group C2, so it should be class 2 there. Big black dots mean nonzero values, and connecting lines mean they should have the same value under Kleinmann symmetry. The tensor up there doesn't obey this symmetry at all!

So, I went back to the paper. They mention that they take their coefficients from an earlier paper, H. Hellwig, J. Liebertz, and L. Bohaty, J.Appl. Phys. 88, 240 (2000), and apply a rotation operation to them. I think this operation is not necessary (but should think about it more).

If we go back to the original paper, and write out their coefficients in contracted form, we get
0 0 0 2.4 0 2.8
2.3 2.53 -1.3 0 2.3 0
0 0 0 -0.9 0 2.4

Now comparing to Boyd's diagrams, we see the right symmetry! It's a little bit off, so just to be fully sure, I'll average the equivalent ones, coming up with this:
0 0 0 2.36 0 2.55
2.55 2.53 -1.1 0 2.36 0
0 0 0 -1.1 0 2.36

Now, let's compare the resulting curves, using your simulation parameters, but just changing the nonlinear tensor in the database:

I think this is the culprit!

Could you try putting these values into your CrystalDatabase.txt and see if the results seem more consistent with what you'd expect? Especially if you've compared with literature, I'd be happy to see it!

Now, for the next part, why is BBO efficiency so high? I'm actually not too sure; I know I have seen lab examples of ~80% conversion, but I don't know an example in BBO (here's a look in LBO, along with some back-conversion oscillations).

I'll have another look at the literature later, and see how consistent the values of the chi(2) of BBO really are - it could be that the stated values are simply too high...

In any case, thanks for the feedback!

[timk110]

Hi Nick,

Thank you for the detailed explanation! I would never have thought to use the crystal classes to check whether the experimentally reported nonlinear coefficients are physically valid. This will be a good check for me to use in the future if I need to add any other crystals to LWE. Its reassuring that the conversion efficiency dependence for BIBO looks a lot more realistic now.

I did find this paper that looked at SHG of nanojoule level fs pulses in BIBO (Galletti M, Pires H, Hariton V, et al. High efficiency second harmonic generation of nanojoule-level femtosecond pulses in the visible based on BiBO. High Power Laser Science and Engineering. 2019;7:e11. doi:10.1017/hpl.2018.72). Their description of their setup is detailed enough that I can reproduce it in LWE. I will let you know how they compare.

I also have some experimental data from my own lab for SHG in BIBO (this is actually what prompted me to use LWE to simualte such pulses) that I will compare to the LWE simulation reuslts with the updated nonlinear tensor.

Thank you also for looking at the BBO literature. I will also have a look to see if I can find anything obvious.

[NickKarpowicz]

Hi Tim,
I have to say I'm still not sure if something else is wrong - If I use the non-rotated chi2 tensor for BiBO, I get optimum angles for OPCPA that don't agree with what people do in the literature, but the original one does. I'm looking around for other things that could be off at the moment...
Best wishes,
Nick

[NickKarpowicz]

Hi again,

So, I think I found a rather embarassing bug on my part in addition to the nonlinear tensor problem. When I refactored the nonlinear polarization calculation, I removed the factor of two on the cross-terms when calculating the elements of the chi(2) tensor in UPPE mode, when I added the FDTD propagation! I think this explains the unexpectedly high BBO conversion efficiency and brings the calculation with BiBO in line with the measurements of Galletti et al.

Here are some plots...

First, plotting the patched and unpatched versions of LWE (patched=after I restored the missing factor 2 last night), with the two versions of the BiBO chi(2) tensor we discussed earlier. I also extracted the data from Galletti at 1000 nm to compare.

This is the LWE input file I used taking the spot size and pulse width from their paper
TestSHG_Galletti_mod_x2.txt

The agreement with both patches applied (LWE and the tensor from the older paper) is promising at least.

If I apply the same patch to your BBO calculation from earlier, you can see that the efficiency is cut in ~half:

You're using Windows, right? I attached the patched version of the LWE executable here:
LightwaveExplorer.zip

After a bit more testing I would actually do a patch update soon, it seems to be an important bug you caught!

[timk110]

Hi Nick,

If this rogue factor of 2 does turn out to be the culprit then I am happy to have at least contributed a little towards this project! Given how useful this code has been to my own work its nice to give something back.

The plot comparing the various LWE versions to Galetti et al. look very promising. I will run the patched version of the LWE execturable and see how it compares to my BIBO SHG data from the lab. I will also see if i can find some published work on SHG in BBO that is detailed enough so that I can reproduce it in LWE and see how it compares.

On a seperate note I did some more tests on SHG in BIBO yesterday (using the old LWE executable) and found some odd behaviour when just varying the crystal angle Phi. From SNLO the phase matching angles for SHG at 1030 nm is theta=166.6 deg, and phi =90 deg, for doing eeo SHG in the YZ plane. However if I do a parameter scan from 85-95 deg for phi I get the following (this is using the original Chi2 tensor).

bibo 1mm shg phi 85-95 old chi 2 unpatched.txt

I am a bit confused why the conversion efficiency should peak for a phi that is slightly off the 90 degrees that SNLO states.

I then used your patched LWE exectuable and your altered Chi2 tensor and repeated the parameter to scan to get the following
bibo 1mm shg phi 85-95 new chi 2 PATCHED.txt

From my understanding of the crystal geometry a phi of 90 means the fundamental wavevector is in the YZ plane. Does this just mean that using BIBO with the fundamental propagation vector not quite being aligned to the YZ plane gives us a better conversion efficiency?

Thanks again for the quick turn around!

Cheers,

Tim

[NickKarpowicz]

Hi Tim,
I think what's going on there is that since BiBO is biaxial, changing phi is also altering the refractive index, plus for perfect phase matching it wants theta ~166.78. If you put in that at theta, the two peaks will merge, but still be skewed.

The refractive index affects the strength of the nonlinearity in LWE - both because it shows up as in front of the nonlinear polarization in the nonlinear wave equation, and because Miller's rule is used to approximate the dispersion of the nonlinear coefficients. My guess is that phi=90 will be the maximum in terms of pure nonlinear coefficient, but going a bit off from that angle changes the refractive index enough to shift the overall efficiency maximum. If that's the case, I'd expect this effect not to happen in BBO since it's uniaxial and phi doesn't change refractive index - is that so?

Best,
Nick

[timk110]

Hi Nick,

That makes sense. As you predicted for the BIBO for theta 166.78 deg there is a single peak in the conversion efficiency as phi varies, but its slightly shifted from 90 deg. Everything was done with the patch and the new Chi2 tensor

vary phi bibo.txt

For BBO, theta =23.376 deg
varying phi bbo.txt

Your argument with regards to BBO being uniaxial makes sense, rotations around the Z axis should not affect the refractive index, however there does seem to be a dependence.
Cheers,

Tim

[NickKarpowicz]

Hi Tim,
ah sorry, I meant that there will be no dependence of the refractive index on phi, there will still be a change in d_eff, so what you show makes sense.

[NickKarpowicz]

also, another update on the nonlinear tensor part. I thought a bit more about what they were doing with the rotations. In the original paper by Hellwig et al., they use the crystal coordinate system, which is why everything is consistent with the symmetry in Boyd. What the 2004 paper did was to relabel the axes, putting things in line with the conventions where the first axis has the smallest index, and they go up in order. The refractive index as was entered in LWE was also within this rotation, taken from a later paper. So, I think it would be best to keep with this convention, since if anyone orders a crystal based on their simulation, it's also what the manufacturers would expect in terms of angles...

I still want to look into this more, I'm now sure that if one uses the tensor I "corrected" above, the refractive index values for the different axes have to be changed as well, so that better agreement with Galletti's measurements must have been coincidental.

Sorry for the slow convergence on this! Crystal symmetry is confusing sometimes!

[NickKarpowicz]

I am starting to think that in the original paper by Hellwig, they /already/ put the tensor in the optical notation. They were certainly aware that people in nonlinear optics were expecting this, since they did put their definitions of refractive index in two different ways, the first one with n_i, then with n_alpha, n_beta, and n_gamma, which they used to define the phase matching angles for nonlinear optics.

I'm starting to think this since I took data from another study, Ghotbi et al., Optics Letter 29, 2530 (2004), and made a similar plot to the one above:

Again, if I use the non-rotated Hellwig tensor, the agreement is good, and if I use the rotated one, it's both bad and falling below the experiment (there are usually more reasons why real life should be worse than simulations than vice versa...)

Unfortunately it's hard to tell from the text which system of coordinates is being used where. Maybe we can figure it out from the data available however; if there are any other cases you see where they specify (or give values that can be translated to) pulse duration, pulse energy, and spot size, and have a reasonable pulse length (<500 fs), I'd be grateful!

[timk110]

So just to clarify with regards to the coordinate systems. The "crystal coordinate system" just uses the normal crystal lattice vectors a,b,c and the angles between them alpha,beta,gamma.

The "optical notation" you mention is the coordinate system where the dielectric tensor is diagonalized (as defined on Crystal nonlinear optics:with SNLO examples pg 23)
snlo.pdf

I am not the best with the various coordiante systems used (all the different notations make my head swim) so apologies if I misunderstood something.

So with regards to the discrpeancy between the Helwig paper and the Ghotbi paper, the Helwig paper was already in the optical notation, and then Ghotbi et al. performed essentially a uneccecary rotation because they mistakenly believed that Helwig had their quantities in the crystal notation?

For further works in the literature on second harmonic generation in BIBO I came across these. They qoute all the neccecary parameters to reproduc in LWE.
Type-I frequency-doubling characteristics of high-power,.pdf

Frequency-Doubling_of_Femtosecond_Pulses_in_Thick_Nonlinear_Crystals_With_Different_Temporal_and_Spatial_Walk-Off_Parameters.pdf

(In this review section 5.2 is the relevant one)
Laser Photonics Reviews - 2010 - Petrov - Femtosecond nonlinear frequency conversion based on BiB3O6 (1).pdf

Hopefully these help.

Cheers,

Tim

[NickKarpowicz]

Exactly - I think the tensor by Hellwig was already in the optical basis (where the dielectric tensor is diagonalized like you said). After some more digging through the literature, I'm actually quite sure of it now. In an earlier publication from 1998, they in fact already published the same tensor, and importantly, they have their refractive index axes laid out in the "optical" order:

In their followup paper, the one cited by the Ghotbi paper, they adopt the rotated system, but only for the refractive index, the chi(2) tensor stays exactly the same.

I did another simulation of a different process, this time an OPCPA, in Kurihara et al., Optics Express 31, 11649 (2023). The non-rotated tensor (Hellwig) predicts an optimized conversion of ~32 microJoules, whereas the rotated one predicts ~170 microJoules - the authors of the paper measure 32-35 microJoules at the output. I think there's no way they would leave behind that much energy if the effective nonlinearity was really as high as the rotated tensor would suggest, and again the original tensor shows better agreement with measurements.

I think I've convinced myself at this point that the tensor as stated in the original paper by Hellwig et al is the correct one to use; I'll do a new release of LWE with the correction later today...

[timk110]

Great! It's nice that there's a satisfying resolution to this. Thanks for your patience.

[NickKarpowicz]

Hi again!
One more update to this... there actually wasn't a factor 2 error in the old code, I was just looking for the 2 in the wrong place when I was worried that something went wrong! I'm fixing it up now and will post a fixed version soon...

[timk110]

Ah I see! Thanks for the info.