BYUCamachoLab / emepy

https://emepy.rtfd.io
MIT License
33 stars 8 forks source link

The dimension of S-parameters matrix for a simple taper is not a 2x2 matrix? #32

Open SprBull opened 1 year ago

SprBull commented 1 year ago

Hi @hammy4815 , much thanks for developing this fantastic tool! I really enjoy diving into it.

But I did run into a bit of confusion when running an example of a simple taper calculation attached here(slightly modified from an example file). For instance, if you use a commercial EME tool(L*) to calculate a taper, it returns a S-matrix with a dimension 2x2 and S21 is directly related to the taper transmission. However, the S-matrix calculated from the example file gives a 4x4 S-matrix, which is confusing. Therefore, I have a few questions:

(1)Shouldn't a simple taper have only 1 input port and 1 output port, resulting in a 2x2 S-matrix(as in the case of some commercial EME tool)?

(2)In the 4x4 S-matrix, which S-parameter is directly related to the taper transmission?

Thanks in advance! taper1.zip

hammy4815 commented 1 year ago

Hi @SprBull ,

Looking at your example, it seems as though you are using 2 modes. The S-Matrix that EMEPy produces is basically the mode s-matrix, rather than just a physical port S-Matrix. This means the Input port is treated to have two inputs (Mode 1 and Mode 2) and the output port also has two ports (Mode 1 and Mode 2). Thus, there are 4 total "ports" for the system.

If you want to know the total power output, multiply your S-Matrix by your input vector (which includes the mode-decomposition in case you want to excite some superposition) and then sum your outputs at the corresponding ports.

For example, for your 2-mode taper, you have some 4x4 S-Matrix M. The input vector X takes the form [left0, left1, right0, right1]. If we send in a single mode source, our input vector takes the form [1, 0, 0, 0]. Y = MX; Y might look like [a, b, c, d]. To find the total output for all modes on the output, take $c^2 + d^2$ since those correspond to the output modes.

Note the total power output will only be perfectly accurate if you include a sufficient mode basis (enough modes). Hope this helps, let me know if you have any other questions!

SprBull commented 1 year ago

@hammy4815 , thanks for such a quick response!

Your answer has exposed my shallow understanding of S-parameters. For that I am deeply grateful. Based on your instructions, I further modified the taper program attached here.

(1)In the program I used 5 modes for each cross-section and I expected a 10x10 S-matrix. But to my confusion, I got a 4x4 S-matrix instead. Can you please tell me what I am doing wrong here?

(2)In this line: eme.propagate(left_coeffs=[1]), the manual says "left_coeffs (list) – A list of floats that represent the mode coefficients for the left side of the full geometry". I don't think I have grasped its meaning. Can you please elaborate on this?

Thanks again in advance! taper2.zip

hammy4815 commented 1 year ago

(1) One thing that is tricky about EME is that the user has to be careful with the mode selections in the solver. This is because once most of the guides modes are found, if the eigensolver keeps searching, it will sometimes find non-physical (spurious) modes based on the numerics of the boundary conditions. If these modes are included in the EME propagation, then the s-parameters will become incorrect.

To avoid this, EMEPy has a very simple method for detecting whether the mode basis contains any spurious modes. It does this basically by looking at the percentage of confinement in the core. If the mode is spread across the entire cladding, EMEPy labels it as spurious and throws it out of the EME propagation. As a result, the mode basis EMEPy uses might include less modes than those specified with num_modes.

Looking at your taper, if we plot the first 5 modes solved we can see:

Screen Shot 2022-09-21 at 3 03 57 PM Screen Shot 2022-09-21 at 3 03 31 PM Screen Shot 2022-09-21 at 3 04 16 PM Screen Shot 2022-09-21 at 3 04 27 PM Screen Shot 2022-09-21 at 3 04 40 PM Screen Shot 2022-09-21 at 3 04 52 PM

It's quite clear that aside from the first two modes, these field profiles don't look physical. EMEPy is finding this as well with it's power calculation and throwing these modes out of the simulation reducing the number of modes here to 2.

(2) These coefficient arrays represent the mode coefficients for in the input to the scattering matrix. left_coeffs specifically represents the inputs to the propagating modes fed into the left port of the device. Thus, [1] represents exciting and sending in the first mode solved for (usually fundamental TE). [0, 1] represents exciting the second mode solved for (either the fundamental TM or TE1 mode depending on the boundaries, symmetries, etc.). If we instead use right_coeffs, then they will be excited from the right port of the device propagation towards the left.

Does that answer your questions?

694856790 commented 1 year ago

(1) One thing that is tricky about EME is that the user has to be careful with the mode selections in the solver. This is because once most of the guides modes are found, if the eigensolver keeps searching, it will sometimes find non-physical (spurious) modes based on the numerics of the boundary conditions. If these modes are included in the EME propagation, then the s-parameters will become incorrect.

To avoid this, EMEPy has a very simple method for detecting whether the mode basis contains any spurious modes. It does this basically by looking at the percentage of confinement in the core. If the mode is spread across the entire cladding, EMEPy labels it as spurious and throws it out of the EME propagation. As a result, the mode basis EMEPy uses might include less modes than those specified with num_modes.

Looking at your taper, if we plot the first 5 modes solved we can see:

Screen Shot 2022-09-21 at 3 03 57 PM Screen Shot 2022-09-21 at 3 03 31 PM Screen Shot 2022-09-21 at 3 04 16 PM Screen Shot 2022-09-21 at 3 04 27 PM Screen Shot 2022-09-21 at 3 04 40 PM Screen Shot 2022-09-21 at 3 04 52 PM

It's quite clear that aside from the first two modes, these field profiles don't look physical. EMEPy is finding this as well with it's power calculation and throwing these modes out of the simulation reducing the number of modes here to 2.

(2) These coefficient arrays represent the mode coefficients for in the input to the scattering matrix. left_coeffs specifically represents the inputs to the propagating modes fed into the left port of the device. Thus, [1] represents exciting and sending in the first mode solved for (usually fundamental TE). [0, 1] represents exciting the second mode solved for (either the fundamental TM or TE1 mode depending on the boundaries, symmetries, etc.). If we instead use right_coeffs, then they will be excited from the right port of the device propagation towards the left.

Does that answer your questions?

Hi, @hammy4815. I have a question that I want to discuss with you. Is there a certain standard for the percentage standard for judging whether it is a spurious mode in EMEpy? Can we customize this standard, can you help me find where this standard can be modified in eme.py. Also, I was using Lumerical MODE EME and found that they don't seem to exclude spurious modes, but let the user use “Mode convergence sweep” to evaluate the amount of modes selected (propagate all modes selected). The function of excluding spurious modes may help to improve the accuracy of calculation results based on the original FDE solver, but on the basis of Tidy3d's FDE, it will make the results inaccurate, because Tidy3d's FDE is highly consistent with Lumerical's results, Excluding some modes may cause errors in the results. Thank you very much!

hammy4815 commented 1 year ago

Hi @694856790 ,

On line 133 of mode.py, there is a method:


  def check_spurious(self, threshold_power: float = 0.05, threshold_neff: float = 0.9) -> bool:
      """Takes in a mode and determine whether the mode is likely spurious based on the ratio of confined to not confined power

      Parameters
      ----------
      threshold_power : float
          threshold of power percentage of Pz in the core to total
      threshold_neff : float
          threshold of real to abs neff percentage

      Returns
      -------
      boolean
          True if likely spurious
      """

      power_bool = self.get_confined_power() < threshold_power
      neff_bool = (np.real(self.neff) / np.abs(self.neff)) < threshold_neff
      return power_bool or neff_bool

Basically, power_bool calculates the ratio of the integrated Poynting vector in just the core, to the integrated Poynting vector of the entire profile. A high power_bool implies high confinement, whereas a low value implies low confinement. You can see here the default value is 0.05. This means if there is less than 5% confinement, the mode is taken to be spurious. You can change the default value of the threshold_power argument in the method from 0.05 to a value of your choice to change this.

In the example provided, it seems that the electromagnetic python solver (EMpy) is being used rather than Tidy3D, but both solvers are consistent. Mode convergence sweep functionality would be useful, it would just require implementation.

I'm slightly confused by your last statement, however. Are you suggesting if using Tidy3D, we should be including the spurious modes in the EME propagation?

Thanks, Ian

694856790 commented 1 year ago

Hi @694856790 ,

On line 133 of mode.py, there is a method:

  def check_spurious(self, threshold_power: float = 0.05, threshold_neff: float = 0.9) -> bool:
      """Takes in a mode and determine whether the mode is likely spurious based on the ratio of confined to not confined power

      Parameters
      ----------
      threshold_power : float
          threshold of power percentage of Pz in the core to total
      threshold_neff : float
          threshold of real to abs neff percentage

      Returns
      -------
      boolean
          True if likely spurious
      """

      power_bool = self.get_confined_power() < threshold_power
      neff_bool = (np.real(self.neff) / np.abs(self.neff)) < threshold_neff
      return power_bool or neff_bool

Basically, power_bool calculates the ratio of the integrated Poynting vector in just the core, to the integrated Poynting vector of the entire profile. A high power_bool implies high confinement, whereas a low value implies low confinement. You can see here the default value is 0.05. This means if there is less than 5% confinement, the mode is taken to be spurious. You can change the default value of the threshold_power argument in the method from 0.05 to a value of your choice to change this.

In the example provided, it seems that the electromagnetic python solver (EMpy) is being used rather than Tidy3D, but both solvers are consistent. Mode convergence sweep functionality would be useful, it would just require implementation.

I'm slightly confused by your last statement, however. Are you suggesting if using Tidy3D, we should be including the spurious modes in the EME propagation?

Thanks, Ian

@hammy4815 Thank you! I'll test the accuracy after adjusting the limit. I state that if you use Tidy3d, it is better to propagate all the modes including spurious modes, which may be better. Because through my test, I found that the higher order modes calculated by EMpy is obviously different from that calculated by Lumerical and Comsol, while the mode calculated by Tidy3d is highly consistent with that calculated by Lumerical and Comsol.

SprBull commented 1 year ago

@hammy4815 Understood sir! Thanks again for your detailed explanation! I have learned a great deal.

In the mean time, I have been doing my homework by digging deeper into EME fundamentals. According to literatures, EME is capable of not only straight waveguides, but also bend waveguides, even grating couplers. I know that EMEpy allows the cascade of layers from left to right, thus enabling the simulation of a straight waveguide. What about a bend waveguide then? Does EMEpy allow an additional layer to be attached to a previous layer with an angle? And what about grating couplers(waveguide input and vertical output)? Does EMEpy allow a new layer to be added on top of previous layers? I am very curious.

hammy4815 commented 1 year ago

Thank you! I'll test the accuracy after adjusting the limit. I state that if you use Tidy3d, it is better to propagate all the modes including spurious modes, which may be better. Because through my test, I found that the higher order modes calculated by EMpy is obviously different from that calculated by Lumerical and Comsol, while the mode calculated by Tidy3d is highly consistent with that calculated by Lumerical and Comsol.

@694856790 , This is a good observation. In theory, EMpy and Tidy3D should be finding similar higher order modes. I am not sure why they disagree, but we will look more into this. Thanks!

hammy4815 commented 1 year ago

@hammy4815 Understood sir! Thanks again for your detailed explanation! I have learned a great deal.

In the mean time, I have been doing my homework by digging deeper into EME fundamentals. According to literatures, EME is capable of not only straight waveguides, but also bend waveguides, even grating couplers. I know that EMEpy allows the cascade of layers from left to right, thus enabling the simulation of a straight waveguide. What about a bend waveguide then? Does EMEpy allow an additional layer to be attached to a previous layer with an angle? And what about grating couplers(waveguide input and vertical output)? Does EMEpy allow a new layer to be added on top of previous layers? I am very curious.

Hi @SprBull ,

With the introduction of Tidy3D, EMEPy now has the framework to solve for mode profiles of cross sections with bends. However, introducing bends into the EME algorithm is more complex and will take a bit more work. Ideally we would like to get there eventually, but there is not a plan for the immediate future to implement bends. As for vertical devices, at the moment if you want any vertical coupling, you need to include the device in the same profile as the device below, which might be infeasible for many devices. We have not previously considered adding functionality for vertical layers, but that is something to keep in mind. Thanks for the suggestion!

694856790 commented 1 year ago

@hammy4815 Understood sir! Thanks again for your detailed explanation! I have learned a great deal. In the mean time, I have been doing my homework by digging deeper into EME fundamentals. According to literatures, EME is capable of not only straight waveguides, but also bend waveguides, even grating couplers. I know that EMEpy allows the cascade of layers from left to right, thus enabling the simulation of a straight waveguide. What about a bend waveguide then? Does EMEpy allow an additional layer to be attached to a previous layer with an angle? And what about grating couplers(waveguide input and vertical output)? Does EMEpy allow a new layer to be added on top of previous layers? I am very curious.

Hi @SprBull ,

With the introduction of Tidy3D, EMEPy now has the framework to solve for mode profiles of cross sections with bends. However, introducing bends into the EME algorithm is more complex and will take a bit more work. Ideally we would like to get there eventually, but there is not a plan for the immediate future to implement bends. As for vertical devices, at the moment if you want any vertical coupling, you need to include the device in the same profile as the device below, which might be infeasible for many devices. We have not previously considered adding functionality for vertical layers, but that is something to keep in mind. Thanks for the suggestion!

@hammy4815 Hi, do the current EME algorithm support the calculation of Sbend waveguide? I want to use EME to calculate some waveguide structures with curve shape.

hammy4815 commented 1 year ago

Hi @694856790,

The EME framework is not set up for any bends. The direction of propagation is fixed in z for simplicity, but eventually we might expand and allow propagation in multiple dimensions with bends.

You can, however, still perform overlaps between straight waveguide modes and bend modes for a simpler analysis with the Tidy3D mode solver. You might be able to analytically simplify the EME process for this particular problem doing so. If you wanted a field profile you might be able to analytically propagate your meshed field profiles along some curved segment if the bend is adiabatic enough. But this would still be a little tricky working in euclidean coordinates.

SprBull commented 1 year ago

@hammy4815 Mr Hammond, it took me a while to fully digest your instructions. I thought I completely got it but my attempt says otherwise.

In the attachment, a simple taper is simulated. In each segment, 10 modes are calculated. With spurious modes dumped, a (8*8) S-matrix is created. Surprisingly, the absolute values of many elements are greater than 1. Then transmission of the 1st mode is calculated but the result shows a very large number ~2500.

I checked your instructions again but could not find errors in the program. What have I done wrong this time? Thanks again in advance!

taper3.zip

694856790 commented 1 year ago

Hi @694856790,

The EME framework is not set up for any bends. The direction of propagation is fixed in z for simplicity, but eventually we might expand and allow propagation in multiple dimensions with bends.

You can, however, still perform overlaps between straight waveguide modes and bend modes for a simpler analysis with the Tidy3D mode solver. You might be able to analytically simplify the EME process for this particular problem doing so. If you wanted a field profile you might be able to analytically propagate your meshed field profiles along some curved segment if the bend is adiabatic enough. But this would still be a little tricky working in euclidean coordinates.

@hammy4815 Hi, I feel that I may not understand what you mean. Is it possible to solve Sbend by dividing enough units?